zxw4688 commented 3 years ago

Hello @lululxvi ， I would like to solve the Wave euqation as ： In my code, the initial value of the Wave equation is defined as
My code is : ` from future import absolute_import from future import division from future import print_function

import numpy as np

import deepxde as dde from deepxde.backend import tf

def main(): beta=4.0 c=(beta**2)/8

定义pde的形式，默认右端项为零

def pde(x, y):
    dy_x = tf.gradients(y, x)[0]
    dy_x, dy_y,dy_t = dy_x[:, 0:1], dy_x[:, 1:2], dy_x[:, 2:3]
    dy_xx = tf.gradients(dy_x, x)[0][:, 0:1]
    dy_yy = tf.gradients(dy_y, x)[0][:, 1:2]
    dy_tt = tf.gradients(dy_t, x)[0][:, 2:3]
    return dy_tt - c*(dy_xx+ dy_yy)

def func(x):
    return np.sin(x[:, 0:1]+x[:, 1:2])*np.cos(x[:, 0:1]+x[:, 1:2])*np.sin(beta*x[:,2:3])
def dt_func(x):
    return beta*np.sin(x[:, 0:1]+x[:, 1:2])*np.cos(x[:, 0:1]+x[:, 1:2])*np.cos(beta*x[:,2:3])

确定区域的形状：矩形区域，时间为[0,1]

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

给出区域的初始条件

IBC0 = dde.IC(geomtime, func, lambda _, on_initial: on_initial)
dt_IBC0 = dde.IC(geomtime, dt_func, lambda _, on_initial: on_initial)
DBC = dde.DirichletBC(geomtime, func, lambda _, on_boundary: on_boundary)

data = dde.data.TimePDE(
    geomtime,
    pde,
    [DBC, IBC0, dt_IBC0],
    num_domain=1000,
    num_boundary=10,
    num_initial=30,
    solution=func,
    num_test=1000,
)

layer_size = [3] + [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.0001, metrics=["l2 relative error"])
losshistory, train_state = model.train(epochs=50000)

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

X = geomtime.random_points(100000)
y_true = func(X)
y_pred = model.predict(X)
print("L2 relative error:", dde.metrics.l2_relative_error(y_true, y_pred))
np.savetxt("test0921.dat", np.hstack((X, y_true, y_pred)))

if name == "main": main()`

zxw4688 commented 3 years ago

And the result of NN as follow : Initializing variables... Training model...

Step Train loss 0 [1.73e-01, 1.93e-01, 6.99e-02, 2.95e+00] 1000 [8.44e-03, 8.57e-02, 4.36e-01, 8.60e-01] 2000 [7.03e-03, 6.50e-02, 4.74e-01, 7.95e-01] 3000 [8.54e-03, 3.05e-02, 4.97e-01, 7.47e-01] 4000 [2.65e-03, 2.85e-02, 5.14e-01, 7.07e-01] 5000 [3.09e-03, 2.68e-02, 5.34e-01, 6.55e-01] 6000 [1.75e-03, 1.42e-02, 5.63e-01, 6.13e-01] 7000 [1.12e-03, 7.77e-03, 5.73e-01, 6.02e-01] 8000 [1.16e-03, 4.75e-03, 5.76e-01, 5.98e-01] 9000 [1.10e-03, 3.00e-03, 5.78e-01, 5.95e-01] 10000 [8.10e-04, 2.07e-03, 5.80e-01, 5.93e-01] 11000 [5.92e-04, 1.46e-03, 5.81e-01, 5.92e-01] 12000 [4.78e-04, 1.08e-03, 5.82e-01, 5.91e-01] 13000 [4.18e-04, 8.56e-04, 5.83e-01, 5.90e-01] 14000 [3.70e-04, 6.94e-04, 5.84e-01, 5.89e-01] 15000 [3.34e-04, 5.54e-04, 5.84e-01, 5.88e-01] 16000 [3.02e-04, 4.32e-04, 5.84e-01, 5.88e-01] 17000 [2.75e-04, 3.33e-04, 5.84e-01, 5.88e-01] 18000 [2.52e-04, 2.46e-04, 5.84e-01, 5.88e-01] 19000 [2.33e-04, 1.85e-04, 5.84e-01, 5.87e-01] 20000 [2.14e-04, 1.36e-04, 5.85e-01, 5.87e-01] 21000 [1.97e-04, 9.94e-05, 5.85e-01, 5.87e-01] 22000 [1.79e-04, 7.19e-05, 5.85e-01, 5.87e-01] 23000 [1.65e-04, 5.28e-05, 5.85e-01, 5.87e-01] 24000 [1.50e-04, 3.82e-05, 5.85e-01, 5.87e-01] 25000 [1.40e-04, 2.77e-05, 5.85e-01, 5.87e-01] 26000 [1.30e-04, 2.07e-05, 5.85e-01, 5.87e-01] 27000 [1.18e-04, 1.80e-05, 5.86e-01, 5.86e-01] 28000 [1.10e-04, 1.34e-05, 5.86e-01, 5.86e-01] 29000 [1.03e-04, 9.94e-06, 5.85e-01, 5.86e-01] 30000 [1.03e-04, 8.42e-06, 5.86e-01, 5.86e-01] 31000 [9.26e-05, 6.88e-06, 5.86e-01, 5.86e-01] 32000 [8.70e-05, 5.63e-06, 5.86e-01, 5.86e-01] 33000 [8.28e-05, 5.39e-06, 5.86e-01, 5.86e-01] 34000 [7.89e-05, 3.72e-06, 5.85e-01, 5.87e-01] 35000 [7.54e-05, 3.65e-06, 5.86e-01, 5.86e-01] 36000 [7.19e-05, 3.43e-06, 5.86e-01, 5.86e-01] 37000 [6.88e-05, 2.86e-06, 5.86e-01, 5.86e-01] 38000 [6.60e-05, 2.55e-06, 5.86e-01, 5.86e-01] 39000 [8.24e-05, 2.37e-06, 5.85e-01, 5.87e-01] 40000 [6.11e-05, 2.07e-06, 5.86e-01, 5.86e-01] 41000 [5.89e-05, 1.89e-06, 5.86e-01, 5.86e-01] 42000 [5.69e-05, 1.59e-06, 5.85e-01, 5.86e-01] 43000 [5.49e-05, 1.58e-06, 5.86e-01, 5.86e-01] 44000 [5.31e-05, 1.50e-06, 5.86e-01, 5.86e-01] 45000 [5.14e-05, 1.36e-06, 5.86e-01, 5.86e-01] 46000 [4.98e-05, 1.23e-06, 5.86e-01, 5.86e-01] 47000 [5.07e-05, 1.66e-06, 5.86e-01, 5.86e-01] 48000 [4.69e-05, 1.13e-06, 5.86e-01, 5.86e-01] 49000 [4.55e-05, 1.01e-06, 5.86e-01, 5.86e-01] 50000 [4.55e-05, 1.30e-06, 5.84e-01, 5.87e-01] Test loss Test metric
[1.76e-01, 0.00e+00, 0.00e+00, 0.00e+00] [1.61e+00]
[1.13e-02, 0.00e+00, 0.00e+00, 0.00e+00] [1.28e+00]
[1.20e-02, 0.00e+00, 0.00e+00, 0.00e+00] [1.26e+00]
[1.76e-02, 0.00e+00, 0.00e+00, 0.00e+00] [1.20e+00]
[6.44e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.20e+00]
[1.06e-02, 0.00e+00, 0.00e+00, 0.00e+00] [1.24e+00]
[6.20e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.34e+00]
[3.81e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.37e+00]
[4.95e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.42e+00]
[5.18e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.48e+00]
[4.28e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.54e+00]
[3.73e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.61e+00]
[3.56e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.69e+00]
[3.44e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.75e+00]
[3.36e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.81e+00]
[3.23e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.86e+00]
[3.06e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.91e+00]
[2.86e-03, 0.00e+00, 0.00e+00, 0.00e+00] [1.96e+00]
[2.66e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.01e+00]
[2.46e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.06e+00]
[2.28e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.11e+00]
[2.10e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.15e+00]
[1.95e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.19e+00]
[1.82e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.22e+00]
[1.68e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.25e+00]
[1.58e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.28e+00]
[1.49e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.30e+00]
[1.39e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.31e+00]
[1.31e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.33e+00]
[1.25e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.34e+00]
[1.18e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.35e+00]
[1.13e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.36e+00]
[1.09e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.37e+00]
[1.04e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.38e+00]
[1.00e-03, 0.00e+00, 0.00e+00, 0.00e+00] [2.39e+00]
[9.65e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.39e+00]
[9.34e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.40e+00]
[9.02e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.40e+00]
[8.74e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.41e+00]
[8.44e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.41e+00]
[8.23e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.41e+00]
[8.00e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.42e+00]
[7.78e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.42e+00]
[7.58e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.42e+00]
[7.37e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.43e+00]
[7.19e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.43e+00]
[7.02e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.43e+00]
[6.93e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.43e+00]
[6.69e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.44e+00]
[6.54e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.44e+00]
[6.39e-04, 0.00e+00, 0.00e+00, 0.00e+00] [2.44e+00]

Best model at step 49000: train loss: 1.17e+00 test loss: 6.54e-04 test metric: [2.44e+00]

'train' took 410.831801 s

Saving loss history to loss.dat ... Saving training data to train.dat ... Saving test data to test.dat ... Predicting... 'predict' took 0.187535 s

L2 relative error: 2.43178430

zxw4688 commented 3 years ago

From the result, it is not good!, I doubt that the definition of initial value is not right in my code, but I can not know how to revise it. Please help me to revise it, Thank you very much!

zxw4688 commented 3 years ago

The figure of result is :

smao-astro commented 3 years ago

dt_IBC0 = dde.IC(geomtime, dtfunc, lambda , on_initial: on_initial)

I do not think IC implemented first order derivative initial condition, what do you think of it? @lululxvi

zxw4688 commented 3 years ago

Hello @smao-astro Then, how to revise it ? Thank you very much！

zxw4688 commented 3 years ago

Hello @smao-astro And I agree with your view! but I do not know how to do it.

zxw4688 commented 3 years ago

Hello @lululxvi , on the above the problem, how do I revise it ?, thank you very much!

lululxvi commented 3 years ago

IC means Dirichlet IC. In your case, you can assume it is Neumann BC, and define a new on_boundary. see #134

zxw4688 commented 3 years ago

@lululxvi 你好，对于波动方程，我的真解是u=sin(x+y)sin(betat), c=beta^2/2, beta=4. 麻烦你看一下，现在已经按照说的改过了，可以运行，但是依然是一阶时间导数的初值条件对应的Loss 不会下降，我认为还是这个一阶时间导数的初始条件给错了，麻烦LuLu博士有空的时候看一下，具体的代码如下：

`from future import absolute_import from future import division from future import print_function

import numpy as np

import deepxde as dde from deepxde.backend import tf

def main(): beta=np.square(2) c=(beta**2)/2

定义pde的形式，默认右端项为零

def pde(x, y):
    dy_x = tf.gradients(y, x)[0]
    dy_x, dy_y,dy_t = dy_x[:, 0:1], dy_x[:, 1:2], dy_x[:, 2:3]
    dy_xx = tf.gradients(dy_x, x)[0][:, 0:1]
    dy_yy = tf.gradients(dy_y, x)[0][:, 1:2]
    dy_tt = tf.gradients(dy_t, x)[0][:, 2:3]
    return dy_tt - c*(dy_xx+ dy_yy)

 #给出真解函数:u=sin(x+y)*sin(beta*t)
def func(x):
    return np.sin(x[:, 0:1]+x[:, 1:2])*np.sin(beta*x[:,2:3])
def NDfunc(x):
    return np.cos(x[:, 0:1]+x[:, 1:2])*np.sin(beta*x[:,2:3])

#求出真解关于时间的一阶导数
def dt_func(x):
    return beta*np.sin(x[:, 0:1]+x[:, 1:2])*np.cos(beta*x[:,2:3])
#def dt_func(x):
   # return    tf.gradients(func(x),x)[0][:,2:3]

def boundary_initial(x, _):
    return  np.isclose(x[-1], 0 )

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

确定区域的形状：矩形区域，时间为[0,1]

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

给出区域的初始条件

IBC0 = dde.IC(geomtime, func, lambda _, on_initial: on_initial)
dt_IBC0 = dde.NeumannBC(geomtime, dt_func, lambda _, on_initial: boundary_initial)
DBC = dde.DirichletBC(geomtime, func, lambda _, on_boundary: on_boundary)

data = dde.data.TimePDE(
    geomtime,
    pde,
    [DBC, IBC0, dt_IBC0],
    num_domain=1000,
    num_boundary=100,
    num_initial=100,
    solution=func,
    num_test=1000,

    )

layer_size = [3] + [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=3000)

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

X = geomtime.random_points(100000)
y_true = func(X)
y_pred = model.predict(X)
print("L2 relative error:", dde.metrics.l2_relative_error(y_true, y_pred))
np.savetxt("test0921.dat", np.hstack((X, y_true, y_pred)))

if name == "main": main() `

lululxvi commented 3 years ago

Use dt_IBC0 = dde.NeumannBC(geomtime, dt_func, boundary_initial). You have the same errors in your other codes.

zxw4688 commented 3 years ago

非常感谢，这种方法，我试过了，可以运行，但是一阶时间导数的初值训练误差没有减少，或者几乎没有下降，都是大于1 的，这个不对啊！

zxw4688 commented 3 years ago

而且我认为这句代码和原代码dt_IBC0 = dde.NeumannBC(geomtime, dtfunc, lambda , on_initial: boundary_initial)应该表达的时同一个意思。

zxw4688 commented 3 years ago

其他的代码我也换了为：计算结果：依然是一阶时间导数的初值条件对应的Loss 不会下降，我认为还是这个一阶时间导数的初始条件给错了。

具体代码为： `from future import absolute_import from future import division from future import print_function

import numpy as np

import deepxde as dde from deepxde.backend import tf

def main(): beta=np.square(2) c=(beta**2)/2

定义pde的形式，默认右端项为零

def pde(x, y):
    dy_x = tf.gradients(y, x)[0]
    dy_x, dy_y,dy_t = dy_x[:, 0:1], dy_x[:, 1:2], dy_x[:, 2:3]
    dy_xx = tf.gradients(dy_x, x)[0][:, 0:1]
    dy_yy = tf.gradients(dy_y, x)[0][:, 1:2]
    dy_tt = tf.gradients(dy_t, x)[0][:, 2:3]
    return dy_tt - c*(dy_xx+ dy_yy)

 #给出真解函数:u=sin(x+y)*sin(beta*t)
def func(x):
    return np.sin(x[:, 0:1]+x[:, 1:2])*np.sin(beta*x[:,2:3])

#求出真解关于时间的一阶导数
def dt_func(x):
    return beta*np.sin(x[:, 0:1]+x[:, 1:2])*np.cos(beta*x[:,2:3])

def boundary_initial(x, _):
    return  np.isclose(x[-1], 0 )

确定区域的形状：矩形区域，时间为[0,1]

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

给出区域的初始条件

IBC0 = dde.IC(geomtime, func, lambda _, on_initial: on_initial)
dt_IBC0 = dde.NeumannBC(geomtime, dt_func, boundary_initial)
DBC = dde.DirichletBC(geomtime, func, lambda _, on_boundary: on_boundary)
#DBC_L= dde.NeumannBC(geomtime, NDfunc, boundary_l)

def data_input_2D_T(nx,ny,nt,Lx,Rx,Ly,Ry,Lt,Rt):
    xx=np.linspace(Lx,Rx,nx)
    yy=np.linspace(Ly,Ry,ny)
    tt=np.linspace(Lt,Rt,nt)
    X,Y=np.meshgrid(xx,yy)
    X=np.reshape(X,[nx*ny,1])
    Y=np.reshape(Y,[nx*ny,1])
    XX,T=np.meshgrid(X,tt)
    YY,T=np.meshgrid(Y,tt)
    X_x=np.reshape(XX,[nx*ny*nt,1])
    Y_y=np.reshape(YY,[nx*ny*nt,1])
    T_t=np.reshape(T,[nx*ny*nt,1])
    XYT_xyt=np.concatenate((X_x,Y_y,T_t),1)
    #print(XT,np.shape(XT))
    return X_x,Y_y,T_t,XYT_xyt
X_x,Y_y,T_t,XYT_xyt=data_input_2D_T(50,50,5,0,1,0,1,0,1)

data = dde.data.TimePDE(
    geomtime,
    pde,
    [DBC, IBC0, dt_IBC0],
    anchors=XYT_xyt,
    solution=func,

    )

layer_size = [3] + [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=3000)

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

X = geomtime.random_points(100000)
y_true = func(X)
y_pred = model.predict(X)
print("L2 relative error:", dde.metrics.l2_relative_error(y_true, y_pred))
np.savetxt("test0921.dat", np.hstack((X, y_true, y_pred)))

if name == "main": main() `

lululxvi commented 3 years ago

No. lambda _, on_initial: boundary_initial and boundary_initial are totally different.
Check the example Euler_beam.py where there is first-order, second-order and third order BC.
It may be good to start from 1D case to test.

zxw4688 commented 3 years ago

你的意思是将“lambda _, on_initial:”这个去掉，改为：

zxw4688 commented 3 years ago

非常感谢LuLu博士的回答，但是我不是很清楚，这个语句： DBC_u_F = dde.DirichletBC(geom, funtrue, lambda , on_boundary: boundary_Front) 和这个语句： DBC_u_F = dde.DirichletBC(geom, fun_true, boundary_Front) 有什么区别？

lululxvi commented 3 years ago

Read https://www.w3schools.com/python/python_lambda.asp for lamba expression.

zxw4688 commented 3 years ago

@lululxvi 非常感谢LuLu博士的准确、耐心、细致的回答!,但是这里lambda像表述的就是一个匿名函数的而已，可是我没有看出DBC_u_F = dde.DirichletBC(geom, funtrue, lambda , on_boundary: boundary_Front) 和这个语句： DBC_u_F = dde.DirichletBC(geom, fun_true, boundary_Front)的区别。实在不好意思啊。才学python不久。非常感谢。

lululxvi commented 3 years ago

You may try the following:

def f(x):
    return x ** 2

g = lambda x: f

f(1)
g(1)

zxw4688 commented 3 years ago

@lululxvi 博士，你好！非常感谢你提供的软件包，我分别尝试了流体方程和固体方程，均可以得到相对比较不错的解。在这里非常感谢！但是我现在想求解一个简单的界面问题，但是第一步如何表示方程都困难，因为方程是在不同区域上的，不知道你的deepxde能否求解该问题，具体问题如下：。也就是如何书写函数pde(x,y). 非常感谢！

lululxvi commented 3 years ago

Yes, it is doable. Because you don't have a global PDE, then you need to set the pde as None , see https://deepxde.readthedocs.io/en/latest/modules/deepxde.data.html#deepxde.data.pde.PDE. Next, you use two OperatorBC to define the two PDEs in two domains. Note, OperatorBC not only works for BC, since you can define the on_boundary as the whole domain of Omage_1 or Omage_2.

zxw4688 commented 3 years ago

非常感谢，LuLu博士的快速准确的回答，我进去看了看，不是很明白，有没有相应的例子啊，非常感谢!

lululxvi commented 3 years ago

I don't have online examples. In OperatorBC(geom, func, on_boundary), you define the pde in func and the corresponding domain in on_boundary.

zxw4688 commented 3 years ago

哦，明白了，我试试，非常感谢lulu博士

	朱兴文

邮箱：zxw4688@126.com |

签名由网易邮箱大师定制

On 10/21/2020 17:07, Lu Lu wrote:

I don't have online examples. In OperatorBC(geom, func, on_boundary), you define the pde in func and the corresponding domain in on_boundary.

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

zxw4688 commented 3 years ago

@lululxvi 博士，你好! 我按照你的建议尝试着写了一个计算椭圆界面问题：的代码，但是出错了，麻烦lulu博士有空的时候帮我看看，主要的问题就是不同区域上不同方程的定义，应该如何描述，以及真解应该如何表示。非常感谢！具体错误为：代码的图片为：具体的代码为： `from future import absolute_import from future import division from future import print_function

import numpy as np

from tensorflow.keras.backend import set_floatx

import sys sys.path.insert( 0, "E:/DNN_matlab/myrebp/solve_PDE_0811/DeepXDE/deepxde-master")

import deepxde as dde from deepxde.geometry.csg import CSGDifference from deepxde.geometry.csg import CSGUnion

import tensorflow as tf

def main(): set_floatx("float64") beta1=1 beta2=2 def pde_u1(x,y): u1 = y[:, 0:1] u_x = tf.gradients( u1 , x )[0] du1_x, du1_y = u_x[:, 0:1], u_x[:, 1:2 ] beta1_du1_xx = tf.gradients( beta1du1_x, x )[0][:, 0:1 ] beta1_du1_yy = tf.gradients( beta1du1_y, x )[0][:, 1:2 ]
f1=-4beta1(x[:,0:1]2+x[:,1:2]2+1)tf.exp(x[:,0:1]2+x[:,1:2]2)
return - beta1_du1_xx - beta1_du1_yy - f1 def pde_u2(x,y):
u2 = y[:, 0:1 ] u_x = tf.gradients( u2 , x )[0] du2_x, du2_y = u_x[:, 0:1], u_x[:, 1:2 ]
beta2_du2_xx = tf.gradients( beta2du2_x, x )[0][:, 0:1 ] beta2_du2_yy = tf.gradients( beta2du2_y, x )[0][:, 1:2 ]
f2=-2beta2tf.exp(x[:,0:1]+x[:,1:2])
return - beta2_du2_xx - beta2_du2_yy - f2
def interface_equ_1(x,y): u1 = y[:, 0:1] u2 = y[:, 0:1] h1 = tf.exp(x[:,0:1]2+x[:,1:2]2) - tf.exp(x[:,0:1]+x[:,1:2]) return u1 - u2 - h1
def interface_equ_2(x,y): u1 = y[:, 0:1] u2 = y[:, 0:1] u1_x = tf.gradients( u1 , x )[0]
du1_x, du1_y = u1_x[:, 0:1], u1_x[:, 1:2 ]
u2_x = tf.gradients( u2 , x )[0]
du2_x, du2_y = u2_x[:, 0:1], u2_x[:, 1:2 ]
h2 =beta2tf.exp(x[:,0:1]+x[:,1:2]) (x[:,0:1]+x[:,1:2]-2)- 2beta1tf.exp(x[:,0:1]2+x[:,1:2]2) (x[:,0:1]2+x[:,1:2]2-(x[:,0:1]+x[:,1:2]))
return (1-x[:,0:1])(beta1du1_x-beta2du2_x)+(1-x[:,1:2])(beta1du1_y-beta2du2_y) - h2 circle = dde.geometry.Disk([1, 1], 0.5) rectangular = dde.geometry.Rectangle([0, 0], [2, 2]) geom1 = CSGDifference(rectangular, circle) geom2 = circle geom = CSGUnion(geom1,geom2)
def true_fun(x): """ y1 = exp(x^2+y^2) notes u1 y2 = exp(x+y) notes u2 """ if x in geom1: return np.exp(x[:,0:1]2+x[:,1:2]2) if x in geom2: return np.exp(x[:,0:1]+x[:,1:2])
def fun_u1(x): return np.exp(x[:,0:1]2+x[:,1:2]2) def fun_u2(x): return np.exp(x[:,0:1]+x[:,1:2]) def boundarygeom1(, on_boundary): return on_boundary and geom1
def boundarygeom2(, on_boundary): return on_boundary and geom2 def boundary_rect(x, on_boundary): return on_boundary and rectangular.on_boundary(x)
def boundary_circle(x, on_boundary): return on_boundary and circle.onboundary(x)
pde1 = dde.OperatorBC(geom, lambda x, y, : pde_u1(x,y), boundarygeom1 ) pde2 = dde.OperatorBC(geom, lambda x, y, : pde_u2(x,y), boundary_geom2 ) DBC_u1 = dde.DirichletBC(geom, funu1, lambda , on_boundary: boundary_rect ) InterBCu12 = dde.OperatorBC(geom, lambda x, y, : interface_equ_1(x,y), boundary_circle ) N_InterBCu12 = dde.OperatorBC(geom, lambda x, y, : interface_equ_2(x,y), boundary_circle ) data = dde.data.PDE( geom, None, [pde1,pde2,DBC_u1,InterBC_u12,N_InterBC_u12], num_boundary=10000,

solution=fun_u1,

                 num_test=10000
                 )
layer_size = [2] + [50] * 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)
losshistory, train_state = model.train(epochs=50000)
dde.saveplot(losshistory, train_state, issave=True, isplot=True)
X = geom.random_points(10000)
y_true = true_fun(X)
y_pred = model.predict(X)
print("L2 relative error:", dde.metrics.l2_relative_error(y_true, y_pred))
np.savetxt("test0921.dat", np.hstack((X, y_true, y_pred)))

if name == "main": main() `

lululxvi commented 3 years ago

There are a few issues in your code:

Line 62, geom is just the rectangle, right?
In true_fun, x in geom1 is not correct. First, x is a Nxd array; second, use geom1.inside(x). Check what it returns.
boundary_geom1 is wrong. I guess you want to use return geom1.inside(x).
Line 89, specify the points inside the domain, not only on boundary

Point 3 and 4 are more about Python grammar. You may start from simple examples to debug, e.g., 1D boundary value problem, and split the domain into left domain and right domain.

zxw4688 commented 3 years ago

@lululxvi 博士，非常感谢你的回答，我按照你说的改过了还是不能运行，麻烦LuLu博士有空的时候帮我看看，谢谢！具体错误为：具体代码为：界面问题的方程为：我的代码为： `from future import absolute_import from future import division from future import print_function

import numpy as np

from tensorflow.keras.backend import set_floatx

import sys sys.path.insert( 0, "E:/DNN_matlab/myrebp/solve_PDE_0811/DeepXDE/deepxde-master")

import deepxde as dde from deepxde.geometry.csg import CSGDifference from deepxde.geometry.csg import CSGUnion

import tensorflow as tf

def main(): set_floatx("float64") beta1=1 beta2=2 def pde_u1(x,y): u1 = y[:, 0:1] u_x = tf.gradients( u1 , x )[0] du1_x, du1_y = u_x[:, 0:1], u_x[:, 1:2 ] beta1_du1_xx = tf.gradients( beta1du1_x, x )[0][:, 0:1 ] beta1_du1_yy = tf.gradients( beta1du1_y, x )[0][:, 1:2 ]
f1=-4beta1(x[:,0:1]2+x[:,1:2]2+1)tf.exp(x[:,0:1]2+x[:,1:2]2)
return - beta1_du1_xx - beta1_du1_yy - f1 def pde_u2(x,y):
u2 = y[:, 0:1 ] u_x = tf.gradients( u2 , x )[0] du2_x, du2_y = u_x[:, 0:1], u_x[:, 1:2 ]
beta2_du2_xx = tf.gradients( beta2du2_x, x )[0][:, 0:1 ] beta2_du2_yy = tf.gradients( beta2du2_y, x )[0][:, 1:2 ]
f2=-2beta2tf.exp(x[:,0:1]+x[:,1:2])
return - beta2_du2_xx - beta2_du2_yy - f2
def interface_equ_1(x,y): u1 = y[:, 0:1] u2 = y[:, 0:1] h1 = tf.exp(x[:,0:1]2+x[:,1:2]2) - tf.exp(x[:,0:1]+x[:,1:2]) return u1 - u2 - h1
def interface_equ_2(x,y): u1 = y[:, 0:1] u2 = y[:, 0:1] u1_x = tf.gradients( u1 , x )[0]
du1_x, du1_y = u1_x[:, 0:1], u1_x[:, 1:2 ]
u2_x = tf.gradients( u2 , x )[0]
du2_x, du2_y = u2_x[:, 0:1], u2_x[:, 1:2 ]
h2 =beta2tf.exp(x[:,0:1]+x[:,1:2]) *(x[:,0:1]+x[:,1:2]-2)

2beta1tf.exp(x[:,0:1]2+x[:,1:2]2) (x[:,0:1]2+x[:,1:2]2-(x[:,0:1]+x[:,1:2]))
return (1-x[:,0:1])(beta1du1_x-beta2du2_x)+(1-x[:,1:2])(beta1du1_y-beta2du2_y) - h2 circle = dde.geometry.Disk([1, 1], 0.5) rectangular = dde.geometry.Rectangle([0, 0], [2, 2]) geom1 = CSGDifference(rectangular, circle) geom2 = circle geom = CSGUnion(geom1,geom2)
def true_fun(x): """ y1 = exp(x^2+y^2) notes u1 y2 = exp(x+y) notes u2 """ if geom1.inside(x): return np.exp(x[:,0:1]2+x[:,1:2]2) if geom2.inside(x): return np.exp(x[:,0:1]+x[:,1:2])
def fun_u1(x): return np.exp(x[:,0:1]2+x[:,1:2]2) def fun_u2(x): return np.exp(x[:,0:1]+x[:,1:2]) def boundary_geom1(x, on_boundary): return on_boundary and geom1.inside(x)
def boundary_geom2(x, on_boundary): return on_boundary and geom2.inside(x) def boundary_rect(x, on_boundary): return on_boundary and rectangular.on_boundary(x)
def boundary_circle(x, on_boundary): return on_boundary and circle.onboundary(x)
pde1 = dde.OperatorBC(geom, lambda x, y, : pde_u1(x,y), boundarygeom1 ) pde2 = dde.OperatorBC(geom, lambda x, y, : pde_u2(x,y), boundary_geom2 ) DBC_u1 = dde.DirichletBC(geom, funu1, lambda , on_boundary: boundary_rect ) InterBCu12 = dde.OperatorBC(geom, lambda x, y, : interface_equ_1(x,y), boundary_circle ) N_InterBCu12 = dde.OperatorBC(geom, lambda x, y, : interface_equ_2(x,y), boundary_circle ) data = dde.data.PDE( geom, None, [pde1,pde2,DBC_u1,InterBC_u12,N_InterBC_u12], num_domain=10000, num_boundary=1000, solution=true_fun, num_test=10000 ) layer_size = [2] + [50] 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) losshistory, train_state = model.train(epochs=50000) dde.saveplot(losshistory, train_state, issave=True, isplot=True) X = geom.random_points(10000) y_true = true_fun(X) y_pred = model.predict(X) print("L2 relative error:", dde.metrics.l2_relative_error(y_true, y_pred)) np.savetxt("test0921.dat", np.hstack((X, y_true, y_pred)))
if name == "main": main() `

lululxvi commented 3 years ago

Install v0.8.6

lululxvi / deepxde

How to define the initial value of Wave equation in your package Deepxde? #127

定义pde的形式，默认右端项为零

确定区域的形状：矩形区域，时间为[0,1]

给出区域的初始条件

定义pde的形式，默认右端项为零

确定区域的形状：矩形区域，时间为[0,1]

给出区域的初始条件

定义pde的形式，默认右端项为零

确定区域的形状：矩形区域，时间为[0,1]

给出区域的初始条件

solution=fun_u1,