Inverse problem for the Navier-Stokes equation of incompressible flow around cylinder#
The problem we are solving#
We will solve an inverse problem for the Navier-Stokes equation of incompressible flow around a cylinder. Let us consider the Navier–Stokes equations in two dimensions2 (2D) given explicitly by
\[
\rho u_t + \rho λ_1(uu_x + vu_y) = −p_x + \rho λ_2(u_{xx} + u_{yy}),
\]
\[
\rho v_t + \rho λ_1(uv_x + vv_y) = −p_y + \rho λ_2(v_{xx} + v_{yy}),
\]
for \(x \in [1, 8]\), \(y \in [-2, 2]\), and \(t \in [0, 7]\) Spatio-temporal domain
where \(u(t, x, y)\) denotes the \(x\)-component of the velocity field, \(v(t, x, y)\) the y-component, \(\rho\) is the fluid density, and \(p(t, x, y)\) the pressure. Here, \(λ = (λ_1, λ_2)\) are the unknown parameters. Solutions to the Navier–Stokes equations are searched in the set of divergence-free functions; i.e.,
This extra equation is the continuity equation for incompressible fluids that describes the conservation of mass of the fluid. We make the assumption that
\[
u = ψ_y, v = −ψ_x,
\]
for some latent function \(ψ(t, x, y)\) Under this assumption, the continuity equation will be automatically satisfied. Given noisy measurements
\[
\{t^i, x^i , y^i, u^i , v^i \}^N_{i=1}
\]
of the velocity field, we are interested in learning the parameters \(λ\) as well as the pressure \(p(t, x, y)\). We define \(f (t, x, y)\) and \(g(t, x, y)\) to be given by
\[
f := u_t + λ_1(uu_x + vu_y) + p_x − λ_2(u_{xx} + u_{yy}),
\]
\[
g := v_t + λ_1(uv_x + vv_y) + p_y − λ_2(v_{xx} + v_{yy}),
\]
Dimensional Analysis#
Step 1: Assign Dimensions to Variables#
Spatial Coordinates \(x\) and \(y\):
The dimensions of \(x\) and \(y\) are length:
\[ [x] = [y] = L. \]
Time \(t\):
The dimension of time is:
\[ [t] = T. \]
Velocity Components \(u\) and \(v\):
Velocity components have dimensions of length per unit time:
\[ [u] = [v] = L / T. \]
Pressure \(p\):
Pressure has dimensions of force per unit area, which can be expressed as mass per unit length per unit time squared:
\[ [p] = M / (L T^2), \]where \(M\) represents mass.
Density \(\rho\):
Density has dimensions of mass per unit volume:
\[ [\rho] = M / L^3. \]
Parameters \(\lambda_1\) and \(\lambda_2\):
The term \(\rho \lambda_1(uu_x + vu_y)\) involves the advection term, which must have the same dimensions as the time derivative \(\rho u_t\).
The term \(\rho \lambda_2(u_{xx} + u_{yy})\) involves the diffusion term, which must have the same dimensions as the time derivative \(\rho u_t\).
Step 2: Analyze the Dimensions of Each Term#
Time Derivative Terms:
The time derivatives \(\rho u_t\) and \(\rho v_t\) have dimensions:
\[ [\rho u_t] = [\rho v_t] = [\rho] \cdot \frac{[u]}{[t]} = \frac{M}{L^3} \cdot \frac{L}{T^2} = \frac{M}{L^2 T^2}. \]
Advection Terms:
The advection terms \(\rho \lambda_1(uu_x + vu_y)\) and \(\rho \lambda_1(uv_x + vv_y)\) involve spatial derivatives of velocity:
\[ [uu_x] = [vu_y] = [uv_x] = [vv_y] = \frac{[u]^2}{[x]} = \frac{(L / T)^2}{L} = \frac{L}{T^2}. \]Therefore, the advection terms have dimensions:
\[ [\rho \lambda_1(uu_x + vu_y)] = [\rho \lambda_1(uv_x + vv_y)] = [\rho] \cdot \frac{L}{T^2} = \frac{M}{L^3} \cdot \frac{L}{T^2} = \frac{M}{L^2 T^2}. \]
Pressure Gradient Terms:
The pressure gradient terms \(-p_x\) and \(-p_y\) have dimensions:
\[ [p_x] = [p_y] = \frac{[p]}{[x]} = \frac{M / (L T^2)}{L} = \frac{M}{L^2 T^2}. \]
Diffusion Terms:
The diffusion terms \(\rho \lambda_2(u_{xx} + u_{yy})\) and \(\rho \lambda_2(v_{xx} + v_{yy})\) involve second spatial derivatives of velocity:
\[ [u_{xx}] = [u_{yy}] = [v_{xx}] = [v_{yy}] = \frac{[u]}{[x]^2} = \frac{L / T}{L^2} = \frac{1}{L T}. \]Therefore, the diffusion terms have dimensions:
\[ [\rho \lambda_2(u_{xx} + u_{yy})] = [\rho \lambda_2(v_{xx} + v_{yy})] = [\rho] \cdot [\lambda_2] \cdot \frac{1}{L T} = \frac{M}{L^3} \cdot \frac{L^2}{T} \cdot \frac{1}{L T} = \frac{M}{L^2 T^2}. \]
Step 3: Determine the Dimensions of \(\lambda_1\) and \(\lambda_2\)#
Advection Parameter \(\lambda_1\):
The advection term \(\rho \lambda_1(uu_x + vu_y)\) must have the same dimensions as \(\rho u_t\):
\[ [\rho \lambda_1] \cdot \frac{L}{T^2} = \frac{M}{L^2 T^2} \implies [\lambda_1] = \frac{M}{L^2 T^2} \cdot \frac{L^2 T^2}{M} = 1. \]Therefore, \(\lambda_1\) has no dimensions.
Diffusion Parameter \(\lambda_2\):
The diffusion term \(\rho \lambda_2(u_{xx} + u_{yy})\) must have the same dimensions as \(\rho u_t\):
\[ [\rho \lambda_2] \cdot \frac{1}{L T} = \frac{M}{L^2 T^2} \implies [\lambda_2] = \frac{M}{L^2 T^2} \cdot \frac{L^4 T}{M} = \frac{L^2}{T}. \]Therefore, \(\lambda_2\) has dimensions of mass per unit length (\(L^2 / T\)).
Step 4: Summary of Dimensions#
Variable/Parameter |
Physical Meaning |
Dimensions |
|---|---|---|
\(x, y\) |
Spatial coordinates |
\(L\) |
\(t\) |
Time |
\(T\) |
\(u, v\) |
Velocity components |
\(L / T\) |
\(p\) |
Pressure |
\(M / (L T^2)\) |
\(\rho\) |
Density |
\(M / L^3\) |
\(\lambda_1\) |
Advection parameter |
\(1\) |
\(\lambda_2\) |
Diffusion parameter |
\(L^2 / T\) |
Step 5: Given Values#
\(\lambda_1\) has no dimensions (\(M / L^3\)).
\(\lambda_2\) has dimensions of length \(^2\) per unit time (\(L^2 / T\)).
These dimensions are consistent with the Navier-Stokes equations for incompressible flow, incorporating the density term \(\rho\).
Implementation#
Below is the implementation of the inverse problem for the Navier-Stokes equation of incompressible flow around a cylinder using PINNx.
We first import the required libraries.
import re
import brainstate
import brainunit as u
import numpy as np
from scipy.io import loadmat
import jax
import matplotlib.pyplot as plt
import pinnx
Define the physical units for the problem.
unit_of_x = u.meter
unit_of_y = u.meter
unit_of_t = u.second
unit_of_rho = u.kilogram / u.meter**3
unit_of_c1 = u.UNITLESS
unit_of_c2 = u.meter2 / u.second
unit_of_u = u.meter / u.second
unit_of_v = u.meter / u.second
unit_of_p = u.pascal
Define Navier Stokes Equations:
def Navier_Stokes_Equation(x, y):
jacobian = net.jacobian(x)
x_hessian = net.hessian(x, y=['u', 'v'], xi=['x'], xj=['x'])
y_hessian = net.hessian(x, y=['u', 'v'], xi=['y'], xj=['y'])
u = y['u']
v = y['v']
p = y['p']
du_x = jacobian['u']['x']
du_y = jacobian['u']['y']
du_t = jacobian['u']['t']
dv_x = jacobian['v']['x']
dv_y = jacobian['v']['y']
dv_t = jacobian['v']['t']
dp_x = jacobian['p']['x']
dp_y = jacobian['p']['y']
du_xx = x_hessian['u']['x']['x']
du_yy = y_hessian['u']['y']['y']
dv_xx = x_hessian['v']['x']['x']
dv_yy = y_hessian['v']['y']['y']
continuity = du_x + dv_y
x_momentum = unit_of_rho * du_t + unit_of_rho * C1.value * (u * du_x + v * du_y) + dp_x - unit_of_rho * C2.value * (du_xx + du_yy)
y_momentum = unit_of_rho * dv_t + unit_of_rho * C1.value * (u * dv_x + v * dv_y) + dp_y - unit_of_rho * C2.value * (dv_xx + dv_yy)
return [continuity, x_momentum, y_momentum]
Define the neural network model:
net = pinnx.nn.Model(
pinnx.nn.DictToArray(x=unit_of_x, y=unit_of_y, t=unit_of_t),
pinnx.nn.FNN([3] + [50] * 6 + [3], "tanh"),
pinnx.nn.ArrayToDict(u=unit_of_u, v=unit_of_v, p=unit_of_p),
)
Define Spatio-temporal domain
# Rectangular
Lx_min, Lx_max = 1.0, 8.0
Ly_min, Ly_max = -2.0, 2.0
# Spatial domain: X × Y = [1, 8] × [−2, 2]
space_domain = pinnx.geometry.Rectangle([Lx_min, Ly_min], [Lx_max, Ly_max])
# Time domain: T = [0, 7]
time_domain = pinnx.geometry.TimeDomain(0, 7)
# Spatio-temporal domain
geomtime = pinnx.geometry.GeometryXTime(space_domain, time_domain)
geomtime = geomtime.to_dict_point(x=unit_of_x, y=unit_of_y, t=unit_of_t)
Define the true values of the \(\lambda_1, \lambda_2\) and the function to load training data.
# true values
C1true = 1.0
C2true = 0.01 * unit_of_c2
# Load training data
def load_training_data(num):
data = loadmat("../dataset/cylinder_nektar_wake.mat")
U_star = data["U_star"] # N x 2 x T
P_star = data["p_star"] # N x T
t_star = data["t"] # T x 1
X_star = data["X_star"] # N x 2
N = X_star.shape[0]
T = t_star.shape[0]
# Rearrange Problem
XX = np.tile(X_star[:, 0:1], (1, T)) # N x T
YY = np.tile(X_star[:, 1:2], (1, T)) # N x T
TT = np.tile(t_star, (1, N)).T # N x T
UU = U_star[:, 0, :] # N x T
VV = U_star[:, 1, :] # N x T
PP = P_star # N x T
x = XX.flatten()[:, None] # NT x 1
y = YY.flatten()[:, None] # NT x 1
t = TT.flatten()[:, None] # NT x 1
u = UU.flatten()[:, None] # NT x 1
v = VV.flatten()[:, None] # NT x 1
p = PP.flatten()[:, None] # NT x 1
# training domain: X × Y = [1, 8] × [−2, 2] and T = [0, 7]
data1 = np.concatenate([x, y, t, u, v, p], 1)
data2 = data1[:, :][data1[:, 2] <= 7]
data3 = data2[:, :][data2[:, 0] >= 1]
data4 = data3[:, :][data3[:, 0] <= 8]
data5 = data4[:, :][data4[:, 1] >= -2]
data_domain = data5[:, :][data5[:, 1] <= 2]
# choose number of training points: num =7000
idx = np.random.choice(data_domain.shape[0], num, replace=False)
x_train = data_domain[idx, 0]
y_train = data_domain[idx, 1]
t_train = data_domain[idx, 2]
u_train = data_domain[idx, 3]
v_train = data_domain[idx, 4]
p_train = data_domain[idx, 5]
# return [x_train, y_train, t_train, u_train, v_train, p_train]
return [x_train * unit_of_x, y_train * unit_of_y, t_train * unit_of_t,
u_train * unit_of_u, v_train * unit_of_v, p_train * unit_of_p]
Define the Parameters to be identified:
C1 = brainstate.ParamState(0.0)
C2 = brainstate.ParamState(0.0 * unit_of_c2)
Get the training data:
[ob_x, ob_y, ob_t, ob_u, ob_v, ob_p] = load_training_data(num=7000)
ob_xyt = {"x": ob_x, "y": ob_y, "t": ob_t}
ob_yv = {"u": ob_u, "v": ob_v, }
observe_bc = pinnx.icbc.PointSetBC(ob_xyt, ob_yv)
Define the problem and solve the inverse problem:
problem = pinnx.problem.TimePDE(
geomtime,
Navier_Stokes_Equation,
[observe_bc],
net,
num_domain=700,
num_boundary=200,
num_initial=100,
anchors=ob_xyt,
)
Setup the model trainer:
model = pinnx.Trainer(problem, external_trainable_variables=[C1, C2])
Define the callbacks for storing results:
fnamevar = "variables.dat"
variable = pinnx.callbacks.VariableValue([C1, C2], period=100, filename=fnamevar)
Compile, train and save trainer:
model.compile(braintools.optim.Adam(1e-3)).train(iterations=10000, callbacks=[variable],
display_every=1000, disregard_previous_best=True)
model.saveplot(issave=True, isplot=True)
model.compile(braintools.optim.Adam(1e-4)).train(iterations=10000, callbacks=[variable],
display_every=1000, disregard_previous_best=True)
model.saveplot(issave=True, isplot=True)
Compiling trainer...
'compile' took 0.096312 s
Training trainer...
Step Train loss Test loss Test metric
0 [0.28003943 * becquerel2, [0.28003943 * becquerel2, []
0.08384456 * (kilogram / klitre * (meter / second) / second) ** 2, 0.08384456 * (kilogram / klitre * (meter / second) / second) ** 2,
0.14824082 * (kilogram / klitre * (meter / second) / second) ** 2, 0.14824082 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.7059334 * meter / second, {'ibc0': {'u': 0.7059334 * meter / second,
'v': 1.947519 * meter / second}}] 'v': 1.947519 * meter / second}}]
1000 [0.00180862 * becquerel2, [0.00180862 * becquerel2, []
0.00279982 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00279982 * (kilogram / klitre * (meter / second) / second) ** 2,
0.00574982 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00574982 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.01001318 * meter / second, {'ibc0': {'u': 0.01001318 * meter / second,
'v': 0.02458837 * meter / second}}] 'v': 0.02458837 * meter / second}}]
2000 [0.00110196 * becquerel2, [0.00110196 * becquerel2, []
0.001189 * (kilogram / klitre * (meter / second) / second) ** 2, 0.001189 * (kilogram / klitre * (meter / second) / second) ** 2,
0.00160737 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00160737 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00394361 * meter / second, {'ibc0': {'u': 0.00394361 * meter / second,
'v': 0.00928425 * meter / second}}] 'v': 0.00928425 * meter / second}}]
3000 [0.00075236 * becquerel2, [0.00075236 * becquerel2, []
0.00077494 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00077494 * (kilogram / klitre * (meter / second) / second) ** 2,
0.00097048 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00097048 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00139976 * meter / second, {'ibc0': {'u': 0.00139976 * meter / second,
'v': 0.00430239 * meter / second}}] 'v': 0.00430239 * meter / second}}]
4000 [0.00042566 * becquerel2, [0.00042566 * becquerel2, []
0.0004626 * (kilogram / klitre * (meter / second) / second) ** 2, 0.0004626 * (kilogram / klitre * (meter / second) / second) ** 2,
0.00049891 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00049891 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00073693 * meter / second, {'ibc0': {'u': 0.00073693 * meter / second,
'v': 0.00118131 * meter / second}}] 'v': 0.00118131 * meter / second}}]
5000 [0.00027099 * becquerel2, [0.00027099 * becquerel2, []
0.00032169 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00032169 * (kilogram / klitre * (meter / second) / second) ** 2,
0.00030988 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00030988 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00051542 * meter / second, {'ibc0': {'u': 0.00051542 * meter / second,
'v': 0.0007421 * meter / second}}] 'v': 0.0007421 * meter / second}}]
6000 [0.00019304 * becquerel2, [0.00019304 * becquerel2, []
0.00024095 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00024095 * (kilogram / klitre * (meter / second) / second) ** 2,
0.00022211 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00022211 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00037054 * meter / second, {'ibc0': {'u': 0.00037054 * meter / second,
'v': 0.00047546 * meter / second}}] 'v': 0.00047546 * meter / second}}]
7000 [0.00014662 * becquerel2, [0.00014662 * becquerel2, []
0.00019332 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00019332 * (kilogram / klitre * (meter / second) / second) ** 2,
0.00016124 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00016124 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00027189 * meter / second, {'ibc0': {'u': 0.00027189 * meter / second,
'v': 0.00031798 * meter / second}}] 'v': 0.00031798 * meter / second}}]
8000 [0.00011527 * becquerel2, [0.00011527 * becquerel2, []
0.00016216 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00016216 * (kilogram / klitre * (meter / second) / second) ** 2,
0.00014464 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00014464 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00026401 * meter / second, {'ibc0': {'u': 0.00026401 * meter / second,
'v': 0.0002146 * meter / second}}] 'v': 0.0002146 * meter / second}}]
9000 [9.65281e-05 * becquerel2, [9.65281e-05 * becquerel2, []
0.00012925 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00012925 * (kilogram / klitre * (meter / second) / second) ** 2,
0.00011425 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00011425 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.0002264 * meter / second, {'ibc0': {'u': 0.0002264 * meter / second,
'v': 0.00015709 * meter / second}}] 'v': 0.00015709 * meter / second}}]
10000 [8.203531e-05 * becquerel2, [8.203531e-05 * becquerel2, []
0.00010743 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00010743 * (kilogram / klitre * (meter / second) / second) ** 2,
8.810719e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 8.810719e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00013619 * meter / second, {'ibc0': {'u': 0.00013619 * meter / second,
'v': 0.00012556 * meter / second}}] 'v': 0.00012556 * meter / second}}]
Best trainer at step 10000:
train loss: 5.39e-04
test loss: 5.39e-04
test metric: []
'train' took 91.177866 s
Saving loss history to /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/loss.dat ...
Saving checkpoint into /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/loss.dat
Saving training data to /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/train.dat ...
Saving checkpoint into /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/train.dat
Saving test data to /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/test.dat ...
Saving checkpoint into /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/test.dat
Compiling trainer...
'compile' took 0.068254 s
Training trainer...
Step Train loss Test loss Test metric
10000 [8.203531e-05 * becquerel2, [8.203531e-05 * becquerel2, []
0.00010743 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00010743 * (kilogram / klitre * (meter / second) / second) ** 2,
8.810719e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 8.810719e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00013619 * meter / second, {'ibc0': {'u': 0.00013619 * meter / second,
'v': 0.00012556 * meter / second}}] 'v': 0.00012556 * meter / second}}]
11000 [7.96016e-05 * becquerel2, [7.96016e-05 * becquerel2, []
0.00010362 * (kilogram / klitre * (meter / second) / second) ** 2, 0.00010362 * (kilogram / klitre * (meter / second) / second) ** 2,
8.5908185e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 8.5908185e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00012997 * meter / second, {'ibc0': {'u': 0.00012997 * meter / second,
'v': 0.00011991 * meter / second}}] 'v': 0.00011991 * meter / second}}]
12000 [7.62613e-05 * becquerel2, [7.62613e-05 * becquerel2, []
9.886667e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 9.886667e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
8.261804e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 8.261804e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00012325 * meter / second, {'ibc0': {'u': 0.00012325 * meter / second,
'v': 0.00011304 * meter / second}}] 'v': 0.00011304 * meter / second}}]
13000 [7.160572e-05 * becquerel2, [7.160572e-05 * becquerel2, []
9.2203576e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 9.2203576e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
7.7911514e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 7.7911514e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00011437 * meter / second, {'ibc0': {'u': 0.00011437 * meter / second,
'v': 0.00010453 * meter / second}}] 'v': 0.00010453 * meter / second}}]
14000 [6.619561e-05 * becquerel2, [6.619561e-05 * becquerel2, []
8.488368e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 8.488368e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
7.294095e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 7.294095e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 0.00010421 * meter / second, {'ibc0': {'u': 0.00010421 * meter / second,
'v': 9.4588875e-05 * meter / second}}] 'v': 9.4588875e-05 * meter / second}}]
15000 [6.082786e-05 * becquerel2, [6.082786e-05 * becquerel2, []
7.7803525e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 7.7803525e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
6.7992296e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 6.7992296e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 9.468454e-05 * meter / second, {'ibc0': {'u': 9.468454e-05 * meter / second,
'v': 8.5897096e-05 * meter / second}}] 'v': 8.5897096e-05 * meter / second}}]
16000 [5.6113255e-05 * becquerel2, [5.6113255e-05 * becquerel2, []
7.168196e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 7.168196e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
6.354173e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 6.354173e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 8.729596e-05 * meter / second, {'ibc0': {'u': 8.729596e-05 * meter / second,
'v': 7.8841156e-05 * meter / second}}] 'v': 7.8841156e-05 * meter / second}}]
17000 [5.2216485e-05 * becquerel2, [5.2216485e-05 * becquerel2, []
6.676306e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 6.676306e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
5.993753e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 5.993753e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 7.945047e-05 * meter / second, {'ibc0': {'u': 7.945047e-05 * meter / second,
'v': 7.2727926e-05 * meter / second}}] 'v': 7.2727926e-05 * meter / second}}]
18000 [4.8839163e-05 * becquerel2, [4.8839163e-05 * becquerel2, []
6.259088e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 6.259088e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
5.6764642e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 5.6764642e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 7.3541385e-05 * meter / second, {'ibc0': {'u': 7.3541385e-05 * meter / second,
'v': 6.803875e-05 * meter / second}}] 'v': 6.803875e-05 * meter / second}}]
19000 [4.608007e-05 * becquerel2, [4.608007e-05 * becquerel2, []
5.8930185e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 5.8930185e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
5.403013e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 5.403013e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 6.8380396e-05 * meter / second, {'ibc0': {'u': 6.8380396e-05 * meter / second,
'v': 6.40664e-05 * meter / second}}] 'v': 6.40664e-05 * meter / second}}]
20000 [4.3594566e-05 * becquerel2, [4.3594566e-05 * becquerel2, []
5.553491e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 5.553491e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
5.129298e-05 * (kilogram / klitre * (meter / second) / second) ** 2, 5.129298e-05 * (kilogram / klitre * (meter / second) / second) ** 2,
{'ibc0': {'u': 6.40729e-05 * meter / second, {'ibc0': {'u': 6.40729e-05 * meter / second,
'v': 6.0597904e-05 * meter / second}}] 'v': 6.0597904e-05 * meter / second}}]
Best trainer at step 20000:
train loss: 2.75e-04
test loss: 2.75e-04
test metric: []
'train' took 76.299114 s
Saving loss history to /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/loss.dat ...
Saving checkpoint into /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/loss.dat
Saving training data to /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/train.dat ...
Saving checkpoint into /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/train.dat
Saving test data to /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/test.dat ...
Saving checkpoint into /home/brainpy/codes/SichaoHe/pinnx/docs/unit-examples-inverse/test.dat
Calculate the mean residual error:
f = model.predict(ob_xyt, operator=Navier_Stokes_Equation)
print("Mean residual:", jax.tree.map(lambda x: u.math.mean(u.math.abs(x)), f))
Mean residual: [0.00478853 * becquerel, 0.00537193 * kilogram / klitre * (meter / second) / second, 0.0051398 * kilogram / klitre * (meter / second) / second]
Plot Variables:
# reopen saved data using callbacks in fnamevar
lines = open(fnamevar, "r").readlines()
# read output data in fnamevar
Chat = np.array(
[
np.fromstring(
min(re.findall(re.escape("[") + "(.*?)" + re.escape("]"), line), key=len),
sep=",",
)
for line in lines
]
)
l, c = Chat.shape
plt.semilogy(range(0, l * 100, 100), Chat[:, 0], "r-")
plt.semilogy(range(0, l * 100, 100), Chat[:, 1], "k-")
plt.semilogy(range(0, l * 100, 100), np.ones(Chat[:, 0].shape) * C1true, "r--")
plt.semilogy(range(0, l * 100, 100), np.ones(Chat[:, 1].shape) * C2true, "k--")
plt.legend(["C1hat", "C2hat", "True C1", "True C2"], loc="right")
plt.xlabel("Epochs")
plt.title("Variables")
plt.show()
/tmp/ipykernel_2729976/4284654551.py:6: DeprecationWarning: string or file could not be read to its end due to unmatched data; this will raise a ValueError in the future.
np.fromstring(
Plot the velocity distribution of the flow field:
for t in range(0, 8):
t = t * unit_of_t
[ob_x, ob_y, ob_t, ob_u, ob_v, ob_p] = load_training_data(num=140000)
xyt_pred = {"x": ob_x, "y": ob_y, "t": t * np.ones((len(ob_x),))}
uvp_pred = model.predict(xyt_pred)
x_pred, y_pred, t_pred = xyt_pred['x'], xyt_pred['y'], xyt_pred['t']
u_pred, v_pred, p_pred = uvp_pred['u'], uvp_pred['v'], uvp_pred['p']
x_true = ob_x[ob_t == t]
y_true = ob_y[ob_t == t]
u_true = ob_u[ob_t == t]
fig, ax = plt.subplots(2, 1)
cntr0 = ax[0].tricontourf(x_pred.mantissa, y_pred.mantissa, u_pred.mantissa, levels=80, cmap="rainbow")
cb0 = plt.colorbar(cntr0, ax=ax[0])
cntr1 = ax[1].tricontourf(x_true.mantissa, y_true.mantissa, u_true.mantissa, levels=80, cmap="rainbow")
cb1 = plt.colorbar(cntr1, ax=ax[1])
ax[0].set_title("u-PINN " + "(t=" + str(t) + ")", fontsize=9.5)
ax[0].axis("scaled")
ax[0].set_xlabel("X", fontsize=7.5, family="Arial")
ax[0].set_ylabel("Y", fontsize=7.5, family="Arial")
ax[1].set_title("u-Reference solution " + "(t=" + str(t) + ")", fontsize=9.5)
ax[1].axis("scaled")
ax[1].set_xlabel("X", fontsize=7.5, family="Arial")
ax[1].set_ylabel("Y", fontsize=7.5, family="Arial")
fig.tight_layout()
plt.show()
