空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

Turbulence and PIML-Based Modeling

Turbulence

AI4S

TurbulenceAI4S

JiaweiMiao

发布于 2023-10-15

推荐镜像 :Third-party software:ai4s-cup-0.1

推荐机型 :c2_m4_cpu

数据集

NB-Nov-23(v1)

Turbulence and PIML-Based Modeling

代码

文本

[1]

from IPython.display import Image

代码

文本

Introduce about Turbulence

代码

文本

Turbulence Simulation Intro

代码

文本

Overview

本Notebook部分内容搬运自github/xiaoh/turbulence-modeling-PIML

This code implements the PIML regresson procedure described in Wang et al. 2017, with the case of flow over periodic hills as example.

J.-X. Wang, J.-L. Wu, and H. Xiao. Physics informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data. Physical Review Fluids. 2(3), 034603, 1-22, 2017. https://doi.org/10.1103/PhysRevFluids.2.034603[DOI:10.1103/PhysRevFluids.2.034603]

代码

文本

Algorithm of PIML-Based Turbulence Modeling

The overall procedure can be summarized as follows:

Perform baseline RANS simulations on both the training flows and the test flow.
Compute the input feature field $q (x)$ based on the local RANS flow variables.
Compute the discrepancies field $Δ τ (x)$ in the RANS-modeled Reynolds stresses for the training flows based on the high-fidelity data.
Construct regression functions $f : q \mapsto Δ τ$ for the discrepancies based on the training data prepared in Step 3, using machine learning algorithms.
Compute the Reynolds stress discrepancies for the test flow by querying the regression functions. The Reynolds stresses can subsequently be obtained by correcting the baseline RANS predictions with the evaluated discrepancies.
Propagate the corrected Reynolds stresses to the mean velocity field by solving the RANS equations with the corrected Reynolds stress field.

This code only performs Step 4. (see the green-shaded box below in the flow chart). The training data prepared in Steps 1-3 are saved in database folder.

代码

文本

[2]

Image(filename='figs/PIML-algorithm.png')

代码

文本

Machine learning algorithms

The procedure implented here consists of three parts:

load training and test data
construct regression function $Δ τ (q)$ (detailed below)
plot the anisotropy parameters $ξ$ and $η$ (componebnts of $Δ τ$ ) and compare with ground truth (DNS)

We used two algorithms to build the regression function:

Random Forests (based on scikit-learn). This is what was used in Wang et al.
Neural networks (based on Tensorflow)

Both algorithms yielded similar results, but the former is cheaper computationally.

代码

文本

The input features consist of 12 variables (see Table 1 below and also Wang et al.)

代码

文本

[3]

Image(filename='/bohr/nbnov23-tuvx/v1/figs/features.png')

代码

文本

[5]

!pip install scikit-learn tensorflow

Looking in indexes: https://pypi.org/simple, https://pypi.ngc.nvidia.com
Requirement already satisfied: scikit-learn in /opt/conda/lib/python3.8/site-packages (0.24.2)
Requirement already satisfied: tensorflow in /opt/conda/lib/python3.8/site-packages (2.13.0)
Requirement already satisfied: joblib>=0.11 in /opt/conda/lib/python3.8/site-packages (from scikit-learn) (1.1.0)
Requirement already satisfied: numpy>=1.13.3 in /opt/conda/lib/python3.8/site-packages (from scikit-learn) (1.22.4)
Requirement already satisfied: scipy>=0.19.1 in /opt/conda/lib/python3.8/site-packages (from scikit-learn) (1.6.3)
Requirement already satisfied: threadpoolctl>=2.0.0 in /opt/conda/lib/python3.8/site-packages (from scikit-learn) (3.1.0)
Requirement already satisfied: google-pasta>=0.1.1 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (0.2.0)
Requirement already satisfied: wrapt>=1.11.0 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (1.15.0)
Requirement already satisfied: tensorflow-io-gcs-filesystem>=0.23.1 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (0.34.0)
Requirement already satisfied: keras<2.14,>=2.13.1 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (2.13.1)
Requirement already satisfied: tensorflow-estimator<2.14,>=2.13.0 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (2.13.0)
Requirement already satisfied: absl-py>=1.0.0 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (1.0.0)
Requirement already satisfied: opt-einsum>=2.3.2 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (3.3.0)
Requirement already satisfied: six>=1.12.0 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (1.16.0)
Requirement already satisfied: grpcio<2.0,>=1.24.3 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (1.58.0)
Requirement already satisfied: packaging in /opt/conda/lib/python3.8/site-packages (from tensorflow) (21.3)
Requirement already satisfied: typing-extensions<4.6.0,>=3.6.6 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (4.5.0)
Requirement already satisfied: flatbuffers>=23.1.21 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (23.5.26)
Requirement already satisfied: h5py>=2.9.0 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (3.9.0)
Requirement already satisfied: gast<=0.4.0,>=0.2.1 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (0.4.0)
Requirement already satisfied: protobuf!=4.21.0,!=4.21.1,!=4.21.2,!=4.21.3,!=4.21.4,!=4.21.5,<5.0.0dev,>=3.20.3 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (4.24.3)
Requirement already satisfied: libclang>=13.0.0 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (16.0.6)
Requirement already satisfied: tensorboard<2.14,>=2.13 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (2.13.0)
Requirement already satisfied: setuptools in /opt/conda/lib/python3.8/site-packages (from tensorflow) (59.5.0)
Requirement already satisfied: astunparse>=1.6.0 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (1.6.3)
Requirement already satisfied: termcolor>=1.1.0 in /opt/conda/lib/python3.8/site-packages (from tensorflow) (2.3.0)
Requirement already satisfied: wheel<1.0,>=0.23.0 in /opt/conda/lib/python3.8/site-packages (from astunparse>=1.6.0->tensorflow) (0.37.1)
Requirement already satisfied: werkzeug>=1.0.1 in /opt/conda/lib/python3.8/site-packages (from tensorboard<2.14,>=2.13->tensorflow) (2.1.1)
Requirement already satisfied: markdown>=2.6.8 in /opt/conda/lib/python3.8/site-packages (from tensorboard<2.14,>=2.13->tensorflow) (3.3.6)
Requirement already satisfied: google-auth<3,>=1.6.3 in /opt/conda/lib/python3.8/site-packages (from tensorboard<2.14,>=2.13->tensorflow) (2.23.0)
Requirement already satisfied: google-auth-oauthlib<1.1,>=0.5 in /opt/conda/lib/python3.8/site-packages (from tensorboard<2.14,>=2.13->tensorflow) (1.0.0)
Requirement already satisfied: tensorboard-data-server<0.8.0,>=0.7.0 in /opt/conda/lib/python3.8/site-packages (from tensorboard<2.14,>=2.13->tensorflow) (0.7.1)
Requirement already satisfied: requests<3,>=2.21.0 in /opt/conda/lib/python3.8/site-packages (from tensorboard<2.14,>=2.13->tensorflow) (2.31.0)
Requirement already satisfied: pyasn1-modules>=0.2.1 in /opt/conda/lib/python3.8/site-packages (from google-auth<3,>=1.6.3->tensorboard<2.14,>=2.13->tensorflow) (0.2.8)
Requirement already satisfied: urllib3<2.0 in /opt/conda/lib/python3.8/site-packages (from google-auth<3,>=1.6.3->tensorboard<2.14,>=2.13->tensorflow) (1.26.16)
Requirement already satisfied: cachetools<6.0,>=2.0.0 in /opt/conda/lib/python3.8/site-packages (from google-auth<3,>=1.6.3->tensorboard<2.14,>=2.13->tensorflow) (5.0.0)
Requirement already satisfied: rsa<5,>=3.1.4 in /opt/conda/lib/python3.8/site-packages (from google-auth<3,>=1.6.3->tensorboard<2.14,>=2.13->tensorflow) (4.8)
Requirement already satisfied: requests-oauthlib>=0.7.0 in /opt/conda/lib/python3.8/site-packages (from google-auth-oauthlib<1.1,>=0.5->tensorboard<2.14,>=2.13->tensorflow) (1.3.1)
Requirement already satisfied: importlib-metadata>=4.4 in /opt/conda/lib/python3.8/site-packages (from markdown>=2.6.8->tensorboard<2.14,>=2.13->tensorflow) (6.8.0)
Requirement already satisfied: zipp>=0.5 in /opt/conda/lib/python3.8/site-packages (from importlib-metadata>=4.4->markdown>=2.6.8->tensorboard<2.14,>=2.13->tensorflow) (3.16.2)
Requirement already satisfied: pyasn1<0.5.0,>=0.4.6 in /opt/conda/lib/python3.8/site-packages (from pyasn1-modules>=0.2.1->google-auth<3,>=1.6.3->tensorboard<2.14,>=2.13->tensorflow) (0.4.8)
Requirement already satisfied: certifi>=2017.4.17 in /opt/conda/lib/python3.8/site-packages (from requests<3,>=2.21.0->tensorboard<2.14,>=2.13->tensorflow) (2023.7.22)
Requirement already satisfied: charset-normalizer<4,>=2 in /opt/conda/lib/python3.8/site-packages (from requests<3,>=2.21.0->tensorboard<2.14,>=2.13->tensorflow) (2.0.12)
Requirement already satisfied: idna<4,>=2.5 in /opt/conda/lib/python3.8/site-packages (from requests<3,>=2.21.0->tensorboard<2.14,>=2.13->tensorflow) (3.3)
Requirement already satisfied: oauthlib>=3.0.0 in /opt/conda/lib/python3.8/site-packages (from requests-oauthlib>=0.7.0->google-auth-oauthlib<1.1,>=0.5->tensorboard<2.14,>=2.13->tensorflow) (3.2.0)
Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in /opt/conda/lib/python3.8/site-packages (from packaging->tensorflow) (3.0.8)
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv

代码

文本

[7]

%matplotlib inline

## Import system modules

# sci computing

import numpy as np

# sklearn importing

from sklearn.ensemble import RandomForestRegressor

# plotting

import matplotlib.pyplot as plt # for plotting

#import matplotlib as mp

# keras importing

from keras.models import Sequential

from keras.layers import Dense

import time

代码

文本

[8]

def loadTrainingData(caseName, ReNum):

trainFeaturesFile = '/bohr/nbnov23-tuvx/v1//database/' + caseName + '/markers/' + ReNum + '/markerFile'

trainResponsesFile = '/bohr/nbnov23-tuvx/v1//database/' + caseName + '/deltaFields/' + ReNum + '/deltaField'

trainFeatures = np.loadtxt(trainFeaturesFile)

trainResponses = np.loadtxt(trainResponsesFile)

return trainFeatures, trainResponses

代码

文本

[9]

def loadTestData(caseName, ReNum):

testFeaturesFile = '/bohr/nbnov23-tuvx/v1//database/' + caseName + '/markers/' + ReNum + '/markerFile'

testResponsesFile = '/bohr/nbnov23-tuvx/v1//database/' + caseName + '/deltaFields/' + ReNum + '/deltaField'

testFeatures = np.loadtxt(testFeaturesFile)

testResponses = np.loadtxt(testResponsesFile)

return testFeatures, testResponses

代码

文本

[10]

def randomForest(trainFeatures, trainResponses, testFeatures, maxFeatures = 'log2', nTree=100):

## Settings of random forests regressor

regModel = RandomForestRegressor(n_estimators=nTree, max_features=maxFeatures)

## Train the random forests regressor

regModel.fit(trainFeatures, trainResponses)

## Prediction

testResponsesPred = regModel.predict(testFeatures)

return testResponsesPred

代码

文本

[11]

def keras_nn(trainFeatures, trainResponses, testFeatures, Nepochs = 100):

'''

This function is to construct neural network based on the training data and predict the

response given test data.

Two hidded layers are used, the number of neurals are 64 and 32, respectively.

'''

model = Sequential()

# The first hidder layer of NN

model.add(Dense(64, input_dim=trainFeatures.shape[1], activation='relu'))

# The second hidder layer of NN

model.add(Dense(32, activation='relu'))

model.add(Dense(2, activation='tanh'))

model.compile(loss='mean_squared_error', optimizer='adam')

# Training

model.fit(trainFeatures, trainResponses, epochs=Nepochs, batch_size=200, verbose=0)

# Prediction

testResponsesPred = model.predict(testFeatures)

return testResponsesPred

代码

文本

[12]

def plotXiEta(XiEta_RANS, testResponses, testResponsesPred, name, symbol='r^'):

# Reconstruct Barycentric coordinates

XiEta_DNS = XiEta_RANS + testResponses

XiEta_ML = XiEta_RANS + testResponsesPred

# Plot Reynolds stress anisotropy in Barycentric triangle

interval = 2

pointsNum = int(XiEta_RANS.shape[0])

plt.figure()

plt.plot([0,1,0.5,0.5,0],[0,0,3**0.5/2.0,3**0.5/2.0,0],'g-')

p1, = plt.plot(XiEta_RANS[:pointsNum:interval,0],XiEta_RANS[:pointsNum:interval,1],

'bo', markerfacecolor='none', markeredgecolor='b',

markeredgewidth=2, markersize=10)

p2, = plt.plot(XiEta_DNS[:pointsNum:interval,0],XiEta_DNS[:pointsNum:interval,1],

'ks', markerfacecolor='none', markeredgecolor='k',

markeredgewidth=2, markersize=10)

p3, = plt.plot(XiEta_ML[:pointsNum:interval,0],XiEta_ML[:pointsNum:interval,1],

symbol, markerfacecolor='none', #markeredgecolor='r',

markeredgewidth=2, markersize=10)

lg = plt.legend([p1,p2,p3], ['RANS', 'DNS', name], loc = 0)

lg.draw_frame(False)

plt.ylim([0,3**0.5/2.0])

plt.show()

代码

文本

[13]

def comparePlotRFNN(XiEta_RANS, testResponses, testResponsesPred_RF, testResponsesPred_NN):

XiEta_DNS = XiEta_RANS + testResponses

XiEta_RF = XiEta_RANS + testResponsesPred_RF

XiEta_NN = XiEta_RANS + testResponsesPred_NN

# Plot Reynolds stress anisotropy in Barycentric triangle

interval = 2

pointsNum = int(XiEta_RANS.shape[0])

plt.figure()

plt.plot([0,1,0.5,0.5,0],[0,0,3**0.5/2.0,3**0.5/2.0,0],'g-')

p1, = plt.plot(XiEta_RANS[:pointsNum:interval,0],XiEta_RANS[:pointsNum:interval,1],

'bo', markerfacecolor='none', markeredgecolor='b',

markeredgewidth=1.5, markersize=8)

p2, = plt.plot(XiEta_DNS[:pointsNum:interval,0],XiEta_DNS[:pointsNum:interval,1],

'ks', markerfacecolor='none', markeredgecolor='k',

markeredgewidth=1.5, markersize=8)

p3, = plt.plot(XiEta_RF[:pointsNum:interval,0],XiEta_RF[:pointsNum:interval,1],

'r^', markerfacecolor='none', markeredgecolor='r',

markeredgewidth=1.5, markersize=8)

p4, = plt.plot(XiEta_NN[:pointsNum:interval,0],XiEta_NN[:pointsNum:interval,1],

'r+', markerfacecolor='none', markeredgecolor='g',

markeredgewidth=1.5, markersize=8)

lg = plt.legend([p1,p2,p3, p4], ['RANS', 'DNS', 'RF', 'NN'], loc = 0)

lg.draw_frame(False)

plt.ylim([0,3**0.5/2.0])

plt.show()

代码

文本

[14]

def iterateLines(dataFolderRANS, testResponses, testResponsesPred, name, symbol='r^'):

# Start index of different sample lines

indexList = [0, 98, 191, 287, 385, 483, 581, 679, 777, 875, 971]

# Make plots at x=2 and x=4

for iterN in [3,5]:

XiEta = np.loadtxt(dataFolderRANS + 'line' + str(iterN) + '_XiEta.xy')

startIndex = indexList[iterN-1]

endIndex = indexList[iterN]

plotXiEta(XiEta, testResponses[startIndex:endIndex,:],

testResponsesPred[startIndex:endIndex,:], name, symbol)

#plt.show()

代码

文本

[15]

def compareResults(dataFolderRANS, testResponses, testResponsesPred_RF, testResponsesPred_NN):

## compare the results in one plot

# Start index of different sample lines

indexList = [0, 98, 191, 287, 385, 483, 581, 679, 777, 875, 971]

# Make plots at x=2 and x=4

for iterN in [3,5]:

XiEta = np.loadtxt(dataFolderRANS + 'line' + str(iterN) + '_XiEta.xy')

startIndex = indexList[iterN-1]

endIndex = indexList[iterN]

comparePlotRFNN(XiEta, testResponses[startIndex:endIndex,:],

testResponsesPred_RF[startIndex:endIndex,:],

testResponsesPred_NN[startIndex:endIndex,:])

代码

文本

Now, plot the anisotropy at the two locations $x / H = 2$ and 4:

代码

文本

[16]

Image(filename='/bohr/nbnov23-tuvx/v1/figs/locations.png')

代码

文本

[17]

# if __name__== "__main__":

# Load data

trainFeatures, trainResponses = loadTrainingData('pehill', 'Re5600')

testFeatures, testResponses = loadTestData('pehill', 'Re10595')

time_begin_RF = time.time()

# Make prediction via the random forest regressor

testResponsesPred_RF = randomForest(trainFeatures, trainResponses, testFeatures, 6, 100)

time_end_RF = time.time()

# Make plots of Reynolds stress anisotropy

dataFolderRANS = '/bohr/nbnov23-tuvx/v1/database/pehill/XiEta-RANS/Re10595/'

iterateLines(dataFolderRANS, testResponses, testResponsesPred_RF, name='RF')

plt.show()

代码

文本

[18]

Nepochs = 1000

time_begin_NN = time.time()

# Make prediction via the neural network

testResponsesPred_NN = keras_nn(trainFeatures, trainResponses, testFeatures, Nepochs)

time_end_NN = time.time()

31/31 [==============================] - 0s 1ms/step

代码

文本

Make plots of Reynolds stress anisotropy (NN results)

代码

文本

[19]

dataFolderRANS = '/bohr/nbnov23-tuvx/v1/database/pehill/XiEta-RANS/Re10595/'

symbol = 'g+'

iterateLines(dataFolderRANS, testResponses, testResponsesPred_NN, name='NN', symbol='m+')

plt.show()

代码

文本

Compare the results of random forest and neural network

代码

文本

[20]

compareResults(dataFolderRANS, testResponses, testResponsesPred_RF, testResponsesPred_NN)

plt.show()

代码

文本

Comparison of computational cost between RF and NN

The cost depends on the number of epoches, which is written in the title of the plot.

代码

文本

[21]

cost_time_RF = time_end_RF - time_begin_RF

cost_time_NN = time_end_NN - time_begin_NN

xlabel = np.arange(2)

plt.bar(xlabel, [cost_time_RF, cost_time_NN], 0.4)

plt.ylabel('CPU time (sec')

plt.xticks(xlabel, ('RF', 'NN'))

plt.title('Epoches = ' + str(Nepochs))

plt.show()

代码

文本

[ ]

代码

文本

Turbulence

AI4S

TurbulenceAI4S

已赞1