Inverse problem for the diffusion-reaction system

Inverse problem for the diffusion-reaction system#

The problem we are solving#

We will solve an inverse problem for diffusion-reaction systems for unknowns \(D\) and \(k_f\):

\[\frac{\partial C_A}{\partial t}=D\frac{\partial^2C_A}{\partial x^2}-k_fC_AC_B^2,\]
\[\frac{\partial C_B}{\partial t}=D\frac{\partial^2C_B}{\partial x^2}-2k_fC_AC_B^2\]

for \(x\in[0,1]\) and \(t\in[0,10]\) initial conditions

\[C_A(x,0)=C_B(x,0)=e^{-20x}\]

and boundary conditions

\[\begin{split}C_A(0,t)=C_B(0,t)=1\\C_A(1,t)=C_B(1,t)=0.\end{split}\]

The expected values of \(D\) and \(k_f\) are \(2\cdot10^{-3}\) and 0.1, respectively.

Below is the dimensional analysis for the given diffusion-reaction system with unknowns \( D \) and \( k_f \):

The equations are:
$\( \frac{\partial C_A}{\partial t} = D \frac{\partial^2 C_A}{\partial x^2} - k_f C_A C_B^2, \quad \frac{\partial C_B}{\partial t} = D \frac{\partial^2 C_B}{\partial x^2} - 2k_f C_A C_B^2. \)$

Dimensional Analysis#

Step 1: Assign Dimensions to Variables#

  1. Spatial Coordinate \(x\):

    • The dimension of \(x\) is length:

      \[ [x] = L. \]

  2. Time \(t\):

    • The dimension of time is:

      \[ [t] = T. \]

  3. Concentrations \(C_A\) and \(C_B\):

    • Concentrations typically have dimensions of moles per unit volume:

      \[ [C_A] = [C_B] = N / L^3, \]

      where \(N\) represents the amount of substance (moles).

  4. Diffusion Coefficient \(D\):

    • The diffusion term involves \(\partial^2 C_A / \partial x^2\), which has dimensions:

      \[ \left[\frac{\partial^2 C_A}{\partial x^2}\right] = \frac{[C_A]}{[x]^2} = \frac{N}{L^5}. \]

    • The time derivative \(\partial C_A / \partial t\) has dimensions:

      \[ \left[\frac{\partial C_A}{\partial t}\right] = \frac{[C_A]}{[t]} = \frac{N}{L^3 T}. \]

    • For dimensional consistency in the diffusion term \(D \frac{\partial^2 C_A}{\partial x^2}\), \(D\) must have dimensions:

      \[ [D] \cdot \frac{N}{L^5} = \frac{N}{L^3 T} \implies [D] = \frac{L^2}{T}. \]

  5. Reaction Rate Constant \(k_f\):

    • The reaction term \(k_f C_A C_B^2\) must have the same dimensions as \(\partial C_A / \partial t\), which is \(N / (L^3 T)\).

    • The product \(C_A C_B^2\) has dimensions:

      \[ [C_A C_B^2] = \left(\frac{N}{L^3}\right)\left(\frac{N}{L^3}\right)^2 = \frac{N^3}{L^9}. \]

    • To balance dimensions, \(k_f\) must satisfy:

      \[ [k_f] \cdot \frac{N^3}{L^9} = \frac{N}{L^3 T} \implies [k_f] = \frac{L^6}{N^2 T}. \]

Step 2: Initial and Boundary Conditions#

  • Initial Conditions: \(C_A(x,0) = C_B(x,0) = e^{-20x}\)

    • The exponential function is dimensionless, so the initial condition is consistent with \(C_A\) and \(C_B\) having dimensions \(N / L^3\).

  • Boundary Conditions: \(C_A(0,t) = C_B(0,t) = 1\) and \(C_A(1,t) = C_B(1,t) = 0\)

    • These conditions specify dimensionless values for \(C_A\) and \(C_B\) scaled appropriately.

Step 3: Summary of Dimensions#

Variable/Parameter

Physical Meaning

Dimensions

\(x\)

Spatial coordinate

\(L\)

\(t\)

Time

\(T\)

\(C_A, C_B\)

Concentrations (solute)

\(N / L^3\)

\(D\)

Diffusion coefficient

\(L^2 / T\)

\(k_f\)

Reaction rate constant

\(L^6 / (N^2 T)\)

Step 4: Given Values#

  • \(D = 2 \cdot 10^{-3}\): The value is consistent with \([D] = L^2 / T\).

  • \(k_f = 0.1\): The value is consistent with \([k_f] = L^6 / (N^2 T)\).

These values and dimensions are appropriate for the given system.

Implementation#

Below is the implementation of the inverse problem for the diffusion-reaction system using PINNx.

We first load the dataset and define the PDE and neural network model. We then define the boundary and initial conditions, create the geometry, and set up the problem. Finally, we generate training data and solve the inverse problem to estimate the diffusion coefficient \(D\) and reaction rate constant \(k_f\).

import brainstate
import brainunit as u
import numpy as np

import pinnx

Define the physical units for the problem:

unit_of_x = u.meter
unit_of_t = u.second
unit_of_c = u.mole / u.meter ** 3

Define the parameters and the PDE:


kf = brainstate.ParamState(0.05 * u.meter ** 6 / u.mole ** 2 / u.second)
D = brainstate.ParamState(1.0 * u.meter ** 2 / u.second)


def pde(x, y):
    jacobian = net.jacobian(x, x='t')
    hessian = net.hessian(x)
    ca, cb = y['ca'], y['cb']
    dca_t = jacobian['ca']['t']
    dcb_t = jacobian['cb']['t']
    dca_xx = hessian['ca']['x']['x']
    dcb_xx = hessian['cb']['x']['x']
    eq_a = dca_t - 1e-3 * D.value * dca_xx + kf.value * ca * cb ** 2
    eq_b = dcb_t - 1e-3 * D.value * dcb_xx + 2 * kf.value * ca * cb ** 2
    return [eq_a, eq_b]

Define the neural network model:

net = pinnx.nn.Model(
    pinnx.nn.DictToArray(x=unit_of_x, t=unit_of_t),
    pinnx.nn.FNN([2] + [20] * 3 + [2], "tanh"),
    pinnx.nn.ArrayToDict(ca=unit_of_c, cb=unit_of_c),
)

Define the geometry:

geom = pinnx.geometry.Interval(0, 1)
timedomain = pinnx.geometry.TimeDomain(0, 10)
geomtime = pinnx.geometry.GeometryXTime(geom, timedomain)
geomtime = geomtime.to_dict_point(x=unit_of_x, t=unit_of_t)

Define the boundary condition function:

def fun_bc(x):
    c = (1 - x['x'] / unit_of_x) * unit_of_c
    return {'ca': c, 'cb': c}


bc = pinnx.icbc.DirichletBC(fun_bc)

Define the initial condition function:

def fun_init(x):
    return {
        'ca': u.math.exp(-20 * x['x'] / unit_of_x) * unit_of_c,
        'cb': u.math.exp(-20 * x['x'] / unit_of_x) * unit_of_c,
    }


ic = pinnx.icbc.IC(fun_init)

Define a function to generate training data:


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 = {'x': X.flatten() * unit_of_x, 't': T.flatten() * unit_of_t}
    y = {'ca': ca.flatten() * unit_of_c, 'cb': cb.flatten() * unit_of_c}
    return x, y


observe_x, observe_y = gen_traindata()
observe_bc = pinnx.icbc.PointSetBC(observe_x, observe_y)

Define the problem and solve the inverse problem:

data = pinnx.problem.TimePDE(
    geomtime,
    pde,
    [bc, ic, observe_bc],
    net,
    num_domain=2000,
    num_boundary=100,
    num_initial=100,
    anchors=observe_x,
    num_test=50000,
)

Warning: 50000 points required, but 50055 points sampled.

variable = pinnx.callbacks.VariableValue([kf, D], period=1000, filename="./variables.dat")
model = pinnx.Trainer(data, external_trainable_variables=[kf, D])
model.compile(braintools.optim.Adam(0.001)).train(iterations=80000, callbacks=[variable])
model.saveplot(issave=True, isplot=True)

Compiling trainer...
'compile' took 0.055027 s

Training trainer...

Step      Train loss                                              Test loss                                               Test metric 
0         [0.09352963 * 10.0^0 * (mmolar / second) ** 2,          [0.08070567 * 10.0^0 * (mmolar / second) ** 2,          []          
           0.30233935 * 10.0^0 * (mmolar / second) ** 2,           0.28232208 * 10.0^0 * (mmolar / second) ** 2,                      
           {'ibc0': {'ca': 2.5078168 * mmolar,                     {'ibc0': {'ca': 2.5078168 * mmolar,                                
                     'cb': 2.447632 * mmolar}},                              'cb': 2.447632 * mmolar}},                               
           {'ibc1': {'ca': 1.8000234 * mmolar,                     {'ibc1': {'ca': 1.8000234 * mmolar,                                
                     'cb': 0.0453864 * mmolar}},                             'cb': 0.0453864 * mmolar}},                              
           {'ibc2': {'ca': 3.1355276 * mmolar,                     {'ibc2': {'ca': 3.1355276 * mmolar,                                
                     'cb': 1.4068981 * mmolar}}]                             'cb': 1.4068981 * mmolar}}]                              
1000      [0.00314551 * 10.0^0 * (mmolar / second) ** 2,          [0.0021539 * 10.0^0 * (mmolar / second) ** 2,           []          
           0.00141168 * 10.0^0 * (mmolar / second) ** 2,           0.00109254 * 10.0^0 * (mmolar / second) ** 2,                      
           {'ibc0': {'ca': 0.0053007 * mmolar,                     {'ibc0': {'ca': 0.0053007 * mmolar,                                
                     'cb': 0.0027947 * mmolar}},                             'cb': 0.0027947 * mmolar}},                              
           {'ibc1': {'ca': 0.00446298 * mmolar,                    {'ibc1': {'ca': 0.00446298 * mmolar,                               
                     'cb': 0.00301004 * mmolar}},                            'cb': 0.00301004 * mmolar}},                             
           {'ibc2': {'ca': 0.00534556 * mmolar,                    {'ibc2': {'ca': 0.00534556 * mmolar,                               
                     'cb': 0.00337772 * mmolar}}]                            'cb': 0.00337772 * mmolar}}]                             
2000      [0.00084838 * 10.0^0 * (mmolar / second) ** 2,          [0.00054396 * 10.0^0 * (mmolar / second) ** 2,          []          
           0.00066913 * 10.0^0 * (mmolar / second) ** 2,           0.0005198 * 10.0^0 * (mmolar / second) ** 2,                       
           {'ibc0': {'ca': 0.00092106 * mmolar,                    {'ibc0': {'ca': 0.00092106 * mmolar,                               
                     'cb': 0.00056548 * mmolar}},                            'cb': 0.00056548 * mmolar}},                             
           {'ibc1': {'ca': 0.00122206 * mmolar,                    {'ibc1': {'ca': 0.00122206 * mmolar,                               
                     'cb': 0.00080589 * mmolar}},                            'cb': 0.00080589 * mmolar}},                             
           {'ibc2': {'ca': 0.00065078 * mmolar,                    {'ibc2': {'ca': 0.00065078 * mmolar,                               
                     'cb': 0.00079521 * mmolar}}]                            'cb': 0.00079521 * mmolar}}]                             
3000      [0.00039193 * 10.0^0 * (mmolar / second) ** 2,          [0.00028126 * 10.0^0 * (mmolar / second) ** 2,          []          
           0.00028135 * 10.0^0 * (mmolar / second) ** 2,           0.00023285 * 10.0^0 * (mmolar / second) ** 2,                      
           {'ibc0': {'ca': 0.00026876 * mmolar,                    {'ibc0': {'ca': 0.00026876 * mmolar,                               
                     'cb': 0.00018014 * mmolar}},                            'cb': 0.00018014 * mmolar}},                             
           {'ibc1': {'ca': 0.00043916 * mmolar,                    {'ibc1': {'ca': 0.00043916 * mmolar,                               
                     'cb': 0.0003322 * mmolar}},                             'cb': 0.0003322 * mmolar}},                              
           {'ibc2': {'ca': 0.00018346 * mmolar,                    {'ibc2': {'ca': 0.00018346 * mmolar,                               
                     'cb': 0.00025419 * mmolar}}]                            'cb': 0.00025419 * mmolar}}]                             
4000      [0.00021093 * 10.0^0 * (mmolar / second) ** 2,          [0.00015969 * 10.0^0 * (mmolar / second) ** 2,          []          
           0.00016543 * 10.0^0 * (mmolar / second) ** 2,           0.00013996 * 10.0^0 * (mmolar / second) ** 2,                      
           {'ibc0': {'ca': 0.00010608 * mmolar,                    {'ibc0': {'ca': 0.00010608 * mmolar,                               
                     'cb': 9.404449e-05 * mmolar}},                          'cb': 9.404449e-05 * mmolar}},                           
           {'ibc1': {'ca': 0.00023372 * mmolar,                    {'ibc1': {'ca': 0.00023372 * mmolar,                               
                     'cb': 0.0001386 * mmolar}},                             'cb': 0.0001386 * mmolar}},                              
           {'ibc2': {'ca': 9.680791e-05 * mmolar,                  {'ibc2': {'ca': 9.680791e-05 * mmolar,                             
                     'cb': 0.00014055 * mmolar}}]                            'cb': 0.00014055 * mmolar}}]                             
5000      [0.00014335 * 10.0^0 * (mmolar / second) ** 2,          [0.00011242 * 10.0^0 * (mmolar / second) ** 2,          []          
           0.00011522 * 10.0^0 * (mmolar / second) ** 2,           0.00010192 * 10.0^0 * (mmolar / second) ** 2,                      
           {'ibc0': {'ca': 5.779891e-05 * mmolar,                  {'ibc0': {'ca': 5.779891e-05 * mmolar,                             
                     'cb': 5.5820285e-05 * mmolar}},                         'cb': 5.5820285e-05 * mmolar}},                          
           {'ibc1': {'ca': 0.00013723 * mmolar,                    {'ibc1': {'ca': 0.00013723 * mmolar,                               
                     'cb': 6.1033348e-05 * mmolar}},                         'cb': 6.1033348e-05 * mmolar}},                          
           {'ibc2': {'ca': 7.234841e-05 * mmolar,                  {'ibc2': {'ca': 7.234841e-05 * mmolar,                             
                     'cb': 0.00010172 * mmolar}}]                            'cb': 0.00010172 * mmolar}}]                             
6000      [0.00011602 * 10.0^0 * (mmolar / second) ** 2,          [9.5097355e-05 * 10.0^0 * (mmolar / second) ** 2,       []          
           7.968601e-05 * 10.0^0 * (mmolar / second) ** 2,         7.246479e-05 * 10.0^0 * (mmolar / second) ** 2,                    
           {'ibc0': {'ca': 3.963989e-05 * mmolar,                  {'ibc0': {'ca': 3.963989e-05 * mmolar,                             
                     'cb': 3.690526e-05 * mmolar}},                          'cb': 3.690526e-05 * mmolar}},                           
           {'ibc1': {'ca': 8.915077e-05 * mmolar,                  {'ibc1': {'ca': 8.915077e-05 * mmolar,                             
                     'cb': 3.357489e-05 * mmolar}},                          'cb': 3.357489e-05 * mmolar}},                           
           {'ibc2': {'ca': 5.7701498e-05 * mmolar,                 {'ibc2': {'ca': 5.7701498e-05 * mmolar,                            
                     'cb': 7.855295e-05 * mmolar}}]                          'cb': 7.855295e-05 * mmolar}}]                           
7000      [9.821482e-05 * 10.0^0 * (mmolar / second) ** 2,        [8.3726925e-05 * 10.0^0 * (mmolar / second) ** 2,       []          
           5.896915e-05 * 10.0^0 * (mmolar / second) ** 2,         5.3390148e-05 * 10.0^0 * (mmolar / second) ** 2,                   
           {'ibc0': {'ca': 3.0593033e-05 * mmolar,                 {'ibc0': {'ca': 3.0593033e-05 * mmolar,                            
                     'cb': 2.6512069e-05 * mmolar}},                         'cb': 2.6512069e-05 * mmolar}},                          
           {'ibc1': {'ca': 6.14094e-05 * mmolar,                   {'ibc1': {'ca': 6.14094e-05 * mmolar,                              
                     'cb': 2.1582948e-05 * mmolar}},                         'cb': 2.1582948e-05 * mmolar}},                          
           {'ibc2': {'ca': 4.6279834e-05 * mmolar,                 {'ibc2': {'ca': 4.6279834e-05 * mmolar,                            
                     'cb': 6.889654e-05 * mmolar}}]                          'cb': 6.889654e-05 * mmolar}}]                           
8000      [8.164347e-05 * 10.0^0 * (mmolar / second) ** 2,        [7.186554e-05 * 10.0^0 * (mmolar / second) ** 2,        []          
           4.8710255e-05 * 10.0^0 * (mmolar / second) ** 2,        4.3789907e-05 * 10.0^0 * (mmolar / second) ** 2,                   
           {'ibc0': {'ca': 2.6550331e-05 * mmolar,                 {'ibc0': {'ca': 2.6550331e-05 * mmolar,                            
                     'cb': 3.1876243e-05 * mmolar}},                         'cb': 3.1876243e-05 * mmolar}},                          
           {'ibc1': {'ca': 4.271252e-05 * mmolar,                  {'ibc1': {'ca': 4.271252e-05 * mmolar,                             
                     'cb': 1.5517275e-05 * mmolar}},                         'cb': 1.5517275e-05 * mmolar}},                          
           {'ibc2': {'ca': 4.1019546e-05 * mmolar,                 {'ibc2': {'ca': 4.1019546e-05 * mmolar,                            
                     'cb': 6.213671e-05 * mmolar}}]                          'cb': 6.213671e-05 * mmolar}}]                           
9000      [6.772583e-05 * 10.0^0 * (mmolar / second) ** 2,        [6.0895793e-05 * 10.0^0 * (mmolar / second) ** 2,       []          
           4.3536915e-05 * 10.0^0 * (mmolar / second) ** 2,        3.854337e-05 * 10.0^0 * (mmolar / second) ** 2,                    
           {'ibc0': {'ca': 2.3455248e-05 * mmolar,                 {'ibc0': {'ca': 2.3455248e-05 * mmolar,                            
                     'cb': 1.5344798e-05 * mmolar}},                         'cb': 1.5344798e-05 * mmolar}},                          
           {'ibc1': {'ca': 3.100229e-05 * mmolar,                  {'ibc1': {'ca': 3.100229e-05 * mmolar,                             
                     'cb': 1.36119e-05 * mmolar}},                           'cb': 1.36119e-05 * mmolar}},                            
           {'ibc2': {'ca': 2.8181494e-05 * mmolar,                 {'ibc2': {'ca': 2.8181494e-05 * mmolar,                            
                     'cb': 6.931232e-05 * mmolar}}]                          'cb': 6.931232e-05 * mmolar}}]                           
10000     [5.5303026e-05 * 10.0^0 * (mmolar / second) ** 2,       [5.092607e-05 * 10.0^0 * (mmolar / second) ** 2,        []          
           3.807514e-05 * 10.0^0 * (mmolar / second) ** 2,         3.359254e-05 * 10.0^0 * (mmolar / second) ** 2,                    
           {'ibc0': {'ca': 1.35118025e-05 * mmolar,                {'ibc0': {'ca': 1.35118025e-05 * mmolar,                           
                     'cb': 1.4254373e-05 * mmolar}},                         'cb': 1.4254373e-05 * mmolar}},                          
           {'ibc1': {'ca': 2.082289e-05 * mmolar,                  {'ibc1': {'ca': 2.082289e-05 * mmolar,                             
                     'cb': 1.0991285e-05 * mmolar}},                         'cb': 1.0991285e-05 * mmolar}},                          
           {'ibc2': {'ca': 2.1368152e-05 * mmolar,                 {'ibc2': {'ca': 2.1368152e-05 * mmolar,                            
                     'cb': 5.3593893e-05 * mmolar}}]                         'cb': 5.3593893e-05 * mmolar}}]                          

Best trainer at step 10000:
  train loss: 2.28e-04
  test loss: 2.19e-04
  test metric: []

'train' took 1500.885740 s

Saving loss history to D:\codes\projects\pinnx\docs\unit-examples-inverse\loss.dat ...
Saving checkpoint into D:\codes\projects\pinnx\docs\unit-examples-inverse\loss.dat
Saving training data to D:\codes\projects\pinnx\docs\unit-examples-inverse\train.dat ...
Saving checkpoint into D:\codes\projects\pinnx\docs\unit-examples-inverse\train.dat
Saving test data to D:\codes\projects\pinnx\docs\unit-examples-inverse\test.dat ...
Saving checkpoint into D:\codes\projects\pinnx\docs\unit-examples-inverse\test.dat
../_images/ac3ed15a72f840e973bfa8de83bbdffb4883dbf2be7225991b69390e5b33c695.png ../_images/8b89a1846d87cd02670fd8ab7a632c8c5bd509d361075ab24520b5e7f99165a3.png ../_images/908fa1a6b7ebf2d1bdf8125c486ae6883fee044d8590901af87a0d504bc2c9bd.png

The inverse problem is solved, and the estimated values of \(D\) and \(k_f\) are obtained. The results are as follows:

kf

ParamState(
  value=0.01828403 * 10.0^0 * metre ** 6 * second ** -1 * mole ** -2
)

D

ParamState(
  value=1.4558828 * meter2 / second
)