Unsteady diffusion + reaction 1/2/3D - Multiphase Reactor Modeling using PyMRM

Reaction: $A \rightarrow B\rightarrow C$

$\frac{\partial c_A}{\partial t} + \mathrm{div}(-D \, \mathrm{grad}(c_A)) = -k_1 \, c_A$

$\frac{\partial c_B}{\partial t} + \mathrm{div}(-D \, \mathrm{grad}(c_B)) = k_1 \, c_A - k_2 \, c_B$

$\frac{\partial c_C}{\partial t} + \mathrm{div}(-D \, \mathrm{grad}(c_C)) = k_2 \, c_B$

with boundary conditions: $c_A= c_w , c_B=0, c_C(0)=0$

The setup of this problem is more complex than is stricktly necesarry. The puprose is to demonstrate that models can be build in such a way that they can be used in 1, 2 or 3 dimensions. For 1D and 2D it was found that the direct sparse LU solver works well. For 3D it seems very slow and BiCGStab with an ILU preconditioner is used. Note that the ILU preconditioner takes some time in this case.

All boundaries use a Dirichlet condition, {'a': 0, 'b': 1, 'd': c_b}, applied with the outward normal convention (see the conventions skill’s pitfalls note P1), so the sign of the boundary flux contribution flips between opposite faces of the domain.

Imports¶

import math
import numpy as np
from scipy.sparse import eye_array, linalg as sla
import matplotlib.pyplot as plt
from IPython.display import clear_output, display
from pymrm import construct_grad, construct_div, NumJac

Reaction term and parameters¶

def reaction(c, k_1, k_2, axis=-1):
    shape_c = c.shape
    shape_t = [math.prod(shape_c[:axis]), shape_c[axis], -1]
    c = c.reshape(shape_t)
    r = np.zeros_like(c)
    r[:, 0, :] = -k_1 * c[:, 0, :]
    r[:, 1, :] = k_1 * c[:, 0, :] - k_2 * c[:, 1, :]
    r[:, 2, :] = k_2 * c[:, 1, :]
    return r.reshape(shape_c)

dim = 2
L = [0.2] * dim
L[0] = 0.1
D = 1e-4
k_1 = 2.0
k_2 = 1.0
t_end = 5
num_c = 3
c_b = np.array([1.0] + [0.0] * (num_c - 1)).reshape((1,) * (dim - 1) + (-1,)).reshape((1,) * dim + (-1,))

shape_c = [31] * dim + [num_c]
num_time_steps = 1000
num_inner_iter = 1
output_interval = 100  # num_time_steps
dt = t_end / num_time_steps
c = np.zeros(shape_c)
numjac = NumJac(shape_c)

Grid, operators and Jacobian assembly¶

Build the grid and the diffusion operators per spatial axis, accumulate them into a constant Jacobian jac and boundary term g_const, subtract the (linear) reaction Jacobian, and factorize jac once (LU for dim != 3, ILU-preconditioned for dim == 3).

# All boundaries Dirichlet bc, with c=0. By making sizes=1 broadcasting will be used.
bc = {'a': 0, 'b': 1, 'd': c_b}
x_f = [None] * dim
x_c = [None] * dim
flux_mat = [None] * dim
flux_bc = [None] * dim
jac = eye_array(math.prod(shape_c), format='csc') / dt
g_const = np.zeros((math.prod(shape_c), 1))
for i in range(dim):
    x_f[i] = np.linspace(0, L[i], shape_c[i] + 1)
    x_c[i] = 0.5 * (x_f[i][:-1] + x_f[i][1:])
    grad_mat, grad_bc = construct_grad(shape_c, x_f[i], x_c[i], bc=(bc, bc), axis=i)
    flux_mat[i] = -D * grad_mat
    flux_bc[i] = -D * grad_bc
    div_mat = construct_div(shape_c, x_f[i], nu=0, axis=i)
    jac += div_mat @ flux_mat[i]
    g_const += div_mat @ flux_bc[i]
f_react, jac_react = numjac(lambda c: reaction(c, k_1, k_2, axis=-1), c)
jac_react.eliminate_zeros()
jac -= jac_react
if dim == 3:
    jac_ilu = sla.spilu(jac)
    jac_pc = sla.LinearOperator(jac_ilu.shape, lambda x: jac_ilu.solve(x))
else:
    jac_lu = sla.splu(jac)

Plot setup¶

labels = [r'$A$', r'$B$', r'$C$']
if dim == 1:
    fig, ax = plt.subplots()
    lines = [ax.plot(x_c[0], c[:, j], label=labels[j])[0] for j in range(num_c)]
    plt.ylim(0, 1.1)
    ax.set_xlabel('position')
    ax.set_ylabel('concentration')
    plt.legend()
elif (dim == 2) | (dim == 3):
    fig, ax = plt.subplots(1, num_c, figsize=(8, 2))
    for i in range(num_c):
        ax[i].set_xticks([])
        ax[i].set_yticks([])
        ax[i].set_xticklabels([])
        ax[i].set_yticklabels([])
        ax[i].set_title(f'species {labels[i]}')
    contour = [None] * num_c
    if dim == 2:
        for i in range(num_c):
            contour[i] = ax[i].pcolormesh(x_f[0], x_f[1], c[:, :, i], shading='flat', cmap='viridis', vmin=0, vmax=1)
    else:
        indx_z = math.floor(shape_c[dim - 1] / 2)
        for i in range(num_c):
            contour[i] = ax[i].pcolormesh(x_f[0], x_f[1], c[:, :, indx_z, i], shading='flat', cmap='viridis', vmin=0, vmax=1)
    cbar_ax = fig.add_axes([0.92, ax[0].get_position().y0, 0.02, ax[0].get_position().y1 - ax[0].get_position().y0])  # [left, bottom, width, height]
    plt.colorbar(contour[0], cax=cbar_ax)
plt.show()

Time stepping¶

for i in range(num_time_steps):
    c_old = c.copy().reshape((-1, 1))
    for j in range(num_inner_iter):
        g = g_const + jac @ c.reshape((-1, 1)) - c_old / dt
        if dim == 3:
            dc, exit_code = sla.bicgstab(jac, g, M=jac_pc)
            c -= dc.reshape(c.shape)
        else:
            c -= jac_lu.solve(g).reshape(c.shape)
    if (i + 1) % output_interval == 0:
        clear_output(wait=True)
        if dim == 1:
            for j in range(num_c):
                lines[j].set_ydata(c[:, j])
        elif dim == 2:
            for i in range(num_c):
                contour[i].set_array(c[:, :, i].ravel())
        elif dim == 3:
            for i in range(num_c):
                contour[i].set_array(c[:, :, indx_z, i].ravel())
        display(fig)

Checking the result¶

The reaction term is linear in c, so the whole implicit-Euler step is linear too: jac and g_const were assembled once and reused for every time step, and they define a residual

g(c) = g_\mathrm{const} + \mathrm{jac}\, c - c_\mathrm{old}/dt

(1)

whose exact Jacobian is jac itself. The linear systems are solved iteratively, so evaluating this residual with pymrm.checks.residual_check shows whether the last linear solve converged to a scale-free backward error. It cannot catch a mis-assembled operator: the residual is built from the same matrices.

Note that c_old still holds the state from the last time step actually taken in the loop above (it is reassigned to c.copy() at the top of each iteration), so this checks the step from that state to the final c, not a step starting from c itself.

from pymrm.checks import residual_check


def _resid(c_flat):
    return (g_const + jac @ c_flat.reshape((-1, 1)) - c_old.reshape((-1, 1)) / dt).ravel(), jac


result = residual_check(_resid, c.ravel())
print(result)
if result['ok']:
    print('Residual check passed: the assembled operators are consistent with the last time step taken.')
else:
    print('Residual check failed: inspect the sign or assembly of the transport/reaction operators.')

{'ok': True, 'backward_error': 2.955238379259546e-16, 'max_abs_residual': 8.526512829121202e-14}
Residual check passed: the assembled operators are consistent with the last time step taken.