# 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. Commented Aug 1, 2023 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. Commented Aug 2, 2023 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. Commented Aug 2, 2023 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. Commented Aug 2, 2023 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. Commented Aug 2, 2023 at 1:26

## 2 Answers

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. Commented Aug 4, 2023 at 4:32

To answer your questions:

1. Yes, multiply the equations by test functions to get a weak form. This should give you a set of ordinary differential equations (ODEs) in time. For each variable, there would be as many ODEs as nodes for that variable.

2. To begin with, I wouldn't do what you have suggested. Since the right hand side is non-linear, you'd need Newton's iterations to solve, which is complex. To begin with, just use the time-lagged quantities on the right hand side, to get an explicit method. Don't worry about the stable time-step right now, keep cutting it down until you get a solution that does not blow up. Once you have something that does not blow up, I would compare it to a solution obtained via Newton's iterations.

• thanks for the help. I am getting it now. I missed it before, but now I can see that the first equation has two nonlinear terms, which would call for using Newton's method as you say. Could you explain your suggestion for an explicit method a little more. So do you suggest like using a fixed time-stepping explicit scheme or multi-step scheme to get around the extra time involved for the nonlinear solve. Like experimenting with simple RK45 or Adams-Bashforth at first. I understand the idea in theory, but have not tried to solve nonlinear problems with explicit methods. Commented Aug 2, 2023 at 16:05
• Like I would probably need to setup a simple toy nonlinear problem and then try to experiment with the explicit versus implicit solvers, to get a feel for this timestep tuning. If you know of an article or even a SO post that played with something like this, please send it along. Commented Aug 2, 2023 at 16:07
• Your nonlinearities look pretty tame unless the coefficients are enormous. The big hurdle for explicit methods will be getting the diffusion stable. If you want to jump straight into Newton, you could try solving the steady state equations first to get used to it Commented Aug 2, 2023 at 19:18
• @krishnab Take a look at Hoffman and Chiang's CFD book. It's got a good introduction to explicit schemes.
– NNN
Commented Aug 3, 2023 at 3:22
• @NNN yes, I will take a look at that book. Thanks for the tip. It will be interesting to experiment with this approach, I am kinda excited about it. Commented Aug 3, 2023 at 6:56