search

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

1D simple transient diffusion PDE with 2 Neumann BCs #40

Closed swsyoon closed 4 years ago

swsyoon commented 4 years ago

Hi Lu,

Thank you very much for sharing this amazing tool. I am trying to learn step-by-step process to build PDE and solve it by reviewing examples.

I am trying to write 1D transient diffusion PDE that I have already coded with finite difference method in MATLAB.

dydt = d2ydx2,
Two BCs are : dydx(x=0) = 0 and dydx(x=1) = -2* y IC : y(t=0) = 1 at all the points.

I have the code here:

def pde(x, y):
    dy_x = tf.gradients(y, x)[0]
    dy_x, dy_t = dy_x[:, 0:1], dy_x[:, 1:]
    dy_xx = tf.gradients(dy_x, x)[0][:, 0:1]
    return dy_t-dy_xx

def boundary_l(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0)

def boundary_r(x, on_boundary):
    return on_boundary and np.isclose(x[0], 1)

geom = dde.geometry.Interval(0, 1)
bc_l = dde.NeumannBC(geom, lambda X: 0 * X, boundary_l)
bc_r = dde.NeumannBC(geom, lambda X: -2 * X, boundary_r)
timedomain = dde.geometry.TimeDomain(0, 1)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)
ic = dde.IC(
    geomtime, lambda x: np.full_like(x[:, 0:1], 1), lambda _, on_initial: on_initial
)
bc = [bc_l, bc_r]
data = dde.data.TimePDE(
    geomtime,
    pde,
    [bc, ic],
    num_domain=400,
    num_boundary=40,
    num_initial=100,
    num_test=10000,
)

layer_size = [2] + [32] * 3 + [1]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)

model = dde.Model(data, net)

model.compile("adam", lr=0.001, metrics=["l2 relative error"])
losshistory, train_state = model.train(epochs=50000)

dde.saveplot(losshistory, train_state, issave=True, isplot=True)

Currently I am getting error in 

data = dde.data.TimePDE(
    geomtime,
    pde,
    [bc, ic],
    num_domain=400,
    num_boundary=40,
    num_initial=100,
    num_test=10000,
)

AttributeError: 'list' object has no attribute 'collocation_points'

I think I did something wrong in bc and ic argument but I am not sure what exactly is wrong. Can you give me some comments on this ? Thank you.

lululxvi commented 4 years ago

Use the following:

data = dde.data.TimePDE(
    geomtime,
    pde,
    [bc_l, bc_r, ic],
    num_domain=400,
    num_boundary=40,
    num_initial=100,
    num_test=10000,
)
swsyoon commented 4 years ago

Thank you for your quick response. However I got another error at the same location. 

ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 2 and the array at index 1 has size 1

I am wondering whether my bc_r is right ? How to apply dydx = -2y in the code below ? ` bc_r = dde.NeumannBC(geom, lambda X: -2 X, boundary_r) `

lululxvi commented 4 years ago 
bc_l = dde.NeumannBC(geom, lambda X: np.zeros((len(X), 1)), boundary_l)
bc_r = dde.RobinBC(geom, lambda X, y: -2 * y, boundary_r)
swsyoon commented 4 years ago

Thanks. I modified it, but still get the same error above at the same location of PDE definition line

data = dde.data.TimePDE(
    geomtime,
    pde,
    [bc_l, bc_r, ic],
    num_domain=400,
    num_boundary=40,
    num_initial=100,
    num_test=10000,
)

I added more information about error message below. 

................................................
File "C:\Users\test\AppData\Local\Continuum\anaconda3\envs\tensorflow_env\lib\site-packages\deepxde-0.6.1-py3.7.egg\deepxde\data\pde.py", line 122, in bc_points
  File "<__array_function__ internals>", line 6, in vstack
File "C:\Users\test\AppData\Local\Continuum\anaconda3\envs\tensorflow_env\lib\site-packages\numpy\core\shape_base.py", line 283, in vstack
    return _nx.concatenate(arrs, 0)
  File "<__array_function__ internals>", line 6, in concatenate
ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 2 and the array at index 1 has size 1

I appreciate your help.

lululxvi commented 4 years ago

Use:

bc_l = dde.NeumannBC(geomtime, lambda X: np.zeros((len(X), 1)), boundary_l)
bc_r = dde.RobinBC(geomtime, lambda X, y: -2 * y, boundary_r)

model.compile("adam", lr=0.001)
swsyoon commented 4 years ago

Thank you very much! It works now.

LiHua1213 commented 3 years ago

def pde(x, u): u_vel, v_vel, p = u[:, 0:1], u[:, 1:2], u[:, 2:]............. spatial_domain = dde.geometry.Rectangle(xmin=[-0.5, -0.5], xmax=[1, 1.5]) timedomain = dde.geometry.TimeDomain(0, 1) geomtime = dde.geometry.GeometryXTime(spatial_domain, timedomain)

def boundary_outflow(x, on_boundary): return on_boundary and np.isclose(x[0], 1) def boundary_inlet(x, on_boundary): return on_boundary and np.isclose(x[0], -0.5) def boundary_b(x, on_boundary): return on_boundary and np.isclose(x[1], -0.5) def boundary_t(x, on_boundary): return on_boundary and np.isclose(x[1], 1.5) boundary_condition_in0 = dde.DirichletBC( geomtime, lambda x: np.ones((len(x), 1)) * 8, boundary_inlet, component=0 ) boundary_condition_in1 = dde.DirichletBC( geomtime, lambda x: np.zeros((len(x), 1)) , boundary_inlet, component=1 ) boundary_condition_bot0 = dde.DirichletBC( geomtime, lambda x: np.zeros((len(x), 1)), boundary_b, component=0 ) boundary_condition_bot1 = dde.DirichletBC( geomtime, lambda x: np.zeros((len(x), 1)), boundary_b, component=1 ) boundary_condition_top1 = dde.DirichletBC( geomtime, lambda x: np.zeros((len(x), 1)), boundary_t, component=2 )

@boundary_condition_top0 = dde.NeumannBC( @geomtime, lambda x: np.zeros((len(x), 1)), boundary_t, component=1 @) @boundary_condition_top2 = dde.NeumannBC( @geomtime, lambda x: np.zeros((len(x), 1)), boundary_t, component=0 @) boundary_condition_out2 = dde.DirichletBC( geomtime, lambda x: np.zeros((len(x), 1)), boundary_outflow, component=2 ) @boundary_condition_out0 = dde.NeumannBC( @geomtime, lambda x: np.zeros((len(x), 1)), boundary_outflow, component=0 @) @boundary_condition_out1 = dde.NeumannBC( @geomtime, lambda x: np.zeros((len(x), 1)), boundary_outflow, component=1 @) dear lu : when i have the boundary conditions (including @)above, an error is roperted as follows: ValueError: Rectangle: Method boundary_normal do not accept points on the vertexes. when i have the boundary conditions (without @)above, the code run correctly. So if i want to run correctly with the boundary conditions (including @)above, how can revise the boundary conditions? Thanks

lululxvi commented 3 years ago

See #303

LiHua1213 commented 3 years ago

Thanks for your help! Now,i want to solve a N-S(2D) problem, and the codes are as follows: I tried some methods(such as changing the num_domain, num_test; changing the network size; changing the learning rate and increasing the model train epochs; changing the loss weights ) from "Q: I failed to train the network or get the right solution, e.g., the training loss is large." at FAQ. However, the train loss and test loss are still large. Is there something wrong for my setting? Please give me some advice to solve such problem! Thank you!

import matplotlib.pyplot as plt import numpy as np import matplotlib.ticker as ticker import deepxde as dde from IPython.core.pylabtools import figsize from matplotlib.path import Path import matplotlib.patches as patches

def main(): Re = 20

def pde(x, u):
    u_vel, v_vel, p = u[:, 0:1], u[:, 1:2], u[:, 2:]
    u_vel_x = dde.grad.jacobian(u, x, i=0, j=0)
    u_vel_y = dde.grad.jacobian(u, x, i=0, j=1)
    u_vel_xx = dde.grad.hessian(u, x, component=0, i=0, j=0)
    u_vel_yy = dde.grad.hessian(u, x, component=0, i=1, j=1)

    v_vel_x = dde.grad.jacobian(u, x, i=1, j=0)
    v_vel_y = dde.grad.jacobian(u, x, i=1, j=1)
    v_vel_xx = dde.grad.hessian(u, x, component=1, i=0, j=0)
    v_vel_yy = dde.grad.hessian(u, x, component=1, i=1, j=1)

    p_x = dde.grad.jacobian(u, x, i=2, j=0)
    p_y = dde.grad.jacobian(u, x, i=2, j=1)

    momentum_x = (
        u_vel * u_vel_x + v_vel * u_vel_y + p_x - 1 / Re * (u_vel_xx + u_vel_yy)
    )
    momentum_y = (
        u_vel * v_vel_x + v_vel * v_vel_y + p_y - 1 / Re * (v_vel_xx + v_vel_yy)
    )
    continuity = u_vel_x + v_vel_y

    return [momentum_x, momentum_y, continuity]

spatial_domain1 = dde.geometry.Rectangle(xmin=[0, 0], xmax=[2, 2])
spatial_domain2 = dde.geometry.Rectangle(xmin=[1, 1], xmax=[1.5, 1.5])
spatial_domain = dde.geometry.CSGDifference(spatial_domain1, spatial_domain2)
timedomain = dde.geometry.TimeDomain(0, 1)
geomtime = dde.geometry.GeometryXTime(spatial_domain, timedomain)

def boundary_outflow(x, on_boundary):
    return on_boundary and spatial_domain1.on_boundary(x[:2]) and np.isclose(x[0], 2)
def boundary_inlet(x, on_boundary):
    return on_boundary and spatial_domain1.on_boundary(x[:2]) and np.isclose(x[0], 0)

def boundary_wall1(x, on_boundary):

return on_boundary and np.isclose(x[0], 1)

def boundary_wall2(x, on_boundary):

return on_boundary and np.isclose(x[0], 1.5)

def boundary_wall3(x, on_boundary):

return on_boundary and np.isclose(x[1], 1)

def boundary_wall4(x, on_boundary):

return on_boundary and np.isclose(x[1], 1.5)

def boundary_bottom(x, on_boundary):
    return on_boundary and spatial_domain1.on_boundary(x[:2]) and np.isclose(x[1], 0)
def boundary_top(x, on_boundary):
    return on_boundary and spatial_domain1.on_boundary(x[:2]) and np.isclose(x[1], 2)

def boundary_rect(x, on_boundary):
    return on_boundary and spatial_domain2.on_boundary(x[:2])

def func(x):
    return np.zeros((len(x), 1))

def func_g(x):
    return 8*np.ones((len(x), 1 ))

boundary_condition_in0 = dde.DirichletBC(geomtime, func_g, boundary_inlet, component=0)
boundary_condition_in1 = dde.DirichletBC(geomtime, func, boundary_inlet, component=1)
boundary_condition_bot1 = dde.DirichletBC(geomtime, func, boundary_bottom, component=1)

boundary_condition_bot2 = dde.DirichletBC(geomtime, func, boundary_bottom, component=0)

boundary_condition_bot11 = dde.NeumannBC(geomtime, func, boundary_bottom, component=0)
boundary_condition_bot12 = dde.NeumannBC(geomtime, func, boundary_bottom, component=1)
boundary_condition_top1 = dde.DirichletBC(geomtime, func, boundary_top, component=1)

boundary_condition_top2 = dde.DirichletBC(geomtime, func, boundary_top, component=0)

boundary_condition_top11 = dde.NeumannBC(geomtime, func, boundary_top, component=0)
boundary_condition_top12 = dde.NeumannBC(geomtime, func, boundary_top, component=1)

boundary_condition_wall11 = dde.DirichletBC(geomtime, func, boundary_wall1, component=0)

boundary_condition_wall12 = dde.DirichletBC(geomtime, func, boundary_wall1, component=1)

boundary_condition_wall21 = dde.DirichletBC(geomtime, func, boundary_wall2, component=0)

boundary_condition_wall22 = dde.DirichletBC(geomtime, func, boundary_wall2, component=1)

boundary_condition_wall31 = dde.DirichletBC(geomtime, func, boundary_wall3, component=0)

boundary_condition_wall32 = dde.DirichletBC(geomtime, func, boundary_wall3, component=1)

boundary_condition_wall41 = dde.DirichletBC(geomtime, func, boundary_wall4, component=0)

boundary_condition_wall42 = dde.DirichletBC(geomtime, func, boundary_wall4, component=1)

boundary_condition_out0 = dde.DirichletBC(geomtime, func, boundary_outflow, component=2)
boundary_condition_out1 = dde.NeumannBC(geomtime, func, boundary_outflow, component=0)
boundary_condition_out2 = dde.NeumannBC(geomtime, func, boundary_outflow, component=1)
boundary_condition_wall11 = dde.DirichletBC(geomtime, func, boundary_rect, component=0)
boundary_condition_wall12 = dde.DirichletBC(geomtime, func, boundary_rect, component=1)
ic_u = dde.IC(geomtime, lambda x: np.zeros((len(x), 1)), lambda _, on_initial: on_initial, component=0)
ic_v = dde.IC(geomtime, lambda x: np.zeros((len(x), 1)), lambda _, on_initial: on_initial, component=1)

data = dde.data.TimePDE(
    geomtime,
    pde,
    [boundary_condition_in0, boundary_condition_in1, boundary_condition_bot1, boundary_condition_top1, boundary_condition_out0, boundary_condition_wall11,
     boundary_condition_wall12, boundary_condition_top12, boundary_condition_bot12, boundary_condition_top11, boundary_condition_bot11,
     boundary_condition_out1, boundary_condition_out2, ic_u, ic_v],
    num_domain=4000,
    num_boundary=800,
    num_initial=300,
    num_test=4000
)

net = dde.maps.FNN([3] + 2 * [50] + [3], "tanh", "Glorot normal")

model = dde.Model(data, net)

model.compile("adam", lr=1e-4)
model.train(epochs=50000)
model.compile("L-BFGS-B", )
model.train()

------------------ 原始邮件 ------------------ 发件人: "lululxvi/deepxde" @.>; 发送时间: 2021年6月16日(星期三) 中午12:11 @.>; @.**@.>; 主题: Re: [lululxvi/deepxde] 1D simple transient diffusion PDE with 2 Neumann BCs (#40)

See #303

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

lululxvi commented 3 years ago

It is hard to see any apparent mistake in the code. Did you check this example https://github.com/lululxvi/deepxde/blob/master/examples/Beltrami_flow.py? One issue is that there are some many IC/BCs. Try to enforce hard constraints for some of them, see FAQ. If you still cannot solve this, you may open a new issue with more details.

LiHua1213 commented 3 years ago

Thanks for your reply. I checked the exemples you mentioned. Now the code can be run with a train loss below 10e-4 but the results was not good compared with numerical siulation. I am trying  to adjust the codes from Q: I failed to train the network or get the right solution, e.g., the training loss is large. However,  the results was also not good compared with numerical siulation. So i am trying to solve an inverse problem now. Thank you again, lu!

------------------ 原始邮件 ------------------ 发件人: "lululxvi/deepxde" @.>; 发送时间: 2021年6月19日(星期六) 上午6:03 @.>; @.**@.>; 主题: Re: [lululxvi/deepxde] 1D simple transient diffusion PDE with 2 Neumann BCs (#40)

It is hard to see any apparent mistake in the code. Did you check this example https://github.com/lululxvi/deepxde/blob/master/examples/Beltrami_flow.py? One issue is that there are some many IC/BCs. Try to enforce hard constraints for some of them, see FAQ. If you still cannot solve this, you may open a new issue with more details.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

LiHua1213 commented 3 years ago

Hello, lu.I want to use  double precision and i checked the FAQ#247. But i can not find the FAQ on how to use double precision.please give me some help, thanks!

------------------ 原始邮件 ------------------ 发件人: "lululxvi/deepxde" @.>; 发送时间: 2021年6月19日(星期六) 上午6:03 @.>; @.**@.>; 主题: Re: [lululxvi/deepxde] 1D simple transient diffusion PDE with 2 Neumann BCs (#40)

It is hard to see any apparent mistake in the code. Did you check this example https://github.com/lululxvi/deepxde/blob/master/examples/Beltrami_flow.py? One issue is that there are some many IC/BCs. Try to enforce hard constraints for some of them, see FAQ. If you still cannot solve this, you may open a new issue with more details.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

lululxvi commented 3 years ago

In the begining, use

dde.config.real.set_float64()