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

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.

• 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.
• @Dr_Sam: Thanks for your answer. Both variables are indeed solved in the same grid (although the reaction takes place on only part of the grid). However I think I did not understand very well the part – "you do not need to interpolate the whole solution $u_A$ in the quadrature points, just its restriction to cell in which the quadrature points are." If I want to compute $\int u_A \phi_i dx$ using Gauss quadrature for example (2D), should not I compute it like this: $\sum_{i,j=0}^{n}u_{A}(q_{i})\cdot\phi_{i}(q_{i})\cdot w_{i}\cdot w_{j}.$ ? (Then, I need $u_A$ a quad. points $q_i$) – Londero Mar 17 '15 at 12:51
• The point is that you have $u_A(q_i) = \sum u_A^k \phi_k(q_i)$ because $u_A$ lives in the finite element space. While you are assembly the vector, you loop over the cells. Suppose that you are watching at the cell $K$. The quadrature points are located in that cell and so, you need to know $u_A$ just in that cell, but that's easy: it's $u_A(q_i) = \sum_{\{k \in I \}} u_A^k \phi_k(q_i)$ with $I$ the indices of the basis function in that cell. – Dr_Sam Mar 17 '15 at 13:13