The result of the time-dependent problem is incorrect after the extension of time

[123new-net]

@lululxvi Thank you for your contribution and the answer to my previous doubts. Now I have a new problem, I hope you can help me solve it. As mentioned in the question, the results calculated by my model at 20S are correct, but after extending to 30s, the loss can also reach 10^-4 magnitude, but the results vary greatly and are wrong. I have a very long time requirement, and I want to know how DeepXDE handles problems with a large time span. Here is my code. `"""Backend supported: tensorflow.compat.v1, tensorflow, pytorch""" import deepxde as dde import numpy as np import matplotlib.pyplot as plt from deepxde.backend import tf

a1=1.67e-4 a2=0.45 a3=5e-4 a4=0.45

rectanglar1 = dde.geometry.Rectangle(xmin=[0,0], xmax=[5,5]) rectanglar2 = dde.geometry.Rectangle(xmin=[2.49,3], xmax=[2.51,3.1]) circle = dde.geometry.Disk([2.5,3.5],0.5) geom1 = dde.geometry.CSGUnion(circle,rectanglar2) geom = dde.geometry.CSGDifference(rectanglar1,geom1) timedomain=dde.geometry.TimeDomain(0,30) geomtime=dde.geometry.GeometryXTime(geom,timedomain)

def pde(x, u): p,s = u[:, 0:1], u[:, 1:2]

dp_xx = dde.grad.hessian(u, x, component=0, i=0, j=0)
dp_yy = dde.grad.hessian(u, x, component=0, i=1, j=1)

ds_t = dde.grad.jacobian(u, x, i=1, j=2)

f_o=a1*(dp_xx+dp_yy)-a2*ds_t
f_w=a3*(dp_xx+dp_yy)+a4*ds_t

return [f_o,f_w]

def boundary_inlet(x, on_boundary): return on_boundary and (np.isclose(x[1],3) and x[0]>=2.49 and x[0]<=2.51) def boundary_outlet(x, on_boundary): return on_boundary and np.isclose(x[1],0) def boundary_wall(x, on_boundary): return on_boundary and circle.on_boundary(x[0:2]) def func_0(x): return np.zeros((len(x),1), dtype=np.float32) def func_1(x): return np.ones((len(x),1), dtype=np.float32) def modify_output(X, u): x, y, t = X[:, 0:1], X[:, 1:2], X[:, 2:3] p, s = u[:, 0:1], u[:, 1:2] p_new = pt s_new = st return tf.concat((p_new, s_new), axis=1)

bc_1 = dde.icbc.DirichletBC(geomtime,lambda x:3,boundary_inlet,component=0) bc_2 = dde.icbc.DirichletBC(geomtime,func_1,boundary_inlet,component=1) bc_3 = dde.icbc.DirichletBC(geomtime,func_0,boundary_outlet,component=0) bc_4 = dde.icbc.NeumannBC(geomtime,lambda x:0,boundary_wall,component=0) bc_5 = dde.icbc.NeumannBC(geomtime,lambda x:0,boundary_wall,component=1) ic1= dde.icbc.IC(geomtime,lambda x:0, lambda ,on_initial: on_initial,component=0) ic2= dde.icbc.IC(geomtime,lambda x:0, lambda ,on_initial: on_initial,component=1) bcs = [bc_1, bc_2, bc_3, bc_4, bc_5, ic_1, ic_2]

data = dde.data.TimePDE(geomtime,pde,bcs,num_domain=5000,num_boundary=500,num_initial=500,num_test=10000)

plt.figure(figsize = (32, 18)) plt.scatter(data.train_x_all[:,0], data.train_x_all[:,1], s = 3) plt.axis('equal')

net = dde.nn.FNN([3] + 6 * [60] + [2], "tanh", "Glorot normal") net.apply_output_transform(modify_output)

model = dde.Model(data, net) resampler = dde.callbacks.PDEResidualResampler(period=100) model.compile("adam",lr=1e-3) model.train(epochs=30000) model.compile("L-BFGS") losshistory, train_state = model.train()

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

Thanks!

[praksharma]

TLDR. If your loss is converging and you are getting a wrong solution to the PDE then it means that either the problem is ill-posed i.e. there are multiple solutions or you defined some other problem. Please tell us the problem you are solving and if possible, plz paste the true and predicted solution at certain time steps.

Also your code is not readable. Plz put it in this container.

"""Backend supported: tensorflow.compat.v1, tensorflow, pytorch"""
import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt
from deepxde.backend import tf

a1=1.67e-4
a2=0.45
a3=5e-4
a4=0.45

rectanglar1 = dde.geometry.Rectangle(xmin=[0,0], xmax=[5,5])
rectanglar2 = dde.geometry.Rectangle(xmin=[2.49,3], xmax=[2.51,3.1])
circle = dde.geometry.Disk([2.5,3.5],0.5)
geom1 = dde.geometry.CSGUnion(circle,rectanglar2)
geom = dde.geometry.CSGDifference(rectanglar1,geom1)
timedomain=dde.geometry.TimeDomain(0,30)
geomtime=dde.geometry.GeometryXTime(geom,timedomain)

def pde(x, u):
    p,s = u[:, 0:1], u[:, 1:2]

    dp_xx = dde.grad.hessian(u, x, component=0, i=0, j=0)
    dp_yy = dde.grad.hessian(u, x, component=0, i=1, j=1)

    ds_t = dde.grad.jacobian(u, x, i=1, j=2)

    f_o=a1*(dp_xx+dp_yy)-a2*ds_t
    f_w=a3*(dp_xx+dp_yy)+a4*ds_t

    return [f_o,f_w]
def boundary_inlet(x, on_boundary):
    return on_boundary and (np.isclose(x[1],3) and x[0]>=2.49 and x[0]<=2.51)
def boundary_outlet(x, on_boundary):
    return on_boundary and np.isclose(x[1],0)
def boundary_wall(x, on_boundary):
    return on_boundary and circle.on_boundary(x[0:2])
def func_0(x):
    return np.zeros((len(x),1), dtype=np.float32)
def func_1(x):
    return np.ones((len(x),1), dtype=np.float32)
def modify_output(X, u):
    x, y, t = X[:, 0:1], X[:, 1:2], X[:, 2:3]
    p, s = u[:, 0:1], u[:, 1:2]
    p_new = pt
    s_new = st
    return tf.concat((p_new, s_new), axis=1)

bc_1 = dde.icbc.DirichletBC(geomtime,lambda x:3,boundary_inlet,component=0)
bc_2 = dde.icbc.DirichletBC(geomtime,func_1,boundary_inlet,component=1)
bc_3 = dde.icbc.DirichletBC(geomtime,func_0,boundary_outlet,component=0)
bc_4 = dde.icbc.NeumannBC(geomtime,lambda x:0,boundary_wall,component=0)
bc_5 = dde.icbc.NeumannBC(geomtime,lambda x:0,boundary_wall,component=1)
ic_1= dde.icbc.IC(geomtime,lambda x:0, lambda _,on_initial: on_initial,component=0)
ic_2= dde.icbc.IC(geomtime,lambda x:0, lambda _,on_initial: on_initial,component=1)
bcs = [bc_1, bc_2, bc_3, bc_4, bc_5, ic_1, ic_2]

data = dde.data.TimePDE(geomtime,pde,bcs,num_domain=5000,num_boundary=500,num_initial=500,num_test=10000)

plt.figure(figsize = (32, 18))
plt.scatter(data.train_x_all[:,0], data.train_x_all[:,1], s = 3)
plt.axis('equal')

net = dde.nn.FNN([3] + 6 * [60] + [2], "tanh", "Glorot normal")
net.apply_output_transform(modify_output)

model = dde.Model(data, net)
resampler = dde.callbacks.PDEResidualResampler(period=100)
model.compile("adam",lr=1e-3)
model.train(epochs=30000)
model.compile("L-BFGS")
losshistory, train_state = model.train()

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

[123new-net]

@praksharma Hi. Thank you very much for your reply. What I'm dealing with is the problem of diffusion in porous media, where the fluid diffuses with pressure from a rectangle at the bottom of a circle. And there seems to be no solution to this problem as far as I know. When you said maybe I defined some other problem, Do you mean that the definitions of BC and IC are wrong? Also, I used Ctrl + E for my code and I'm not sure why it became unreadable. I'm sorry for that.

[praksharma]

What do you mean when you say " And there seems to be no solution to this problem as far as I know"? The problem has no solution when you solve it using FEM or FVM?

I think your problem is same as the rotating cone problem a transient benchmark problem in CFD.

For this "When you said maybe I defined some other problem, Do you mean that the definitions of BC and IC are wrong?" yes.

You can use the following style to paste your multi-line python code with indentation.

[123new-net]

@praksharma Thank you for your prompt reply. Maybe I misunderstood what you were saying about the solution. When I say there is no solution, I mean there is no analytical solution. I have obtained numerical solutions for this problem using FVM. If you mean that you want me to import numerical solutions into the model, then that works for me. And I was wondering if you have any code or literature on the rotating cone problem for the transient benchmark problem in CFD. Thanks in advance.

[praksharma]

This is the original paper. You can find a Python code by searching for citations on Google scholar. Original_paper.pdf

I can see your code. It would be better if you write the problem so we can verify if the problem was correctly implemented.

[123new-net]

@praksharma Thanks for sharing. Now I am again confused about the solution. Adding real solution data to DeepXDE should only be used to calculate residuals and relative errors, as in the Burgers Equation case. In addition, in question #155 , Dr. Lulu also mentioned that the reference solutions in the examples are only used to check the accuracy of the solution obtained from DeepXDE. This is not consistent with my understanding of usage. I think a loss can be added to calculate the error between the predicted solution and the real solution, so as to constrain the predicted solution to be close to the real solution.

[praksharma]

Ok. Let us start from basics. When You solve a PDE using FEM, you do not supply the true solution. Because if you already know the solution why would you even solve the PDE in the first place. The same is true for PINNs. PINNs do not require true solution.

I requested you to paste the true solution (from FEM or FVM) here so that we can see the difficulty of the of the problem. DeepXDE uses the true solution, if supplied, to calculate the test loss which is a performance metric and has nothing to do with the training of PINNs. Hope this clarifies your understanding.

We can see into your problem if you could paste the PDE, the domain and the boundary conditions of your problem.

[123new-net]

@praksharma Thank you for your patient explanation. Now I understand what you mean. The problem I'm trying to solve is the diffusion of liquid through the soil after leaking from the bottom of the underground pipe. The PDE is shown below. The boundary conditions are that the pressure at the hole is 3 and the saturation of the liquid is 1. The initial conditions are that the pressure and saturation are 0. The size of the region is not fixed, because I think it gets bigger over time. But the outer boundary of the region is always maintained to be a zero pressure boundary. The solution at time t=100 that I obtained using FEM is shown below.

Also, I read some papers on PINN and some mentioned that a Loss could be added to calculate the error between the predicted value and the true value, and I wondered if DeepXDE could do this.

[praksharma]

Will look into your problem. Why do you need two rectangles rectanglar2 = dde.geometry.Rectangle(xmin=[2.49,3], xmax=[2.51,3.1])?

Yes in some PINN frameworks people use the true solution during the training. We refer to them as the reconstruction error. But than the PINN is no more a PDE solver it becomes a data-driven PDE solver.

[praksharma]

You can apply the reconstruction error using PointSetBC. Here is the source code.

[123new-net]

@praksharma Thank you for your prompt reply! The second rectangle is the one I used to represent the simplified hole. What if my goal is to build a Data-driven PDE Solver? I wonder if DeepXDE supports this kind of operation? I will carefully study the module PointSetBC, thank you for your advice.

[praksharma]

Yes it is possible and that why I told you about PointSetBC. This is how you use it

reconstruction_constraint = dde.icbc.PointSetBC(true_x, true_y)

where true_x and true_y are simply numpy arrays of true solution at particular spatio-temporal locations.

[123new-net]

@praksharma Thank you very much for your help. Now my model is working very well. And when I say it works well, I mean my two-dimensional model. Now I want to extend to the 3D model. I have a problem when I set up the geometry. I hope you can help me out. I tried to build a cylinder in geometry, but I couldn't find the command to build a cylinder in DeepXDE's API. It seems that in DeepXDE.Geometry the 3D geometry is only cuboid and sphere. Thanks again for your help!

[123new-net]

@lululxvi Or is there no source code for the cylinder in the DeepXDE library?

[praksharma]

If the geometry is not defined, you can pass anchor points for collocation points and PointSetBC for BC points.

[123new-net]

@praksharma Thank you for your reply! You mean I can just go to data = dde.data.timepde(...) add anchors for training, without setting geometry?

[praksharma]

yes.

[123new-net]

@praksharma Thank you for your help. This scheme is feasible. I can easily implement time-independent problems when I use the scheme you provide. For example, Poission_Lshape. But I had trouble with time-dependent problems. Because in the anchor you need a column of time data. For example, in the Burgers equation problem, I use the following code to generate anchors. The results worked well. x,t = np.meshgrid(np.linspace(-1,1, 50), np.linspace(0,0.99, 50)) anchors = np.vstack((np.ravel(x), np.ravel(t))).T But this approach to my own problem, even in a two-dimensional model, when you don't take too much x, the number of training points combined with t increases dramatically. Not to mention scaling up to a 3D model. I wonder if you have any advice on that. In addition, I want to know for general problems where geometric regions can be set, what is the relationship between x and t in the set of training points? My understanding is that different x correspond to different t, all randomly distributed in the region.

[123new-net]

By the way, what should geomtime change in the following code if I do not set the geometry region? bc = dde.icbc.DirichletBC(geomtime, ...) data = dde.data.TimePDE(geomtime, ... anchors=X)

[praksharma]

Yes the number of points will increase exponentially with the increase in number of features. Instead of np.meshgrid() I would suggest you to use itertools.product().

For the second comment, I guess you can simply leave the geometry and pde. You need to pass the anchor points and the BCs. I didn't try this. May be you can look into the source code. Also for the initial condition provide a higher weight coefficient.

[123new-net]

@praksharma Thank you very much for your advice. But after my attempt, it seems that it can not solve my problem very well. I can't build a cylinder in geometry, so I can't set the boundary conditions on the surface of the cylinder. So even if I enter any points in anchors that I want, they're not going to be recognized as any of the points in the border that I need to satisfy BC. You can imagine that for Poission_Lshape problems, the points in my input anchors are L-shaped. But for the geometric region, I can only build a full rectangle. The result is shown below.

I know this is a very tricky problem, and thank you again for your help. @lululxvi I know this is presumptuous, but I wonder if you could add the source code for the cylinder.

[praksharma]

I mean just create a geometry in your favourite CAD tool. Mesh the geometry and extract the nodal coordinates. Firstly, you need to put all the nodal coordinates in the anchor points. Then extract the boundary coordinates (there are many ways) for each individual prescribed value and put it in PointSetBC .

I have done this for a 100-times more complicated geometry. It works. You can write your own python code. You will probably need itertools.product.

[123new-net]

@praksharma Thank you for your reply. I can understand the first half of your proposal. But I don't quite understand what you mean when you say that extract boundary coordinates for each individual prescribed value. I was wondering if you could give me an example, for example Poission_Lshape.

I learned the use of itertools.product, which is the Cartesian product of input traversable objects. Since all the coordinates I input are data, doesn't this function function just like meshgrid?

[praksharma]

Suppose you have 2 different BCs. You need to extract the nodal coordinates for each BC and manually apply the prescribed value for each BC. You can use np.meshgrid but the code looks cleaner with itertools.product.

[123new-net]

@praksharma Thank you for your patient explanation. Now I think I understand what you mean. You are saying that all BCS and ICS need to be converted to PointSetBC for representation. As shown below, true_x represents the coordinates I entered and true_y represents the boundary values I need to set. I don't know if I understand correctly. bc=dde.icbc.PointSetBC(true_x, true_y) However, this seems to only implement the DirichletBC modification, which is not easy for NeumannBC to implement.

[praksharma]

Yes, that is correct.

True this is not implemented for NeumannBC.

[123new-net]

@praksharma Thanks for your advice, now my 3D model works as well. However, the results are still somewhat lacking, and I would like to ask some more questions. It's about the loss weights. When I do not set the loss weights, the loss of a BC in the result is large, so I set the weight of this BC to a smaller value. Although the total loss has decreased, the result does not seem to meet the requirements of small weight BC. My understanding is that when I put a small weight on the BC, the network will weaken the ability to satisfy the BC. I wonder if you have any suggestions.

[praksharma]

Can you please put a 0.5 weight for the PDE loss and 1.0 for all the BC losses?

Although the total loss has decreased, the result does not seem to meet the requirements of small weight BC.

I have also experienced the same.

[123new-net]

I tried your suggestion, but the loss of BC is still large.

[praksharma]

How much epochs did you train?

[123new-net]

I trained 30000 with the adam optimizer, and then L-BFGS.

[praksharma]

can you show the expected solution and the predicted solution?

[123new-net]

I used PointSetBC to substitute solutions at several points in time. The predicted results at these time points are correct. But at other time points, the results can not meet the other boundary conditions well.

[praksharma]

This happens. Thats why we have separate PINNs for temporal problems.It is an active field of research.

[yang1805524385]

Do I need to define callbacks in deepxde？