How to use a discrete function in NeumannBC

[rodrigogdourado]

Hi, Dr. Lu Lu and community

I would like to implement the heat flux boundary condition shown below. However q(t) is a vector [q_1, q_2, q_3, ... , q_N] of discrete values associated to the vector of discrete times [t_1, t_2, t_3, ... , t_N]. I don't know how to do it using deepxde. This is similar to defining a discrete function inicial condition using PointSetBC, but in PointSetBC we calculate the difference T-T_0 (dirichlet BC) for example, not dT/dx.

Thank you!

[forxltk]

See #208 (comment)

[rodrigogdourado]

Thank you very much @forxltk ! This is exactly what I was trying to implement.

[rodrigogdourado]

Results of: 1D heat conduction with a heat flux boundary condition as a discrete function of time.

import deepxde as dde
import tensorflow as tf
import numpy as np
import numbers

# parameters

L = 1;
time =1;
q = 1;
T0 = 0;
k = 1.2;
rho = 1;
cp = 1;
alpha = k/(rho*cp)

#anchors
tt = np.arange(0,time+0.001,time/50).reshape(-1,1)
xx = np.zeros([len(tt),1])
ancora_x = np.hstack([xx,tt])

#heat flux function
fluxo = np.zeros([len(tt),1])
for i in range(len(tt)):
    if i>5 and i<=35:
        fluxo[i]=q

# PointSetBC for 1D neumann bc 
class PointSetBCDot(object):

    def __init__(self, points, values, component, j):
        self.points = np.array(points)
        if not isinstance(values, numbers.Number) and values.shape[1] != 1:
            raise RuntimeError(
                "PointSetBC should output 1D values. Use argument 'component' for different components."
            )
        self.j = j
        self.values = values
        self.component = component

    def first_normal_derivative(self, X, inputs, outputs, beg, end):
        dydt = dde.grad.jacobian(outputs, inputs, i=self.component, j=self.j)[beg:end]
        return tf.reduce_sum(dydt, axis=1, keepdims=True)

    def collocation_points(self, X):
        return self.points

    def error(self, X, inputs, outputs, beg, end):
        #return self.first_normal_derivative(X, inputs, outputs, beg, end) - tf.constant(self.values)
        return self.first_normal_derivative(X, inputs, outputs, beg, end) + self.values/k

# define a geometria do espaco e do tempo

geom = dde.geometry.Interval(0, L)
timedomain = dde.geometry.TimeDomain(0, time)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

# define o PDE e funcao de loss

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

# BCs
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)

bc_l = PointSetBCDot(ancora_x,fluxo,component=0,j=0)
bc_r = dde.icbc.NeumannBC(geomtime, lambda x:0, boundary_r)
ic = dde.icbc.IC(geomtime, lambda x:T0, lambda _, on_initial:on_initial)

data = dde.data.TimePDE(
    geomtime,
    pde,
    [bc_l,bc_r,ic],
    num_domain=500,
    num_boundary=100,
    num_initial=80,
    anchors = ancora_x,
    num_test=None,
    )

layer_size = [2]+[128]*2+[1]
net = dde.nn.FNN(layer_size,"tanh","Glorot normal")
#net = dde.nn.FNN(layer_size,"tanh","Glorot normal")
model = dde.Model(data,net)

model.compile("adam",lr=0.1e-3,loss_weights=[1,2,1,3])
losshistory, train_state = model.train(iterations=30000)

model.compile("L-BFGS",loss_weights=[1,2,1,3])
losshistory, train_state = model.train()

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

tt2 = np.arange(0,time+0.01,0.05)
xx = np.ones(len(tt2))*L

A = np.stack([xx,tt2])
A=A.T
print(A)

#%% roda o predict do modelo

resultado = model.predict(A)
print(resultado)

import matplotlib.pyplot as plt

#%% reference COMSOL Multiphysics

test = np.array([[0.0000,	0.0000],
[0.10000,	0.0000],
[0.20000,	0.012359],
[0.30000,	0.077224],
[0.40000,	0.16663],
[0.50000,	0.26301],
[0.60000,	0.36185],
[0.70000,	0.46133],
[0.80000,	0.54365],
[0.90000,	0.57987],
[1.0000,	0.59044]])
plt.plot(tt2,resultado)
plt.plot(test[:,0],test[:,1],'o',color='red')

The result is ok, but can be improved