lululxvi / deepxde

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

Is this the correct way/syntax to write the problem ? #1253

Closed athtareq closed 1 year ago

athtareq commented 1 year ago

Hello, I'm trying to solve the problem

Screenshot from 2023-04-17 12-11-25

where $\Omega$ is a [0,2]x[0,1] rectangle and the solutions are $c_1$, $c_2$ and $\phi$

I wrote the following code and it's working, but the errors keep fluctuating between 1e-1 and 1e-3. Can someone please check if the syntax is correct ? Note that the $div$ operator is only comprised of space derivatives.

m1=1
m2=-1
d1=1e-3
d2=2*1e-5
eps=1e-4

def lap(u,x):
  u_xx = dde.grad.hessian(u, x, i=0, j=0)
  u_yy = dde.grad.hessian(u, x, i=1, j=1)
  lap=u_xx+u_yy
  return lap

def pde(x, y):
    c1, c2, phi = y[:, 0:1], y[:, 1:2],y[:, 2:]
    grad_phi = dde.grad.jacobian(phi, x, i=0) # gradient of phi
    #phi_y = dde.grad.jacobian(phi, x, i=0, j=1) # derivative of the the output wrt second input (y_coor)
    c1_t = dde.grad.jacobian(c1, x, i=0, j=2)
    c2_t = dde.grad.jacobian(c2, x, i=0, j=2) 
    div1_1 = dde.grad.jacobian(c1*grad_phi[: , 0:2], x, i=0, j=0) 
    div1_2 = dde.grad.jacobian(c1*grad_phi[: , 0:2], x, i=0, j=1) 
    div1=div1_1+div1_2
    div2_1 = dde.grad.jacobian(c2*grad_phi[: , 0:2], x, i=0, j=0) 
    div2_2 = dde.grad.jacobian(c2*grad_phi[: , 0:2], x, i=0, j=1)
    div2=div2_1+div2_2

    return [c1_t - d1*lap(c1,x)-m1*div1+c1*c2,
            c2_t - d2*lap(c2,x)-m2*div2-c1*c2,
            -eps*lap(phi,x)-c1+c2]

def boundary(x, on_boundary):
    return on_boundary

def func_zero(x) :
    return np.zeros([len(x),1])

geom = dde.geometry.Rectangle([0,0],[1,2])
timedomain = dde.geometry.TimeDomain(0, 1)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

bc1 = dde.icbc.NeumannBC(geomtime, func_zero, boundary,0)
bc2 = dde.icbc.NeumannBC(geomtime, func_zero, boundary,1)
bc3 = dde.icbc.DirichletBC(geomtime, func_zero, boundary,2)
ic1 = dde.icbc.IC(geomtime,lambda x: 5,lambda _, on_initial: on_initial,0)
ic2 = dde.icbc.IC(geomtime,lambda x: 10,lambda _, on_initial: on_initial,1)

data = dde.data.TimePDE(geomtime,pde,[bc1,bc2,bc3,ic1,ic2],num_domain=800,num_boundary=80,num_initial=80)

net = dde.nn.FNN([3] + [70] * 8 + [3], "tanh", "Glorot normal")

model = dde.Model(data, net)
model.compile("adam", lr=1e-3)
losshistory, train_state = model.train(iterations=150000)
dde.saveplot(losshistory, train_state, issave=True, isplot=True)
athtareq commented 1 year ago

@lululxvi sorry for the disturbance, any ideas ?

lululxvi commented 1 year ago

If the code can run, then the syntax should be correct.

athtareq commented 1 year ago

@lululxvi thank you!