Reaction-Diffusion problem A->B, solving for B

Question

I need to solve a Reaction-Diffusion using Finite Elements, piecewise linear elements. In this problem, a reaction $A \rightarrow B$, with rate law $ r_A = - k_A \cdot u_A $, takes part, where $u_i$ denotes concentration. Initially, $u_A = u_B = 0$. The time dependent formulation for the conservation of $A$ and $B$ is:

$ \begin{matrix} \frac{\partial \ u_A}{\partial \ t} -\Delta u_A - k \cdot u_A = f_A \end{matrix} $

$ \begin{matrix} \frac{\partial \ u_B}{\partial \ t} -\Delta u_B + k \cdot u_A = 0 \end{matrix}$

My question is: How is the best way to solve for $u_B$? Solving only $u_A$ seems trivial, using Crank-Nicholson as time discretization and finding a weak formulation that looks like:

$U^n[M+\delta_t\theta (A- kM)] = U^{n-1}[M+\delta_t(1-\theta) (-A+kM)] + \delta_t (F+ N)$

Where $U^n$ denotes the solution at time step $n$; $M, A, F, N$ are the mass matrix, stiffness matrix, load vector and B.C. vector (Neumann); $\delta_t$ the time step and $\theta$ the time discretization parameter ($1/2$ for C-N).

The first approach that I thought of is to solve for $u_A$ at each time step and then, using the value of $u_A$ at each time step, solve for $u_B$ using similar weak formulation:

$U^n[M+\delta_t\theta (A + kC)] = U^{n-1}[M+\delta_t(1-\theta) (-A-kC)]$

Where $kC$ takes the place of $kM$ in the first equation, and $C = \int u_A \ \phi_i dx$.

This would imply evaluating the matrix $C$ at each time step, using the $u_A$ calculated on the previous time step, interpolated on the quadrature points, which will certainly delay the process.

Is this the standard procedure or is there a different approach, more efficient? Thank you in advance.

Answer 1

I would do it exactly how you describe, with the two following remarks.

• I would solve the two evolutions simultaneously. In the way you describe, you will need to store and retrieve all the time steps for the solution $u_A$. Instead, you could solve sequentially for $u_A$ and $u_B$ at each time step, since $u_B$ just requires two solutions for $u_A$.
• If you assume that the solutions for $u_A$ and $u_B$ are solved on the same grid, then the computation of the vector $C$ is not that much expensive: because you loop on the cells to assemble the vector, you do not need to interpolate the whole solution $u_A$ in the quadrature points, but just its restriction to cell in which the quadrature points are.

Hope it helps!