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

Question

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.

Answer 1

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

Answer 2

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.