Inverse problem for the diffusion-reaction system
Problem setup
We will solve an inverse problem for diffusion-reaction systems for unknowns \(D\) and \(k_f\):
for \(x \in [0, 1]\) and \(t \in [0, 10]\) initial conditions
and boundary conditions
The training dataset is here, and the expected values of \(D\) and \(k_f\) are \(2 \cdot 10^{-3}\) and \(0.1\), respectively.
Implementation
We first import DeepXDE and numpy (np):
import deepxde as dde
import numpy as np
Now, we define the unknown variables \(D\) and \(k_f\) with initial guesses of 1 and 0.05, respectively:
kf = dde.Variable(0.05)
D = dde.Variable(1.0)
We define the computational geometries by using the built-in Interval and TimeDomain classes and combining them with GeometryXTime:
geom = dde.geometry.Interval(0, 1)
timedomain = dde.geometry.TimeDomain(0, 10)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)
Now, we create the reaction-inverse PDE:
def pde(x, y):
ca, cb = y[:, 0:1], y[:, 1:2]
dca_t = dde.grad.jacobian(y, x, i=0, j=1)
dca_xx = dde.grad.hessian(y, x, component=0, i=0, j=0)
dcb_t = dde.grad.jacobian(y, x, i=1, j=1)
dcb_xx = dde.grad.hessian(y, x, component=1, i=0, j=0)
eq_a = dca_t - 1e-3 * D * dca_xx + kf * ca * cb ** 2
eq_b = dcb_t - 1e-3 * D * dcb_xx + 2 * kf * ca * cb ** 2
return [eq_a, eq_b]
Here, the first parameter is the \(t\)-coordinate, and it is represented by x. The second parameter has \(C_A\) and \(C_B\), which are represented by y. Then, we use dde.grad.jacobian and dde.grad.hessian to represent the desired first and second order partial derivatives.
Next, we consider the Dirichlet boundary conditions:
def fun_bc(x):
return 1 - x[:, 0:1]
Now, we define the initial conditions:
def fun_init(x):
return np.exp(-20 * x[:, 0:1])
Now that these are defined, we apply them:
bc_a = dde.icbc.DirichletBC(
geomtime, fun_bc, lambda _, on_boundary: on_boundary, component=0
)
bc_b = dde.icbc.DirichletBC(
geomtime, fun_bc, lambda _, on_boundary: on_boundary, component=1
)
ic1 = dde.icbc.IC(geomtime, fun_init, lambda _, on_initial: on_initial, component=0)
ic2 = dde.icbc.IC(geomtime, fun_init, lambda _, on_initial: on_initial, component=1)
Now, we generate the training data by getting it from here:
def gen_traindata():
data = np.load("dataset/reaction.npz")
t, x, ca, cb = data["t"], data["x"], data["Ca"], data["Cb"]
X, T = np.meshgrid(x, t)
X = np.reshape(X, (-1, 1))
T = np.reshape(T, (-1, 1))
Ca = np.reshape(ca, (-1, 1))
Cb = np.reshape(cb, (-1, 1))
return np.hstack((X, T)), Ca, Cb
After generating the data, we organize it:
observe_x, Ca, Cb = gen_traindata()
observe_y1 = dde.icbc.PointSetBC(observe_x, Ca, component=0)
observe_y2 = dde.icbc.PointSetBC(observe_x, Cb, component=1)
Now, we can define the TimePDE problem as follows:
data = dde.data.TimePDE(
geomtime,
pde,
[bc_a, bc_b, ic1, ic2, observe_y1, observe_y2],
num_domain=2000,
num_boundary=100,
num_initial=100,
anchors=observe_x,
num_test=50000,
)
We have 2000 training residual points in the domain, 100 points on the boundary, 100 points for the initial conditions, and 50000 to test the PDE residual. anchors specifies the training points as well.
Now, we create the network:
net = dde.nn.FNN([2] + [20] * 3 + [2], "tanh", "Glorot uniform")
This network has two inputs, one for the \(t\)-coordinate and one for the \(x\)-coordinate, and three hidden layers with 20 neurons each. The output layer has two outputs, one for \(C_A\) and one for \(C_B\). We also choose tanh to be the activation function, and the initializer is Glorot uniform.
Now, we create the Model and specify the optimizer, learning rate, and external_trainable_variables. We also output the values of \(D\) and \(k_f\) every 1000 iterations.
model = dde.Model(data, net)
model.compile("adam", lr=0.001, external_trainable_variables=[kf, D])
variable = dde.callbacks.VariableValue([kf, D], period=1000, filename="variables.dat")
Lastly, we train this network for 80000 iterations:
losshistory, train_state = model.train(iterations=80000, callbacks=[variable])
dde.saveplot(losshistory, train_state, issave=True, isplot=True)
Complete Code
"""Backend supported: tensorflow.compat.v1, tensorflow, pytorch, paddle"""
import deepxde as dde
import numpy as np
def gen_traindata():
data = np.load("../dataset/reaction.npz")
t, x, ca, cb = data["t"], data["x"], data["Ca"], data["Cb"]
X, T = np.meshgrid(x, t)
X = np.reshape(X, (-1, 1))
T = np.reshape(T, (-1, 1))
Ca = np.reshape(ca, (-1, 1))
Cb = np.reshape(cb, (-1, 1))
return np.hstack((X, T)), Ca, Cb
kf = dde.Variable(0.05)
D = dde.Variable(1.0)
def pde(x, y):
ca, cb = y[:, 0:1], y[:, 1:2]
dca_t = dde.grad.jacobian(y, x, i=0, j=1)
dca_xx = dde.grad.hessian(y, x, component=0, i=0, j=0)
dcb_t = dde.grad.jacobian(y, x, i=1, j=1)
dcb_xx = dde.grad.hessian(y, x, component=1, i=0, j=0)
eq_a = dca_t - 1e-3 * D * dca_xx + kf * ca * cb ** 2
eq_b = dcb_t - 1e-3 * D * dcb_xx + 2 * kf * ca * cb ** 2
return [eq_a, eq_b]
def fun_bc(x):
return 1 - x[:, 0:1]
def fun_init(x):
return np.exp(-20 * x[:, 0:1])
geom = dde.geometry.Interval(0, 1)
timedomain = dde.geometry.TimeDomain(0, 10)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)
bc_a = dde.icbc.DirichletBC(
geomtime, fun_bc, lambda _, on_boundary: on_boundary, component=0
)
bc_b = dde.icbc.DirichletBC(
geomtime, fun_bc, lambda _, on_boundary: on_boundary, component=1
)
ic1 = dde.icbc.IC(geomtime, fun_init, lambda _, on_initial: on_initial, component=0)
ic2 = dde.icbc.IC(geomtime, fun_init, lambda _, on_initial: on_initial, component=1)
observe_x, Ca, Cb = gen_traindata()
observe_y1 = dde.icbc.PointSetBC(observe_x, Ca, component=0)
observe_y2 = dde.icbc.PointSetBC(observe_x, Cb, component=1)
data = dde.data.TimePDE(
geomtime,
pde,
[bc_a, bc_b, ic1, ic2, observe_y1, observe_y2],
num_domain=2000,
num_boundary=100,
num_initial=100,
anchors=observe_x,
num_test=50000,
)
net = dde.nn.FNN([2] + [20] * 3 + [2], "tanh", "Glorot uniform")
model = dde.Model(data, net)
model.compile("adam", lr=0.001, external_trainable_variables=[kf, D])
variable = dde.callbacks.VariableValue([kf, D], period=1000, filename="variables.dat")
losshistory, train_state = model.train(iterations=80000, callbacks=[variable])
dde.saveplot(losshistory, train_state, issave=True, isplot=True)