lululxvi / deepxde

A library for scientific machine learning and physics-informed learning
https://deepxde.readthedocs.io
GNU Lesser General Public License v2.1
2.64k stars 739 forks source link

my case can't converge #629

Open q769855234 opened 2 years ago

q769855234 commented 2 years ago

Dear sir, I try to sove a case, flowing around the square column, but the solution is incorrect, I hope you can provide some help, my code is as follows: import numpy as np from deepxde.gradients import jacobian, hessian from deepxde.geometry import Rectangle, TimeDomain, GeometryXTime, CSGDifference from deepxde import IC, OperatorBC, DirichletBC, NeumannBC from deepxde.data import TimePDE from deepxde import maps, Model from matplotlib.pyplot import streamplot, contourf, colorbar,savefig from deepxde import saveplot from matplotlib.pyplot import plot, pause, clf, figure import torch import tensorflow as tf nu = 0.1 rho = 1 u_max = 1

def pde(x, y): u, v, p = y[:, 0:1], y[:, 1:2], y[:, 2:3] ux = jacobian(y, x, i=0, j=0) uy = jacobian(y, x, i=0, j=1) ut = jacobian(y, x, i=0, j=2) vx = jacobian(y, x, i=1, j=0) vy = jacobian(y, x, i=1, j=1) vt = jacobian(y, x, i=1, j=2) px = jacobian(y, x, i=2, j=0) py = jacobian(y, x, i=2, j=1) uxx = hessian(y, x, component=0, i=0, j=0) uyy = hessian(y, x, component=0, i=1, j=1) vxx = hessian(y, x, component=1, i=0, j=0) vyy = hessian(y, x, component=1, i=1, j=1) momentum1 = ut+uux+vuy+1/rhopx-nu(uxx+uyy) momentum2 = vt+uvx+vvy+1/rhopy-nu(vxx+vyy) continuity = ux+vy return [momentum1, momentum2, continuity]

lx = 5 ly = 2 geom1 = Rectangle(xmin=[0, 0], xmax=[lx, ly]) l = ly/5 xa = lx/4-l/2 xb = lx/4+l/2 ya = ly/2-l/2 yb = ly/2+l/2

geom2 = Rectangle(xmin=[xa, ya], xmax=[xb, yb]) geom = CSGDifference(geom1, geom2)

timedomain = TimeDomain(0, 2) geomtime = GeometryXTime(geom, timedomain)

def init_0(X):

return 0

def init_umax(X): return u_max

ic_u = IC(geomtime, initumax, lambda , on_initial: on_initial, component=0) ic_v = IC(geomtime, init0, lambda , on_initial: on_initial, component=1) ic_p = IC(geomtime, init0, lambda , on_initial: on_initial, component=2)

def boundary_column(X, on_bounday): return on_bounday and geom2.on_boundary((X[0], X[1]))

def boundary_sides(X, on_boundary): y = X[1] top = np.isclose(y, 2) bottom = np.isclose(y, 0) sides = top or bottom return on_boundary and sides

def boundary_inlet(X, on_boundary): x = X[0] y = X[1] inlet = np.isclose(x, 0) return on_boundary and inlet

def condition_0(X): return 0

def condition_umax(X): return u_max

bc_column_u = DirichletBC(geomtime, condition_0, boundary_column, component=0) bc_column_v = DirichletBC(geomtime, condition_0, boundary_column, component=1) bc_column_p = NeumannBC(geomtime, condition_0, boundary_column, component=2) bc_sides_p = NeumannBC(geomtime, condition_0, boundary_sides, component=2) bc_inlet_u = DirichletBC(geomtime, condition_umax, boundary_inlet, component=0) bc_inlet_v = DirichletBC(geomtime, condition_0, boundary_inlet, component=1)

data = TimePDE(geomtime, pde, [bc_column_u, bc_column_v, bc_column_p, bc_sides_p, bc_inlet_u, bc_inlet_v, ic_u, ic_v, ic_p], num_domain=3000, num_initial=1000, num_boundary=3000) net = maps.FNN([3, 50, 50, 50, 50, 50, 3], "tanh", "Glorot normal")

model = Model(data, net) model.compile("adam", lr=1e-4) losshistory, train_state = model.train(epochs=5000)

model.compile("L-BFGS") losshistory, train_state = model.train()

saveplot(train_state=train_state, loss_history=losshistory, isplot=True, issave=False)

dx = 0.01 dy = 0.01 dt = 0.01 x = np.arange(0, 5+dx, dx) y = np.arange(0, 2+dy, dy) ts = np.arange(0, 2+dt, dt)

X = np.zeros((len(x)len(y), 3)) xs = np.vstack((x,)len(y)).reshape(-1) ys = np.vstack((y,)*len(x)).T.reshape(-1) X[:, 0] = xs X[:, 1] = ys

X[:, 2] = 2 Y = model.predict(X) for i in range(len(X)): if X[i][0] >= xa and X[i][0] <= xb and X[i][1] >= ya and X[i][1] <= yb: Y[i] = 0

u = Y[:, 0].reshape(len(y), len(x)) v = Y[:, 1].reshape(len(y), len(x)) p = Y[:, 2].reshape(len(y), len(x))

K = 3 figure(figsize=(Kx[-1], Ky[-1])) streamplot(x, y, u, v) contourf(x, y, p) colorbar() savefig('pinn.png')

q769855234 commented 2 years ago

and the correction solution shoube be like this: image

q769855234 commented 2 years ago

but I just get the solution like this: image

praksharma commented 2 years ago

It is clear that you defined the wrong problem. Your inlet velocity doesn't even match. Please can you paste the problem statement instead of the simulation figure? Also don't forget to paste the loss values for each iteration.

q769855234 commented 2 years ago

1650758598(1)

q769855234 commented 2 years ago

很明显,您定义了错误的问题。你的入口速度甚至不匹配。请问您可以粘贴问题陈述而不是模拟图吗?也不要忘记粘贴每次迭代的损失值。

the problem is flow around the square column, initial condition:u,v,p=0 or u=1,v,p=0 at domain boundary condition:u,v,∂p/∂n=0 at sides, coloumn u=umax=1, v=0 at inlet

q769855234 commented 2 years ago

很明显,你定义了一个问题。你的入口速度甚至不匹配。请问你可以把问题陈述出来而不是模拟错误吗?

问题是围绕方柱流动, 初始条件:u,v,p=0 或 u=1,v,p=0 在域 边界条件:u,v,∂p/∂n=0 在边上,列 u =umax=1, v=0 在入口处

the result of u=1 and u=0 at domain is not same, but shoule be similar

lululxvi commented 2 years ago

Your problem is similar to this one https://github.com/lululxvi/deepxde/issues/634