# Solving systems of advection-diffusion-reaction equations with finite element methods

I have been doing a lot of self-study on numerically solving PDEs so that I can solve system of linear and nonlinear Advection-Diffusion-Reaction (ADR) systems on complex meshes. I have been watching some wonderful videos by Qiqi Wang, as well as reading the books by Langtangen and by Larson and Bengzon. I have also used Elmer FEM the package for solving some of these problems, but that system relatively self-contained with specific equation and solver setups.

There is of course a lot of machinery that goes into finite element methods, etc. However, I am having trouble understanding one simple question. How does one solve systems of ADRs using Finite Elements? Part of the problem is pedagogical. In many of FEM textbooks, the authors show how to solve the 1d and 2d Poisson equation, and heat equation, but then the focus shifts to solving specific commonly used systems of PDEs, like the Navier-Stokes equations, or the Euler equations, or Maxwell equations, etc. But these specialized systems have very specialized optimizations that have been developed over decades, so the implementation of those systems in FEM is not going to generalize to the wider class of ADR systems. It is kinda like trying to learn the simple principles for ADR systems by looking at really really complicated and optimized implementations of some code.

So let me start with a prototype system of ADR equations, and then identify the numerical questions about solving this system using FEM.I took this example from the article "Pattern Formation and Transition to Chaos in a Chemotaxis Model of Acute Inflammation" by Giunta, et al., 2021. Ideally I would solve this problem on a 3 dimensional mesh, but solving on a 1 or 2 dimensional mesh is instructive enough.

$$\frac{\partial m}{\partial t} = \nabla_x \cdot (D_m \nabla_x m) - \nabla_x \cdot (\psi \frac{m}{(1 + \alpha c)^2}\nabla_x c) + rmc(1 - \frac{m}{\bar{m}})$$

$$\frac{\partial c}{\partial t} = \nabla_x \cdot (D_c\nabla_x c) + \nu_c\frac{m}{1 + \beta a^\rho} - \mu_c c$$

$$\frac{\partial a}{\partial t} = \nabla_x \cdot (D_a\nabla_x a) + \nu_a \frac{m}{1 + \beta a^\rho} - \mu_a a$$

1. So I imagine I need to solve for the weak form of each of the 3 equations. So would I just multiply each equation by the test function, and then use integration by parts to compute the weak form of each equation? Or is there something where I have to convert this system to a matrix equation first, and then compute the weak form of the matrix, etc?

2. Once I compute the weak form, how would I proceed with implementing solve. So I would need to solve all 3 equations simultaneously to match the boundary conditions, etc. So I imagine the linear solve $$Ku = f$$ would have to include the components from the $$[m, c, a]$$ components for each point in the mesh. So to do that, would I need to do a Kronecker product to convert the tensor into something I can solve with a linear solve? The details of which solver to use is not important now--iterative, multigrid, etc., but the setup is what I am trying to understand.

If anyone knows any good resources or articles on solving these types of problems using finite elements or FEM for ADRs, please pass them along.

• Is that second equation really correct? Looks like dimensional errors to me with $\dot{c}$ and $c$ with with no constants in-between. Aug 1 at 23:47
• @MikaelÖhman I checked the equations in the book and I copied them correctly. I think the author converts the system to spherical coordinates, but I just showed the equation in the original coordinates. I get what you mean in terms of the dimensional question. It had not even occurred to me. Aug 2 at 0:30
• Some fields of research really like to handwave off dimensional errors like this for some reason. I never understood why. There are several more issues with units here, between $\nabla^2 a$, $m$, $a$ and things like $(1-m)$. It's worse the more i look at it. I think you would be better off trying to find a simpler example that's also a bit more rigorous. Aug 2 at 0:39
• @MikaelÖhman I will look for a better example. I just wanted to confirm though, like it is really about the constants in front of the variables, like m, c, a. If there are constants in front of the variables, then that is a dimensional form which is okay. If there are no constants in front of the variables, then that is non-dimensional from. But if there is a mix of constants and no constants, then that is a mess. Is that a simple heuristic to look for. Aug 2 at 1:14
• No I just do a quick mental dimensional analysis. Especially if the equations come from chemistry or biology. I just fear trying to solve such a messy system which clearly now depends on the choice of time and length scales (ugh) is going to be extra difficult for no good reason. Aug 2 at 1:26

The conceptual framework is to consider your system of equations as a vector of equations, and multiply (dot product) with a vector of test functions. After integration, you then end up with a weak formulation that has just a single number on the left and right hand side.

For the deal.II library, we have summarized this general approach on this page: https://dealii.org/developer/doxygen/deal.II/group__vector__valued.html

• Oh thanks so much for responding. I love your videos :). Yeah, the link you provided is very helpful. I read through it once and then I will go back and work through the details. But it makes sense. The vector idea made sense in principle, but I needed to see some worked out examples to understand the tricks, etc. Thanks again. Aug 4 at 4:32