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

Convective boundary condition problem #253

Closed thintn222 closed 2 years ago

thintn222 commented 3 years ago

Dear @lululxvi, Thank you for the great and useful library that can help a beginner like me approaching the PINNs techniques. Recently, I start investigating a heat transfer problem that consists of 2 convective boundary conditions (RobinBC) and 1 Dirichlet BC. Here is the domain image The PDE and RobinBC are shown, respectively, as image

First, to test the algorithm, I use only one convection BC (Robin one) one the upper edge CD. The result shows a good prediction with the temperature increases linearly from the bottom to the upper. This really makes sense. However, when I do the real problem with 2 convective BC (the right and upper edges), the results look extremely off as shown below, which the bottom edge seems not to hold the temperature during the simulation even I use Dirichlet in there.

image

Here is the code for implementing the algorithm

k = 52
h = 750
xmin = 0
xmax = 0.6
ymin = 0
ymax = 1
t_min = 0
t_max = 1
T_0 = 0
T_AB = 1
rhoCp = 36 

def pde(x, T):
    dT_x = tf.gradients(T, x)[0]
    dT_x, dT_y, dT_t = dT_x[:,0:1], dT_x[:,1:2], dT_x[:,2:]
    dT_xx = tf.gradients(dT_x, x)[0][:, 0:1]
    dT_yy = tf.gradients(dT_y, x)[0][:, 1:2]
    return k*dT_xx + k*dT_yy - rhoCp*dT_t
###
def boundary_l(x, on_boundary):
    return on_boundary and np.isclose(x[0], xmin)
def boundary_r(x, on_boundary):
    return on_boundary and np.isclose(x[0], xmax)
def boundary_b(x, on_boundary):
    return on_boundary and np.isclose(x[1], ymin)
def boundary_u(x, on_boundary):
    return on_boundary and np.isclose(x[1], ymax)

def func_d_AD(x):
    return np.ones((len(x),1), dtype=np.float32)*T_0

def func_n(x, T):
    return np.ones((len(x),1), dtype=np.float32)*(h/k)*(T-T_0) # convection BC

def func_d_AB(x):
    return np.ones((len(x),1), dtype=np.float32)*T_AB

geom = dde.geometry.Rectangle([xmin, ymin], [xmax, ymax])
timedomain = dde.geometry.TimeDomain(t_min, t_max)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)
##

bc_r = dde.RobinBC(geomtime, func_n, boundary_r)
bc_b = dde.DirichletBC(geomtime, func_d_AB, boundary_b) #func_d_AB, boundary_b)
bc_u = dde.RobinBC(geomtime, func_n, boundary_u)

##
data = dde.data.TimePDE(
    geomtime, 
    pde, 
    [bc_r, bc_b, bc_u], 
    num_domain=5000, 
    num_boundary=5000, 
    num_test= 1000
)

layer_size = [3] + [20] * 8 + [1]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
#net.apply_output_transform(lambda x, T:  (x[:, 1:2] - 1) * x[:, 1:2] * T + T_AB) #Is it right to define the hard constraints for dirichlet BC at the bottom edge?

##
model = dde.Model(data, net)
model.compile("adam", lr=0.001)
model.train(epochs=50000)
model.compile("L-BFGS-B")
losshistory, train_state = model.train()
dde.saveplot(losshistory, train_state, issave=True, isplot=True)

The training and testing are not converged well even I use the loss_weights option. I read all of the Q&A to get some ideas however the problem still existed. One more notice that I tried to use the model to predict the material point in the bottom edges only and see that they are not always equal to the reference value as the Dirichlet should do. Do you think the hard constraints BC will work in this problem?. Is it the right expression for the DirichletBC hard constraints in my case? net.apply_output_transform(lambda x, T: x[:, 1:2] * T + T_AB

Once again, thank you so much for the useful library DeepXDE. Do you have any suggestions that I can try out?

Thank you very much.

thintn222 commented 3 years ago

For updating, I'm pretty sure the error is about the hard constraints BC since the PDE does not satisfy the dirichletBC even I put lots of weight in there. How can I set the Dirichlet BC for this case? For example, T(x, y=0) = 1. I read the Q&A. However, these problems are about, i.e., y(0) = 0, not 1. Thank you, Lu Lu, again for the amazing library.

lululxvi commented 3 years ago

Yes, the hard constraint above is correct for T(x,y=0)=1.

thintn222 commented 3 years ago

Hi @lululxvi, I'm sorry for bothering you again. As your suggestion, since the PINN does not satisfy the BC, I use the hard constraint as net.apply_output_transform(lambda x, T: x[:, 1:2] * T + T_AB #T_AB = 1 However, the results are still completely off, and the RobinBC loss does not converge even the PDE loss converge well to 1e-4

image

I tried using loss weights to put more weight in the RobinBC, but it still cannot converge. All the data and geometry also have been scaled to 0-1. Do you have any suggestions that I can try out? Thank you so much for an amazing library <3

lululxvi commented 3 years ago

To verify your code, you can also use the exact solution as the training data, see the examples of inverse problems. If the loss still cannot do down, then it is highly possible that there might be some issue in the code.

Another suggestion is that you may first try a 1D problem.