how to write non-homogeneous Neumann boundary condition with out calculating by hand in 2D and interface boundary condition

[Karankaren]

Hi LuLu, Thank you so much for sharing all this great tool and for your quick response to all the people on Github. I am trying to implement a PINN to solve a PDE in a 2D rectangular domain using the domain decomposition method. After I read all the FAQ sections, I still haven't been able to write a couple of things correctly,

The normal derivate in Neumann part of boundary i.e. (du/dn = g where the reference function u(x,y)=x^2+y^2 without calculate derivative by hand like you did in (Poisson equation in 1D with Dirichlet/Neumann boundary conditions)). I mean, I want to use some thing like this ( first define derivative ( dy_x = tf.gradients(y, x)[0] dy_x, dy_y = dy_x[:, 0:1], dy_x[:, 1:2]) with respect normal vector and second define (def func(x): return x2+y2) then find some thing like this func_n(x): return (+or - 1)dy_x(func) ) The interface conditions u1=u2 and du_1/dn_1=du_2/dn_2. I hope you will answer me, thank.

[lululxvi]

Do you mean you would like to do something like Eq. 10 in https://doi.org/10.1364/OE.384875?

[Karankaren]

Thanks LuLu for your quick reply. Yes, but with normal derivative on interface of rectangular domain (xmin=[0, 0], xmax=[1, 1]) like this du_1/dn_1=-du_2/dn_2 with (n_1 and n_2 are normal direction on interface between subdomain 1 and subdomain 2 ) and also in interface u_1=u_2 but with line in interface between the two subdomains not exactly that you had in this paper.

subdomain2 subdomain1

---

u_2	u _1		interface	interface		n_ 2	n_ 1

._____. with u_1=u_2 on interface and du_1/dn_1=-du_2/dn_2 on interface

[lululxvi]

You can use OperatorBC to define it.

[Blackrobin15]

Hello,@Karankaren Did you solve this problem? I am facing the same problem now.

[lululxvi]

Related to https://github.com/lululxvi/deepxde/issues/421

[juliosdutra]

Do you mean you would like to do something like Eq. 10 in https://doi.org/10.1364/OE.384875?

Hi, Dr. Lu Lu,

Firstly, I would like to thank you for the remarkable contributions you have made to the community.

I am particularly interested in implementing the interface conditions to represent perfect thermal contact and heat flux continuity at a interface between to different solids, which is kind similar to Eq. 10 of your paper. Therefore, I would like to know if the code for this paper is available online. If you could provide any guidance or information, I would greatly appreciate it.

Problem statement

So far, I've considered a simple heat transfer problem, as if it were a rod composed of two different materials, with Dirichlet conditions at x=0 and x=L:

• from x=0 to x=L/2, we have conductivity coefficient a_1 and temperature T_1(x,t);
• from x=L/2 to x=L, we have coefficient a_2 and temperature T_2(x,t);
• at the interface, we have to meet both T_1 - T_2 = 0 and a_1dT_1/dx + a_2dT_2/dx = 0 for x=L/2 and t>0.

Since I did not find a concise solution I tried to implement this problem. But, I have some questions:

1. I found a way to enforce the heat flux continuity using geom.boundary_normal( ). But, the respective loss is not decreasing (in fact, it is constant). Is it correct to use geom.boundary_normal( )? Is there an alternative way?

2. I am struggling to come up with away to enforce T_1 - T_2 = 0. How could I get the values of T_1 and T_2 at x=L/2 during the model train?

Learning from the discussions I found here, my code is as follows:

#  Libraries to import
import deepxde as dde
import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt

# pamaters and geometries
a1 = 1.14 # conductivity coefficient 1
a2 = 0.01 # conductivity coefficient 2 # fairly different from a1
L = 1.0 # total length
T0 = 0.0 # temperature specified at x=0
TL = 1.0 # temperature specified at x=L
tend = 1.0 # final time for the simulation

geom1 = dde.geometry.Interval(0, L/2) # first solid
geom2 = dde.geometry.Interval(L/2, L) # second solid

timedomain = dde.geometry.TimeDomain(0, tend)
geomtime = dde.geometry.GeometryXTime(geom1|geom2, timedomain)

# Models and the respective domains
def pde_T1(x, y, _):
    dy_t = dde.grad.jacobian(y, x, i=0, j=1)
    dy_xx = dde.grad.hessian(y, x, i=0, j=0)
    return dy_t - a1 * dy_xx

def pde_T2(x, y, _):
    dy_t = dde.grad.jacobian(y, x, i=0, j=1)
    dy_xx = dde.grad.hessian(y, x, i=0, j=0)
    return dy_t - a2 * dy_xx

def on_domain1(x, on_domain):
    return geom1.inside(x)[0]

def on_domain2(x, on_domain):
    return geom2.inside(x)[0]

# Boundary and initial conditions
def on_boundary1(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0)

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

def boundary_initial(x, on_initial):
    return on_initial and np.isclose(x[1], 0)

# interface conditions
def on_interf(x, on_boundary):
    return on_boundary and np.isclose(x[0], L/2)

def flux_int(x,y,X): 
    # I need help here.  
    return (a1*geom1.boundary_normal(X) + a2*geom2.boundary_normal(X)).reshape(-1,1)

def Temp_int(x,y,X):
    # I need help here.
    # T1_int:  how to get from geom1 at x=L/2?
    # T2_int = how to get from geom2 at x=L/2?
    pass

# Setting the IC
def init_func(X):
    x = X[:, 0:1]
    y = X[:, 1:2]
    t = np.zeros((len(X),1))
    for count, x_ in enumerate(x):
        if x_ < L/2:
            t[count] = T0
        else:
            t[count] = T0 + 2*(Ts-T0) * (x_ - L/2)
    return t

ic = dde.IC(geomtime, init_func, boundary_initial)

# Seting the BCs
pde1 = dde.OperatorBC(geomtime1, pde_T1, on_boundary = on_domain1)
pde2 = dde.OperatorBC(geomtime2, pde_T2, on_boundary = on_domain2)
bc1 = dde.icbc.DirichletBC(geomtime1, lambda x: T0*np.ones((len(x),1)), on_boundary1)
bc2 = dde.icbc.DirichletBC(geomtime2, lambda x: TL*np.ones((len(x),1)), on_boundary2) # not used in loss

# Setting the BC at the interface with 500 points
X = np.hstack( (np.full((500), L/2).reshape(-1,1), timedomain.random_points(500))).reshape(-1, 2)
FluxInterf = dde.icbc.PointSetOperatorBC(X,
                                     np.zeros((X.shape[0],1)), # fluxes must add up to zero at x=L/2.
                                     lambda x, y, X : flux_int(x, y, X[:,0]))

# Setting the problem
loss = [pde1, pde2, bc1, ic, FluxInterf]
data = dde.data.TimePDE(
    geomtime,
    None,
    loss,
    num_domain=1000,
    num_boundary=500,
    num_initial=500,
    num_test=500)

loss_weights = [10, 10, 0.1, 0.1, 100]
net = dde.nn.FNN([2] + 4 * [50] + [1], "tanh", "Glorot normal")

# Enforcing BC at x=L
def  output_transform(x, y):    
    xx, t = x[:,0:1], x[:,1:2]
    return (L-xx)*y + Ts
net.apply_output_transform(output_transform)

model = dde.Model(data, net)

model.compile("adam", lr=1.0e-3, loss_weights = loss_weights)
losshistory, train_state = model.train(iterations=25000)

model.compile("L-BFGS")
losshistory, train_state = model.train()

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

I got the train loss and the results below:

[2.50e-04, 4.78e-04, 4.05e-05, 5.16e-04, 1.28e+00] 

Thank you for your time and consideration.

Best regards.

[lululxvi]

I found a way to enforce the heat flux continuity using geom.boundary_normal( ). But, the respective loss is not decreasing (in fact, it is constant). Is it correct to use geom.boundary_normal( )? Is there an alternative way?

You can write the loss explicitly, instead of using geom.boundary_normal()

I am struggling to come up with away to enforce T_1 - T_2 = 0. How could I get the values of T_1 and T_2 at x=L/2 during the model train?

You can use PointSetOperatorBC.

[w-venton]

@lululxvi I am also interested in this question, if there are two materials, can I use OperatorBC to make them meet the continuity of temperature and heat flow? like this: def func(x): return x[:,0:1]-0.5

bc1 = dde.OperatorBC(geom, func, lambda x: np.isclose(x[:, 0],0.5)) bc2 = dde.OperatorBC(geom, func, lambda x: np.isclose(x[:, 0],0.5),derivative_order=1)

then in pde(x,y) ,return k(dT_xx+dT_yy)-source+bc1.penaltybc1.func(x)_(T-bc2.value)+bc2.penaltybc2.func(x)_(kdT_x-k*bc2.value)

can you give me some suggestions? Or is there any other method that can solve the problem?

[lululxvi]

if there are two materials, can I use OperatorBC to make them meet the continuity of temperature and heat flow? like this: def func(x): return x[:,0:1]-0.5

What do you try to do exactly?

[juliosdutra]

I found a way to enforce the heat flux continuity using geom.boundary_normal( ). But, the respective loss is not decreasing (in fact, it is constant). Is it correct to use geom.boundary_normal( )? Is there an alternative way?

You can write the loss explicitly, instead of using geom.boundary_normal()

I am struggling to come up with away to enforce T_1 - T_2 = 0. How could I get the values of T_1 and T_2 at x=L/2 during the model train?

You can use PointSetOperatorBC.

Hi Lulu,

I am still working on the problem of finding a solution to a heat transfer problem in a material with different coefficients of thermal conductivity. I have been using tf.where to determine where T1 and T2 should be, and I have also used the PointSetOperatorBC operator to enforce the two interface conditions of thermal flow and temperature continuities, as you suggested.

At this point, I would like to ask for your feedback on whether you believe this approach (shown below) can be considered a solution to the problem, or if there are any further changes or suggestions you would make to improve the solution. Any help or guidance you can offer would be greatly appreciated.

Thank you.

Code

(just copied the essential parts)

# Parameters and geometry

a1 = 1.14 # conductivity coefficient 1 (left side)
a2 = 0.01 # conductivity coefficient 2 (right side)
L = 1.0  # total length
T0 = 0.0  # temperature specified at x=0
TL = 1.0  # temperature specified at x=L
tend = 1.0  # final time for the simulation

geom1 = dde.geometry.Interval(0, L/2) # first solid
geom2 = dde.geometry.Interval(L/2, L) # second solid

timedomain = dde.geometry.TimeDomain(0, tend)
geomtime = dde.geometry.GeometryXTime(geom1|geom2, timedomain)

# System pde

def pde(x, y):    
    T1 = tf.where(tf.less(x[:,0:1],L/2), y[:,0:1], y[:,1:2])
    T2 = tf.where(tf.greater(x[:,0:1],L/2), y[:,1:2], y[:,0:1])

    dT1_t = dde.grad.jacobian(T1, x, i=0, j=1)
    dT2_t = dde.grad.jacobian(T2, x, i=0, j=1)

    dT1_xx = dde.grad.hessian(T1, x, i=0, j=0)
    dT2_xx = dde.grad.hessian(T2, x, i=0, j=0)      
    return tf.where(x[:,0:1]<L/2, dT1_t - a1 * dT1_xx, dT2_t - a2 * dT2_xx)

# Flux continuity

def continuity_flux(x, y, X):
    T1 = tf.where(tf.less(x[:,0:1],L/2), y[:,0:1], y[:,1:2])
    T2 = tf.where(tf.greater(x[:,0:1],L/2), y[:,1:2], y[:,0:1])

    dT1_x = dde.grad.jacobian(T1, x, i=0, j=0)
    dT2_x = dde.grad.jacobian(T2, x, i=0, j=0)

    return a2*dT2_x-a1*dT1_x

# Temperature continuity

def continuity_temp(x, y, X):
    T1 = tf.where(tf.less(x[:,0:1],L/2), y[:,0:1], y[:,1:2])
    T2 = tf.where(tf.greater(x[:,0:1],L/2), y[:,1:2], y[:,0:1])
    return T2-T1

# Interface condition

x_interface = np.vstack((np.full((100), L/2), np.linspace(0, tend, num=100))).T
bc_continuity1 = dde.icbc.PointSetOperatorBC(x_interface, 0.0 , continuity_flux)
bc_continuity2 = dde.icbc.PointSetOperatorBC(x_interface, 0.0 , continuity_temp)

# Model training

BC = [bc_continuity1, bc_continuity2, ic]
loss_weights = [1, 1, 1, 100]

data = dde.data.TimePDE(
    geomtime, pde, BC,
    num_domain=3000,
    num_boundary=500,
    num_initial=500,
    num_test=5000,
    anchors=x_interface,)

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

def  output_transform(x, y):    
    X = x[:,0:1]
    return X*(L-X)*y + X*Ts + (L-X)*T0

net.apply_output_transform(output_transform)

model = dde.Model(data, net)

model.compile("adam", lr=0.001, loss_weights = loss_weights)
losshistory, train_state = model.train(iterations=20000, display_every=1000)

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

Results

PDE residue and the complet solution.

Main result

Residue at the interface x=L/2.

[w-venton]

@juliosdutra hi, I am also interested in the continuity of heat flux of two materials. Your result looks good, have you tested its accuracy? And I have a question. What do you think about using a large rectangle and then using torch.where for the different coefficients? Then is it feasible to use operatorBC at the interface. In this case, I don't quite know how to get dT/dx of the two materials on interface. In addition, I think that the data trained with a rectangle should already satisfy T continuous, so can it not require additional constraints? Do you have any idea?

[w-venton]

@juliosdutra hi, I want to ask what's the meaning of y[:,0:1] and y[:,1:2]？ What is his physical meaning？

[yuluoxiansen]

@juliosdutra hi, I want to ask what's the meaning of y[:,0:1] and y[:,1:2]？ What is his physical meaning？ 你好，我想问一下 y[:,0:1] 和 y[:,1:2] 是什么意思？他的物理意义是什么？

maybe T1 and T2