# 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.

## 1 Answer

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!

• I've also found method of lines very reliable for reaction diffusion problems, has the advantage that you only need to think about the spatial discretisation and let the ODE solver worry about time. – boyfarrell Mar 16 '15 at 19:29
• @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
• Ah ok, of course it makes much more sense, retrieving the information on the go. Thanks! @boyfarrell: Indeed it is a good method for this type of problem, I had completely forgot to look into it, thanks for the idea. – Londero Mar 17 '15 at 17:08