deepxde
deepxde copied to clipboard
lululxvi
Flow around an airfoil - incompressible NS. Low residuals but non-physical result obtained + some possible not related geometry bugs
I'm facing a problem with the result obtained from the DeepXDE library compared with the literature. It may be probably a mistake of mine, but I think it may be also helpful a discussion for everyone interested in the matter.
As anticipated in the title, I'm solving the non-dimensional steady incompressible NS system of PDEs written in the V-P formulation in order to obtain the behaviour of the flow around an airfoil.
In particular, I'm trying to replicate the results showed in this paper: https://arc.aiaa.org/doi/abs/10.2514/6.2022-0187, following some tips for the training that I found in: https://arxiv.org/pdf/2003.06496.pdf.
The setup of the problem is as follows:
- u component = 1.0 on the inlet, top, bottom of the farfield
- v component = 0.0 on the inlet, top, bottom of the farfield
- u and v components = 0.0 on the airfoil boundary (no-slip condition)
All conditions are of Dirichlet type. No condition was set on p, as done in the second paper I linked (note: the problem persists also if a Dirichlet boundary condition is prescribed for p at the outlet).
In order to try to simplify the learning process, I enforced the boundary conditions on the farfield in the hard way.
Here is the code I used, hopefully well commented:
import deepxde as dde
import numpy as np
import tensorflow as tf
from matplotlib import pyplot as plt
dde.config.set_random_seed(48)
dde.config.set_default_float('float64')
xmin, xmax = -4.0, 6.0
ymin, ymax = -4.0, 4.0
def boundaryNACA4D(M, P, SS, c, n):
"""
Compute the coordinates of a NACA 4-digits airfoil
Args:
M: maximum camber value (*100)
P: position of the maximum camber (*10)
SS: maximum thickness (*100)
c: chord length
n: the total points sampled will be 2*n
"""
m = M / 100
p = P / 10
t = SS / 100
if (m == 0):
p = 1
# Chord discretization (cosine discretization)
xv = np.linspace(0.0, c, n+1)
xv = c / 2.0 * (1.0 - np.cos(np.pi * xv / c))
# Thickness distribution
ytfcn = lambda x: 5 * t * c * (0.2969 * (x / c)**0.5 -
0.1260 * (x / c) -
0.3516 * (x / c)**2 +
0.2843 * (x / c)**3 -
0.1015 * (x / c)**4)
yt = ytfcn(xv)
# Camber line
yc = np.zeros(np.size(xv))
for ii in range(n+1):
if xv[ii] <= p * c:
yc[ii] = c * (m / p**2 * (xv[ii] / c) * (2 * p - (xv[ii] / c)))
else:
yc[ii] = c * (m / (1 - p)**2 * (1 + (2 * p - (xv[ii] / c)) * (xv[ii] / c) - 2 * p))
# Camber line slope
dyc = np.zeros(np.size(xv))
for ii in range(n+1):
if xv[ii] <= p * c:
dyc[ii] = m / p**2 * 2 * (p - xv[ii] / c)
else:
dyc[ii] = m / (1 - p)**2 * 2 * (p - xv[ii] / c)
# Boundary coordinates and sorting
th = np.arctan2(dyc, 1)
xU = xv - yt * np.sin(th)
yU = yc + yt * np.cos(th)
xL = xv + yt * np.sin(th)
yL = yc - yt * np.cos(th)
x = np.zeros(2 * n + 1)
y = np.zeros(2 * n + 1)
for ii in range(n):
x[ii] = xL[n - ii]
y[ii] = yL[n - ii]
x[n : 2 * n + 1] = xU
y[n : 2 * n + 1] = yU
return np.vstack((x, y)).T
def navier_stokes(x, y):
"""
System of PDEs to be minimized: incompressible Navier-Stokes equations (V-P formulation)
"""
Re = 50 # Reynolds number
u, v, p = y[:, 0:1], y[:, 1:2], y[:, 2:3]
u_x = dde.grad.jacobian(y, x, i = 0, j = 0)
u_y = dde.grad.jacobian(y, x, i = 0, j = 1)
v_x = dde.grad.jacobian(y, x, i = 1, j = 0)
v_y = dde.grad.jacobian(y, x, i = 1, j = 1)
u_xx = dde.grad.hessian(y, x, i = 0, j = 0, component = 0)
u_yy = dde.grad.hessian(y, x, i = 1, j = 1, component = 0)
v_xx = dde.grad.hessian(y, x, i = 0, j = 0, component = 1)
v_yy = dde.grad.hessian(y, x, i = 1, j = 1, component = 1)
p_x = dde.grad.jacobian(y, x, i = 2, j = 0)
p_y = dde.grad.jacobian(y, x, i = 2, j = 1)
continuity = u_x + v_y
momentum_x = u * u_x + v * u_y + p_x - 1.0 / Re * (u_xx + u_yy)
momentum_y = u * v_x + v * v_y + p_y - 1.0 / Re * (v_xx + v_yy)
return [continuity, momentum_x, momentum_y]
# Geometry definition
farfield = dde.geometry.Rectangle([xmin, ymin], [xmax, ymax])
airfoil = dde.geometry.Polygon(boundaryNACA4D(0, 0, 12, 1, 100)) # Generating a NACA0012 airfoil
geom = dde.geometry.CSGDifference(farfield, airfoil)
# Boundaries definition
def boundary_farfield_inlet(x, on_boundary):
return on_boundary and np.isclose(x[0], xmin)
def boundary_farfield_top(x, on_boundary):
return on_boundary and np.isclose(x[1], ymax)
def boundary_farfield_outlet(x, on_boundary):
return on_boundary and np.isclose(x[0], xmax)
def boundary_farfield_bottom(x, on_boundary):
return on_boundary and np.isclose(x[1], ymin)
def boundary_airfoil(x, on_boundary):
# Note: return on_boundary and airfoil.on_boundary(x) gives an error (Problem in dde Polygon method on_boundary?)
return on_boundary and (not farfield.on_boundary(x))
# Boundary values definition
def funU(x):
return 1.0
def funV(x):
return 0.0
#def funP(x):
# return 0.0
# Boundary conditions assembly
bc_inlet_u = dde.DirichletBC(geom, funU, boundary_farfield_inlet, component = 0)
bc_inlet_v = dde.DirichletBC(geom, funV, boundary_farfield_inlet, component = 1)
bc_top_u = dde.DirichletBC(geom, funU, boundary_farfield_top, component = 0)
bc_top_v = dde.DirichletBC(geom, funV, boundary_farfield_top, component = 1)
bc_bottom_u = dde.DirichletBC(geom, funU, boundary_farfield_bottom, component = 0)
bc_bottom_v = dde.DirichletBC(geom, funV, boundary_farfield_bottom, component = 1)
#bc_outlet_p = dde.DirichletBC(geom, funP, boundary_farfield_outlet, component = 2)
#bc_outlet_u = dde.NeumannBC(geom, funV, boundary_farfield_outlet, component = 0)
bc_airfoil_u = dde.DirichletBC(geom, funV, boundary_airfoil, component = 0)
bc_airfoil_v = dde.DirichletBC(geom, funV, boundary_airfoil, component = 1)
bcs = [bc_inlet_u, bc_inlet_v, bc_top_u, bc_top_v, bc_bottom_u, bc_bottom_v, bc_airfoil_u, bc_airfoil_v]
# Problem setup
data = dde.data.PDE(geom, navier_stokes, bcs, num_domain = 10000, num_boundary = 7500, num_test = 5000)
# Neural network definition
layer_size = [2] + [50] * 6 + [3]
n = 10
activation = f"LAAF-{n} tanh"
initializer = 'Glorot uniform'
net = dde.nn.FNN(layer_size, activation, initializer)
# Enforcing hard boundary conditions where easily possible
def modify_output(X, Y):
x, y = X[:, 0:1], X[:, 1:2]
u, v, p = Y[:, 0:1], Y[:, 1:2], Y[:, 2:3]
u_new = (x - xmin) * (y - ymin) * (y - ymax) * u / 100 + 1.0 #Note: divided by 100 to balance the effect of pre-multiplications
v_new = (x - xmin) * (y - ymin) * (y - ymax) * v / 20
return tf.concat((u_new, v_new, p), axis=1)
net.apply_output_transform(modify_output)
# Model definition
model = dde.Model(data, net)
model.compile(optimizer = 'adam', lr = 1e-3, loss_weights = [1, 1, 1, 10, 10, 10, 10, 10, 10, 10, 10]) # Giving more weight to bcs (actually no needed for hard imposed ones)
#resampler = dde.callbacks.PDEResidualResampler(period=1000)
#early_stopping = dde.callbacks.EarlyStopping(patience = 40000)
# Training strategy:
# 5000 epochs with adam and lr = 1e-3
# 5000 epochs with adam and lr = 1e-4
# 10000 epochs with adam and lr = 1e-5
# 10000 epochs with adam and lr = 1e-6
# L-BFGS-B at the end to fine tuning the network parameters
losshistory, train_state = model.train(epochs = 5000, display_every = 100, model_save_path = './')
dde.saveplot(losshistory, train_state, issave = True, isplot = True)
model.compile(optimizer = 'adam', lr = 1e-4, loss_weights = [1, 1, 1, 10, 10, 10, 10, 10, 10, 10, 10])
losshistory, train_state = model.train(epochs = 5000, display_every = 100, model_save_path = './')
dde.saveplot(losshistory, train_state, issave = True, isplot = True)
model.compile(optimizer = 'adam', lr = 1e-5, loss_weights = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
losshistory, train_state = model.train(epochs = 10000, display_every = 100, model_save_path = './')
dde.saveplot(losshistory, train_state, issave = True, isplot = True)
model.compile(optimizer = 'adam', lr = 1e-6, loss_weights = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
losshistory, train_state = model.train(epochs = 10000, display_every = 100, model_save_path = './')
dde.saveplot(losshistory, train_state, issave = True, isplot = True)
model.compile(optimizer = 'L-BFGS-B', loss_weights = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
losshistory, train_state = model.train(display_every = 100, model_save_path = './')
dde.saveplot(losshistory, train_state, issave = True, isplot = True)
# Plotting tool: thanks to @q769855234 code snippet in DeepXDE GitHub issues directory
dx = 0.01
dy = 0.01
dt = 0.01
x = np.arange(xmin, xmax + dy, dx)
y = np.arange(ymin, ymax + dy, dy)
X = np.zeros((len(x)*len(y), 2))
xs = np.vstack((x,)*len(y)).reshape(-1)
ys = np.vstack((y,)*len(x)).T.reshape(-1)
X[:, 0] = xs
X[:, 1] = ys
# Model predictions generation
Y = model.predict(X)
u = Y[:, 0].reshape(len(y), len(x))
v = Y[:, 1].reshape(len(y), len(x))
p = Y[:, 2].reshape(len(y), len(x))
for i in range(Y.shape[1]):
plt.figure(figsize = (16, 9))
plt.streamplot(x, y, u, v, density = 1.5)
plt.contourf(x, y, Y[:, i].reshape(len(y), len(x)))
plt.plot(boundaryNACA4D(0, 0, 12, 1, 100)[:, 0], boundaryNACA4D(0, 0, 12, 1, 100)[:, 1])
plt.colorbar()
plt.savefig('NACA0012NS{}.png'.format(i))
I can easily obtain an overall residual lower than 1e-4, so I can say I'm satisfied by the PINN ability to minimize the residuals (I read that as a rule of thumb residuals in the order of 1e-4 are enough to consider the convergence reached).
The problem is that what I can see in the plots of the results is different from what I expect from the physics of the problem (you can look at the second paper as a reference).
Pressure field:
u velocity component:
v velocity component:
As you can see in the attached picture, the point in which the pressure has a maximum is far ahead the leading edge of the airfoil, where, as the physics prescribes there should be null velocity and thus maximum pressure. The behaviour of the velocity is strange as well. It seems like that the presence of the airfoil is felt by the flow far before it arrives close to the airfoil surface.
The maximum of the pressure field, the decrease of u and the increase in the module of v, all of them happen too much ahead of the airfoil leading edge. What I can see is that the airfoil is felt by the flow, but too much, as it was bigger (I'm sure the geometry of the airfoil is correct anyway, since it is also correctly plotted).
What I have already tried:
- Adding more points in the domain and in the boundaries
- Increase the hidden layers of the network up to 10
- Different activation functions
- Different training schedules
- Impose all the boundary conditions in the soft way
Am I missing something? Wrong formulation of boundary conditions? Any problem in the code? Any help is much appreciated!
Thanks @lululxvi for this amazing library!
Not related to previous problem:
As concern DeepXDE, I may have found some bugs in the Polygon class. In particular, if define the airfoil boundary condition as:
def boundary_airfoil(x, on_boundary):
return on_boundary and airfoil.on_boundary(x)
I get the following error:
File "C:\Users\Riccardo\anaconda3\envs\spyder-env\lib\site-packages\deepxde\icbc\boundary_conditions.py", line 47, in filter
return X[self.on_boundary(X, self.geom.on_boundary(X))]
IndexError: arrays used as indices must be of integer (or boolean) type
Moreover, also if for some reasons (not related to the case in exam) I need to use the method inside(x), I obtain an error as well:
def boundary_airfoil(x, on_boundary):
return on_boundary and airfoil.inside(x)
File "C:\Users\Riccardo\anaconda3\envs\spyder-env\lib\site-packages\deepxde\geometry\geometry_2d.py", line 376, in wn_PnPoly
V[i, 1] <= P[:, 1:2], # Start y <= P[1]
IndexError: too many indices for array: array is 1-dimensional, but 2 were indexe
It seems something wrong with the arrays shape inside the Winding number algorithm
Seems good implementation. I didn't solve this problem by myself, so I am not sure what the problem is. Could you try a smaller domain?
Use
def boundary_airfoil(x, on_boundary):
return on_boundary and airfoil.on_boundary(np.array([x]))[0]
Similarly for inside
Thank you very much for your kind reply. I'm going to work on it and keep you and anyone interested updated
Why do you have multiple model compilation with different learning rate? Adam already changes the learning rate based on the momentum values. Is there an advantage of doing this?
You are right, Adam adapts the learning rate of each parameter already on its own. But the value you set when you compile the model is an upper bound. What I can say from my own general experience (I'm new to PINNs, but I followed a university course on "classical" neural networks) is that using some sort of weight decay schedule the risk of overfitting is reduced or postponed. In the specific case of PINNs I noticed a reduction in the residuals instability. Anyway in my opinion is something someone can try while fine tuning a model which already works well in order to increase the performances, it's not that important in my case where I'm facing more fundamental problems
Update:
First of all, thank you @lululxvi for having clarified me how to correctly use the functions on_boundary() and inside() when dealing with a Polygon. I think this may also better answer to issue #573.
As concern the problem itself, I did quite a lot changes to the code, such as:
- Reducing the computational domain as kindly suggested by @lululxvi .
- Changing the formulation of the governing equations, as suggested in https://arxiv.org/pdf/2002.10558.pdf, passing from the classical V-P NS formulation to the continuum-mechanichs based one.
- Using the stream function psi as output variable rather than u and v. In this way the divergence free condition (mass conservation) is exactly satisfied. On the other hand this alternative approach complicates the imposition of boundary conditions in the hard way (I still didn't have a look in depth on how to implement them now), so in the end I'm not sure about the overall net gain of this strategy.
- Increasing the total number of collocation points to 50k + bc ones.
- Refining the collocation points in the portion of the domain close to the airfoil leading edge, where the maximum of velocity gradients are expected, in the way you can see in the figure:
I still didn't fine tuned the model and I trained it for 10k epochs with Adam and 50k epochs with L-BFGS-B. The overall residual I obtain is greater than the one I got in the notebook above, but this time the solution it looks definetly as expected. Here I leave some pictures:
-
Loss:
-
Pressure field:
-
u velocity component:
-
v velocity component:
I will continue working on it to improve the results and comparing them with some CFD solver.
I think a few possible strategies for achieving better results could be:
- Increase the number of training epochs
- Tune the hyperparameters of the network
- Implement an adaptive loss weights algorithm, but I think I still don't have the knowledge to do it on my own.
- Code a gPINN. That paper is really interesting @lululxvi, but with this formulation I would get 25 losses to be minimized, maybe the problem becomes too complicated for the optimizer? I may try anyway.
- Force at least some boundary condition in the hard way
Here I leave also the code I have used to obtain the results if it may help anyone, for example #629:
import deepxde as dde
import numpy as np
import tensorflow as tf
from matplotlib import pyplot as plt
dde.config.set_random_seed(48)
dde.config.set_default_float('float64')
xmin, xmax = 0.0, 1.0
ymin, ymax = 0.0, 0.7
rho = 1.0
mu = 0.02
umax = 1.0
def boundaryNACA4D(M, P, SS, c, n, offset_x, offset_y):
"""
Compute the coordinates of a NACA 4-digits airfoil
Args:
M: maximum camber value (*100)
P: position of the maximum camber alog the chord (*10)
SS: maximum thickness (*100)
c: chord length
n: the total points sampled will be 2*n
"""
m = M / 100
p = P / 10
t = SS / 100
if (m == 0):
p = 1
# Chord discretization (cosine discretization)
xv = np.linspace(0.0, c, n+1)
xv = c / 2.0 * (1.0 - np.cos(np.pi * xv / c))
# Thickness distribution
ytfcn = lambda x: 5 * t * c * (0.2969 * (x / c)**0.5 -
0.1260 * (x / c) -
0.3516 * (x / c)**2 +
0.2843 * (x / c)**3 -
0.1015 * (x / c)**4)
yt = ytfcn(xv)
# Camber line
yc = np.zeros(np.size(xv))
for ii in range(n+1):
if xv[ii] <= p * c:
yc[ii] = c * (m / p**2 * (xv[ii] / c) * (2 * p - (xv[ii] / c)))
else:
yc[ii] = c * (m / (1 - p)**2 * (1 + (2 * p - (xv[ii] / c)) * (xv[ii] / c) - 2 * p))
# Camber line slope
dyc = np.zeros(np.size(xv))
for ii in range(n+1):
if xv[ii] <= p * c:
dyc[ii] = m / p**2 * 2 * (p - xv[ii] / c)
else:
dyc[ii] = m / (1 - p)**2 * 2 * (p - xv[ii] / c)
# Boundary coordinates and sorting
th = np.arctan2(dyc, 1)
xU = xv - yt * np.sin(th)
yU = yc + yt * np.cos(th)
xL = xv + yt * np.sin(th)
yL = yc - yt * np.cos(th)
x = np.zeros(2 * n + 1)
y = np.zeros(2 * n + 1)
for ii in range(n):
x[ii] = xL[n - ii]
y[ii] = yL[n - ii]
x[n : 2 * n + 1] = xU
y[n : 2 * n + 1] = yU
return np.vstack((x + offset_x, y + offset_y)).T
def navier_stokes(x, y):
"""
System of PDEs to be minimized: incompressible Navier-Stokes equation in the
continuum-mechanics based formulations.
"""
psi, p, sigma11, sigma22, sigma12 = y[:, 0:1], y[:, 1:2], y[:, 2:3], y[:, 3:4], y[:, 4:5]
u = dde.grad.jacobian(y, x, i = 0, j = 1)
v = - dde.grad.jacobian(y, x, i = 0, j = 0)
u_x = dde.grad.jacobian(u, x, i = 0, j = 0)
u_y = dde.grad.jacobian(u, x, i = 0, j = 1)
v_x = dde.grad.jacobian(v, x, i = 0, j = 0)
v_y = dde.grad.jacobian(v, x, i = 0, j = 1)
sigma11_x = dde.grad.jacobian(y, x, i = 2, j = 0)
sigma12_x = dde.grad.jacobian(y, x, i = 4, j = 0)
sigma12_y = dde.grad.jacobian(y, x, i = 4, j = 1)
sigma22_y = dde.grad.jacobian(y, x, i = 3, j = 1)
continuumx = rho * (u * u_x + v * u_y) - sigma11_x - sigma12_y
continuumy = rho * (u * v_x + v * v_y) - sigma12_x - sigma22_y
constitutive1 = - p + 2 * mu * u_x - sigma11
constitutive2 = - p + 2 * mu * v_y - sigma22
constitutive3 = mu * (u_y + v_x) - sigma12
constitutive4 = p + (sigma11 + sigma22) / 2
return continuumx, continuumy, constitutive1, constitutive2, constitutive3, constitutive4
# Geometry defintion
farfield = dde.geometry.Rectangle([xmin, ymin], [xmax, ymax])
airfoil = dde.geometry.Polygon(boundaryNACA4D(0, 0, 12, 0.2, 250, 0.20, 0.35))
geom = dde.geometry.CSGDifference(farfield, airfoil)
inner_rec = dde.geometry.Rectangle([0.15, 0.28], [0.28, 0.42])
outer_dom = dde.geometry.CSGDifference(farfield, inner_rec)
outer_dom = dde.geometry.CSGDifference(outer_dom, airfoil)
inner_dom = dde.geometry.CSGDifference(inner_rec, airfoil)
inner_points = inner_dom.random_points(10000)
outer_points = outer_dom.random_points(40000)
farfield_points = farfield.random_boundary_points(1280)
airfoil_points = boundaryNACA4D(0, 0, 12, 0.2, 125, 0.20, 0.35)
points = np.append(inner_points, outer_points, axis = 0)
points = np.append(points, farfield_points, axis = 0)
points = np.append(points, airfoil_points, axis = 0)
# Boundaries definition
def boundary_farfield_inlet(x, on_boundary):
return on_boundary and np.isclose(x[0], xmin)
def boundary_farfield_top_bottom(x, on_boundary):
return on_boundary and (np.isclose(x[1], ymax) or np.isclose(x[1], ymin))
def boundary_farfield_outlet(x, on_boundary):
return on_boundary and np.isclose(x[0], xmax)
def boundary_airfoil(x, on_boundary):
return on_boundary and (not farfield.on_boundary(x))
# Boundary values definition
def fun_u_farfield(x, y, _):
return dde.grad.jacobian(y, x, i = 0, j = 1) - 1.0
def fun_no_slip_u(x, y, _):
return dde.grad.jacobian(y, x, i = 0, j = 1)
def fun_no_slip_v(x, y, _):
return - dde.grad.jacobian(y, x, i = 0, j = 0)
def funP(x):
return 0.0
# Boundary conditions assembly
bc_inlet_u = dde.OperatorBC(geom, fun_u_farfield, boundary_farfield_inlet)
bc_inlet_v = dde.OperatorBC(geom, fun_no_slip_v, boundary_farfield_inlet)
bc_top_bottom_u = dde.OperatorBC(geom, fun_u_farfield, boundary_farfield_top_bottom)
bc_top_bottom_v = dde.OperatorBC(geom, fun_no_slip_v, boundary_farfield_top_bottom)
bc_outlet_p = dde.DirichletBC(geom, funP, boundary_farfield_outlet, component = 1)
bc_airfoil_u = dde.OperatorBC(geom, fun_no_slip_u, boundary_airfoil)
bc_airfoil_v = dde.OperatorBC(geom, fun_no_slip_v, boundary_airfoil)
bcs = [bc_inlet_u, bc_inlet_v, bc_top_bottom_u, bc_top_bottom_v, bc_outlet_p, bc_airfoil_u, bc_airfoil_v]
# Problem setup
data = dde.data.PDE(geom, navier_stokes, bcs, num_domain = 0, num_boundary = 0, num_test = 5000, anchors = points)
plt.figure(figsize = (32, 18))
plt.scatter(data.train_x_all[:,0], data.train_x_all[:,1], s = 0.3)
plt.axis('equal')
dde.config.set_default_float('float64')
# Neural network definition
layer_size = [2] + [40] * 8 + [5]
activation = 'tanh'
initializer = 'Glorot uniform'
net = dde.nn.FNN(layer_size, activation, initializer)
# Model definition
model = dde.Model(data, net)
model.compile(optimizer = 'adam', lr = 5e-4, loss_weights = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2]) # Giving more weight to bcs
losshistory, train_state = model.train(epochs = 10000, display_every = 100, model_save_path = './')
dde.saveplot(losshistory, train_state, issave = True, isplot = True)
model.compile(optimizer = 'L-BFGS-B', loss_weights = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2])
model.train_step.optimizer_kwargs = {'options': {'maxcor': 50,
'ftol': 1.0 * np.finfo(float).eps,
'maxfun': 50000,
'maxiter': 50000,
'maxls': 50}}
losshistory, train_state = model.train(display_every = 100, model_save_path = './')
dde.saveplot(losshistory, train_state, issave = True, isplot = True)
# Plotting tool: thanks to @q769855234 code snippet
dx = 0.01
dy = 0.01
x = np.arange(xmin, xmax + dy, dx)
y = np.arange(ymin, ymax + dy, dy)
X = np.zeros((len(x)*len(y), 2))
xs = np.vstack((x,)*len(y)).reshape(-1)
ys = np.vstack((y,)*len(x)).T.reshape(-1)
X[:, 0] = xs
X[:, 1] = ys
def getU(x, y):
return dde.grad.jacobian(y, x, i = 0, j = 1)
def getV(x, y):
return - dde.grad.jacobian(y, x, i = 0, j = 0)
def getP(x, y):
return y[:, 1:2]
# Model predictions generation
u = model.predict(X, operator = getU)
v = model.predict(X, operator = getV)
p = model.predict(X, operator = getP)
#for i in range(len(X)):
# if airfoil.inside(np.array([X[i]]))[0]:
# u[i] = 0.0
# v[i] = 0.0
u = u.reshape(len(y), len(x))
v = v.reshape(len(y), len(x))
p = p.reshape(len(y), len(x))
airfoil_plot = boundaryNACA4D(0, 0, 12, 0.2, 150, offset_x = 0.20, offset_y = 0.35)
fig1, ax1 = plt.subplots(figsize = (16, 9))
#ax1.streamplot(x, y, u, v, density = 1.5)
clev = np.arange(p.min(), 2, 0.001)
cnt1 = ax1.contourf(x, y, p, clev, cmap = plt.cm.coolwarm, extend='both')
plt.axis('equal')
plt.fill(airfoil_plot[:, 0], airfoil_plot[:, 1])
fig1.colorbar(cnt1)
plt.savefig('NACA0012NS0.png')
fig2, ax2 = plt.subplots(figsize = (16, 9))
#ax2.streamplot(x, y, u, v, density = 1.5)
clev = np.arange(0, u.max(), 0.001)
cnt2 = ax2.contourf(x, y, u, clev, cmap = plt.cm.coolwarm, extend='both')
plt.axis('equal')
plt.fill(airfoil_plot[:, 0], airfoil_plot[:, 1])
fig2.colorbar(cnt2)
plt.savefig('NACA0012NS1.png')
fig3, ax3 = plt.subplots(figsize = (16, 9))
#ax3.streamplot(x, y, u, v, density = 1.5)
clev = np.arange(-0.235, v.max(), 0.001)
cnt3 = ax3.contourf(x, y, v, clev, cmap = plt.cm.coolwarm, extend='both')
plt.axis('equal')
plt.fill(airfoil_plot[:, 0], airfoil_plot[:, 1])
fig3.colorbar(cnt3)
plt.savefig('NACA0012NS2.png')
@riccardotomada Sounds great. Once you finalize it, you are welcome to submit a pull request to add the code in the examples of DeepXDE.