Defining BCs

[akapet00]

Less an issue, more noobish kind of question:

I can not understand how to define boundary conditions for, lets say, a simple analytically solvable ODE. If the equation in its general form looks like:
k y_xx - a y + q = 0

with BCs:
(Dirichlet) T(x=L) = T_c
(Robins) -k T_x(x=0) = c[d - T(x=0)]

how can I define this using DDE?

Thank you in advance!

[lululxvi]

Check the example here https://github.com/lululxvi/deepxde/blob/master/examples/Poisson_Robin_1d.py

[akapet00]

Thank you. In my case the Dirichlet b.c. is not a function of some variable x. It is constant value. I tried something like:

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], L)

geom = dde.geometry.Interval(0, L)
bc_l = dde.RobinBC(geom, lambda x, y: (h_0/k)*y - (h_0*T_f/k), boundary_l)
bc_r = dde.DirichletBC(geom, lambda x: T_c, boundary_r)
data = dde.data.PDE(geom, 1, pde, [bc_l, bc_r], 16, 2, func=None, num_test=100)

but I get the error saying that:

AttributeError: 'int' object has no attribute 'shape'

When I created the array with all values set to zero and only the last element to T_c, I got completely wrong prediction plot. I guess this is a wrong approach?

[engsbk]

Thank you. In my case the Dirichlet b.c. is not a function of some variable x. It is constant value. I tried something like:

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], L)

geom = dde.geometry.Interval(0, L)
bc_l = dde.RobinBC(geom, lambda x, y: (h_0/k)*y - (h_0*T_f/k), boundary_l)
bc_r = dde.DirichletBC(geom, lambda x: T_c, boundary_r)
data = dde.data.PDE(geom, 1, pde, [bc_l, bc_r], 16, 2, func=None, num_test=100)

but I get the error saying that:

AttributeError: 'int' object has no attribute 'shape'

When I created the array with all values set to zero and only the last element to T_c, I got completely wrong prediction plot. I guess this is a wrong approach?

Is there a difference in writing the BC when dealing with geom (time independent) and geomtime (time dependent) domains.

I get different results when I try:

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

bc_l = dde.DirichletBC(geomtime, lambda x: 0, boundary_l)

and

bc_l = dde.DirichletBC(geomtime, lambda x: 0, lambda _, boundary_l: boundary_l)

Which one is the correct way for implementing the left boundary condition to be zero?

[lululxvi]

Could you try:
bc_r = dde.DirichletBC(geom, lambda x: 0 * x + T_c, boundary_r)

[akapet00]

Unfortunately, no luck again.
This is what I'm getting:

when the actual analytical solution is:

Besides that, the error is occurring:

/home/alk/anaconda3/lib/python3.7/site-packages/deepxde/metrics.py:13: RuntimeWarning: divide by zero encountered in double_scalars
return np.linalg.norm(y_true - y_pred) / np.linalg.norm(y_true)

This is, I believe, when the test metrics is calculated. Final result for test metrics is set to inf. Therefore, the test metrics plot is missing:

I presume that boundary conditions for this case are ill defined somehow... Do you have any documentation on how to define boundary conditions (Dirichlet, Neumann or Robins)?

The whole code I am trying to make it work is below :)

def pde(x, y):
    dy_x = tf.gradients(y, x)[0]
    dy_xx = tf.gradients(dy_x, x)[0]
    return k*dy_xx - w_b * y + (w_b*T_a + Q_m)

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], L)

geom = dde.geometry.Interval(0, L)
bc_l = dde.RobinBC(geom, lambda X, y: (h_0/k)*y - (h_0*T_f/k), boundary_l)
bc_r = dde.DirichletBC(geom, lambda x: 0 * x + T_c, boundary_r)
data = dde.data.PDE(geom, 1, pde, [bc_l, bc_r], 16, 2, func=None, num_test=100)

layer_size = [1] + [20] * 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=10000)

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

[lululxvi]

Could me provide the full code, including the values of w_b, Q_m, etc?

[akapet00]

Yes, of course. I copied and pasted relevant parts alongside with the function that calculates the analytical solution. I hope I did it correctly.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from numpy import sinh, cosh
import tensorflow as tf
tf.compat.v1.logging.set_verbosity(tf.compat.v1.logging.ERROR)
import deepxde as dde

# consts 
c_b = 4200 
rho_b = 1e3 
om_b = 5e-4
c = 4200 
h_f = 100
T_c = 37
L = 3e-2
k = 0.5
w_b = c_b*rho_b*om_b
Q_m = 33800
h_0 = 10
T_f = 25
T_a = 37

def analyticSol(x):
    A = np.sqrt(w_b/k)
    lterm = (T_c - T_a - Q_m/w_b) * (A * cosh(A * x)  +\
                +  (h_0/k) * sinh(A * x)) / (A * cosh(A * L) + (h_0/k) * sinh(A * L))
    rterm = h_0/k * (T_f - T_a - Q_m/w_b) * sinh(A * (L * np.ones(shape=x.shape) - x)) / ( A * cosh(A * L) + (h_0/k) * sinh(A * L)) 

    return  T_a + Q_m/w_b + lterm + rterm

def pde(x, y):
    dy_x = tf.gradients(y, x)[0]
    dy_xx = tf.gradients(dy_x, x)[0]
    return k*dy_xx - w_b*y + (w_b*T_a + Q_m)

def boundary_l(x, on_boundary):
    return on_boundary

def boundary_r(x, on_boundary):
    return on_boundary

def main():
    geom = dde.geometry.Interval(0, L)
    bc_l = dde.RobinBC(geom, lambda X, y: (h_0/k)*y - (h_0*T_f/k), boundary_l)
    bc_r = dde.DirichletBC(geom, lambda x: 0*x + T_c, boundary_r)
    data = dde.data.PDE(geom, 1, pde, [bc_l, bc_r], 16, 2, func=None, num_test=100)

    layer_size = [1] + [40] * 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=10000, uncertainty=True)

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

if __name__=='__main__':
    main()

[athtareq]

Hello, I have the same question regarding boundary conditions, and I get the same constant function as a result. Did you find any solutions or workarounds ?

[akapet00]

@tareqath I did not yet.
But the predicted function you see on the first figure is not actually of constant value throughout the domain. It looks like this:

Best,
Ante

[akapet00]

@tareqath @lululxvi Quick update. With the following code I am able to reproduce the overall change in the ODE but still I think that boundary conditions are not well defined:

def pde(x, y):
    dy_x = tf.gradients(y, x)[0]
    dy_xx = tf.gradients(dy_x, x)[0]
    return k*dy_xx - w_b*y
def boundary_l(x, on_boundary):
    return on_boundary
def boundary_r(x, on_boundary):
    return on_boundary
def func(x):
    return 0*x - (w_b*T_a + Q_m)

geom = dde.geometry.Interval(0, L)
bc_l = dde.RobinBC(geom, lambda x, y: (h_0/k)*y - (h_0*T_f/k), boundary_l)
bc_r = dde.DirichletBC(geom, lambda x: 0*x + T_c, boundary_r)
data = dde.data.PDE(geom, 1, pde, [bc_l, bc_r], 16, 2, func=func, num_test=100)

layer_size = [1] + [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, metrics=["l2 relative error"])
losshistory, train_state = model.train(epochs=10000, uncertainty=True)

What I am getting is:

But true results are set to 0 and T values are completely wrong.

Note that the problem is:

where the values of all constants are written here.

[akapet00]

Hi @lululxvi,

sorry for bothering you again. I tried this code for a few other PDE time-invariant equations (including those in your examples directory) and it worked flawlessly. I was going through the code base and I couldn't find what I am doing wrong. I presume it has to be something with initial conditions but even when I set both boundaries to Dirichlet it still produces same invalid results.
I am sending you current script (it is in .txt format, so should just convert it to .py when (and if :) ) you download it. First plot is prediction vs. analytical solution, second plot is showing learning curves.

Thank you!
Best,
Ante

[lululxvi]

@antelk I checked your code. Here are some modifications:

1. Define left and right boundary locations.

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], L)

Note: bc_r = dde.DirichletBC(geom, lambda x: np.full_like(x, T_c), boundary_r)

2. Use more training points, e.g., 100.

data = dde.data.PDE(geom, 1, pde, [bc_l, bc_r], 100, 2, func=analyticSol, num_test=100)

3. Because the solution is of order 10, then let's multiply the network output by 10.

net = dde.maps.FNN(...)
net.outputs_modify(lambda x, y: y * 10)

4. The three loss terms are of different scales when step = 0. Use loss_weights for rescaling, i.e., weighted losses.

model.compile("adam", lr=0.001, metrics=["l2 relative error"], loss_weights=[1e-7, 1e-2, 1])

5. Train with more iterations, e.g., 50000.

losshistory, train_state = model.train(epochs=50000)

Summary: By using the codes above, we can have the solution with the correct right boundary condition, but the solution is still of large error. The difficulty comes from the optimization. You can carefully monitor the loss values during the optimization. Your domain and PDE have different scales. Can you rescale your problem, such that L = 1 and most coefficients are of scale 1?

• By the way, i tried to use Dirichlet BCs for both left and right, and then it works perfectly.

bc_l = dde.DirichletBC(geom, lambda x: analyticSol(x), boundary_l)

[jprince58]

This strikes me as a perfect example was to why PDEs should be made dimensionless before trying to use DeepXDE for solutions. Dimensionless PDEs tend to be better scaled in general.

[akapet00]

Thank you! I am getting perfect results even with less iterations than 10k now. I completely ignored the fact that term w_b*T_a + Q_m is huge and that values on x and y axis are of different scale.
Can I ask you, is the uncertainty somehow connected with dropout as regularization procedure or something else? Do you think that there is some interesting way (other than Monte Carlo method) to propagate uncertainty through this neural network by treating inputs as variables from some statistical distribution?

[lululxvi]

Yes, the uncertainty is based on dropout, see https://www.sciencedirect.com/science/article/pii/S0021999119305340 for details. Treating inputs as distribution is doable, and there are some work on this, but I am not an expert on this.

[akapet00]

Alright. Thank you again.
I will close an issue now.

[LanPeng-94]

Thank you! I am getting perfect results even with less iterations than 10k now. I completely ignored the fact that term w_b*T_a + Q_m is huge and that values on x and y axis are of different scale.
Can I ask you, is the uncertainty somehow connected with dropout as regularization procedure or something else? Do you think that there is some interesting way (other than Monte Carlo method) to propagate uncertainty through this neural network by treating inputs as variables from some statistical distribution?

Hello，I tried to reproduce your case using DEEPXDE, I performed the calculation of the domain and the scaling of the PDE, I was able to calculate it accurately using the Dirichlet boundary condition, unfortunately, when I used the Robin boundary condition you provided I was unable to calculate it, I would like to ask you, how did you implement the Robin boundary condition calculation