The problems we have encountered so far—with the notable exception
of the Navier–Stokes equations—all share a common feature: they all
involve models expressed by a *single* scalar or vector PDE. In many
situations the model is instead expressed as a system of PDEs,
describing different quantities possibly governed by (very) different
physics. As we saw for the Navier–Stokes equations, one way to solve
a system of PDEs in FEniCS is to use a splitting method where we solve
one equation at a time and feed the solution from one equation into
the next. However, one of the strengths with FEniCS is the ease by
which one can instead define variational problems that couple several
PDEs into one compound system. In this section, we will look at how to use
FEniCS to write solvers for such systems of coupled PDEs.
The goal is to demonstrate how easy it is to implement fully implicit,
also known as monolithic, solvers in FEniCS.

Our model problem is the following system of advection–diffusion–reaction equations: $$ \begin{align} \tag{3.37} \frac{\partial u_1}{\partial t} + w \cdot \nabla u_1 - \nabla\cdot(\epsilon\nabla u_1) &= f_1 -K u_1 u_2, \\ \tag{3.38} \frac{\partial u_2}{\partial t} + w \cdot \nabla u_2 - \nabla\cdot(\epsilon\nabla u_2) &= f_2 -K u_1 u_2, \\ \tag{3.39} \frac{\partial u_3}{\partial t} + w \cdot \nabla u_3 - \nabla\cdot(\epsilon\nabla u_3) &= f_3 + K u_1 u_2 - K u_3. \end{align} $$

This system models the chemical reaction between two species \( A \) and \( B \) in some domain \( \Omega \): $$ A + B \rightarrow C. $$

We assume that the reaction is *first-order*, meaning that the
reaction rate is proportional to the concentrations \( [A] \) and \( [B] \) of
the two species \( A \) and \( B \):
$$
\frac{\mathrm{d}}{\mathrm{d}t} [C] = K [A] [B].
$$
We also assume that the formed species \( C \) spontaneously decays with a
rate proportional to the concentration \( [C] \). In the PDE system
(3.37)--(3.39),
we use the variables \( u_1 \), \( u_2 \), and \( u_3 \) to denote the
concentrations of the three species:
$$
u_1 = [A], \quad u_2 = [B], \quad u_3 = [C].
$$
We see that the chemical reactions are accounted for in the
right-hand sides of the PDE system (3.37)--(3.39).

The chemical reactions take part at each point in the domain \( \Omega \). In addition, we assume that the species \( A \), \( B \), and \( C \) diffuse throughout the domain with diffusivity \( \epsilon \) (the terms \( -\nabla\cdot(\epsilon\nabla u_i) \)) and are advected with velocity \( w \) (the terms \( w\cdot\nabla u_i \)).

To make things visually and physically interesting, we shall let the chemical reaction take place in the velocity field computed from the solution of the incompressible Navier–Stokes equations around a cylinder from the previous section. In summary, we will thus be solving the following coupled system of nonlinear PDEs: $$ \begin{align} \tag{3.40} \varrho\left(\frac{\partial w}{\partial t} + w \cdot \nabla w\right) &= \nabla\cdot\sigma(w, p) + f, \\ \nabla \cdot w &= 0, \tag{3.41}\\ \frac{\partial u_1}{\partial t} + w \cdot \nabla u_1 - \nabla\cdot(\epsilon\nabla u_1) &= f_1 - K u_1 u_2, \tag{3.42}\\ \frac{\partial u_2}{\partial t} + w \cdot \nabla u_2 - \nabla\cdot(\epsilon\nabla u_2) &= f_2 - K u_1 u_2, \tag{3.43}\\ \frac{\partial u_3}{\partial t} + w \cdot \nabla u_3 - \nabla\cdot(\epsilon\nabla u_3) &= f_3 + K u_1 u_2 - K u_3. \tag{3.44} \end{align} $$ We assume that \( u_1 = u_2 = u_3 = 0 \) at \( t = 0 \) and inject the species \( A \) and \( B \) into the system by specifying nonzero source terms \( f_1 \) and \( f_2 \) close to the corners at the inflow, and take \( f_3 = 0 \). The result will be that \( A \) and \( B \) are convected by advection and diffusion throughout the channel, and when they mix the species \( C \) will be formed.

Since the system is one-way coupled from the Navier–Stokes subsystem
to the advection–diffusion–reaction subsystem, we do not need to
recompute the solution to the Navier–Stokes equations, but can just
read back the previously computed velocity field \( w \) and feed it into
our equations. But we *do* need to learn how to read and write
solutions for time-dependent PDE problems.

We obtain the variational formulation of our system by multiplying each equation by a test function, integrating the second-order terms \( -\nabla\cdot(\epsilon\nabla u_i) \) by parts, and summing up the equations. When working with FEniCS it is convenient to think of the PDE system as a vector of equations. The test functions are collected in a vector too, and the variational formulation is the inner product of the vector PDE and the vector test function.

We also need introduce some discretization in time. We will use the backward Euler method as before when we solved the heat equation and approximate the time derivatives by \( (u_i^{n+1}-u_i^n) / \dt \). Let \( v_1 \), \( v_2 \), and \( v_3 \) be the test functions, or the components of the test vector function. The inner product results in $$ \begin{align} \tag{3.45} & \int_{\Omega} (\dt^{-1} (u_1^{n+1} - u_1^n) v_1 + w \cdot \nabla u^{n+1}_1 \, v_1 + \epsilon \nabla u^{n+1}_1 \cdot \nabla v_1) \dx \\ + & \int_{\Omega} (\dt^{-1} (u_2^{n+1} - u_2^n) v_2 + w \cdot \nabla u^{n+1}_2 \, v_2 + \epsilon \nabla u^{n+1}_2 \cdot \nabla v_2) \dx \nonumber \\ + & \int_{\Omega} (\dt^{-1} (u_3^{n+1} - u_3^n) v_3 + w \cdot \nabla u^{n+1}_3 \, v_3 + \epsilon \nabla u^{n+1}_3 \cdot \nabla v_3) \dx \nonumber \\ - & \int_{\Omega} (f_1 v_1 + f_2 v_2 + f_3 v_3) \dx \nonumber \\ - & \int_{\Omega} (-K u^{n+1}_1 u^{n+1}_2 v_1 - K u^{n+1}_1 u^{n+1}_2 v_2 + K u^{n+1}_1 u^{n+1}_2 v_3 - K u^{n+1}_3 v_3) \dx = 0. \nonumber \end{align} $$ For this problem it is natural to assume homogeneous Neumann boundary conditions on the entire boundary for \( u_1 \), \( u_2 \), and \( u_3 \); that is, \( \partial u_i/\partial n = 0 \) for \( i = 1, 2, 3 \). This means that the boundary terms vanish when we integrate by parts.

The first step is to read the mesh from file. Luckily, we made sure to save the mesh to file in the Navier–Stokes example and can now easily read it back from file:

```
mesh = Mesh('navier_stokes_cylinder/cylinder.xml.gz')
```

The mesh is stored in the native FEniCS XML format (with additional gzipping to decrease the file size).

Next, we need to define the finite element function space. For this problem, we need to define several spaces. The first space we create is the space for the velocity field \( w \) from the Navier–Stokes simulation. We call this space \( W \) and define the space by

```
W = VectorFunctionSpace(mesh, 'P', 2)
```

It is important that this space is exactly the same as the space we
used for the velocity field in the Navier–Stokes solver. To read the
values for the velocity field, we use a `TimeSeries`

:

```
timeseries_w = TimeSeries('navier_stokes_cylinder/velocity_series')
```

This will initialize the object `timeseries_w`

which we will call
later in the time-stepping loop to retrieve values from the
file `velocity_series.h5`

(in binary HDF5 format).

For the three concentrations \( u_1 \), \( u_2 \), and \( u_3 \), we want to
create a *mixed space* with functions that represent the full system
\( (u_1, u_2, u_3) \) as a single entity. To do this, we need to define a
`MixedElement`

as the product space of three simple finite elements
and then used the mixed element to define the function space:

```
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1, P1])
V = FunctionSpace(mesh, element)
```

```
P2 = VectorElement('P', triangle, 2)
P1 = FiniteElement('P', triangle, 1)
TH = P2 * P1
```

This syntax works great for two elements, but for three or more
elements we meet a subtle issue in how the Python interpreter handles
the `*`

operator. For the reaction system, we create the mixed element
by `element = MixedElement([P1, P1, P1])`

and one would be tempted to
write

```
element = P1 * P1 * P1
```

However, this is equivalent to writing `element = (P1 * P1) * P1`

so
the result will be a mixed element consisting of two subsystems, the
first of which in turn consists of two scalar subsystems.

Finally, we remark that for the simple case of a mixed system consisting of three scalar elements as for the reaction system, the definition is in fact equivalent to using a standard vector-valued element:

```
element = VectorElement('P', triangle, 1, dim=3)
V = FunctionSpace(mesh, element)
```

Once the space has been created, we need to define our test functions
and finite element functions. Test functions for a mixed function
space can be created by replacing `TestFunction`

by `TestFunctions`

:

```
v_1, v_2, v_3 = TestFunctions(V)
```

Since the problem is nonlinear, we need to work with functions rather
than trial functions for the unknowns. This can be done by using the
corresponding `Functions`

construction in FEniCS. However, as we will
need to access the `Function`

for the entire system itself, we first
need to create that function and then access its components:

```
u = Function(V)
u_1, u_2, u_3 = split(u)
```

These functions will be used to represent the unknowns \( u_1 \), \( u_2 \), and \( u_3 \)
at the new time level \( n+1 \). The corresponding values at the previous
time level \( n \) are denoted by `u_n1`

, `u_n2`

, and `u_n3`

in our program.

When now all functions and test functions have been defined, we can express the nonlinear variational problem (3.45):

```
F = ((u_1 - u_n1) / k)*v_1*dx + dot(w, grad(u_1))*v_1*dx \
+ eps*dot(grad(u_1), grad(v_1))*dx + K*u_1*u_2*v_1*dx \
+ ((u_2 - u_n2) / k)*v_2*dx + dot(w, grad(u_2))*v_2*dx \
+ eps*dot(grad(u_2), grad(v_2))*dx + K*u_1*u_2*v_2*dx \
+ ((u_3 - u_n3) / k)*v_3*dx + dot(w, grad(u_3))*v_3*dx \
+ eps*dot(grad(u_3), grad(v_3))*dx - K*u_1*u_2*v_3*dx + K*u_3*v_3*dx \
- f_1*v_1*dx - f_2*v_2*dx - f_3*v_3*dx
```

The time-stepping simply consists of solving this variational problem
in each time step by a call to the `solve`

function:

```
t = 0
for n in range(num_steps):
t += dt
timeseries_w.retrieve(w.vector(), t)
solve(F == 0, u)
u_n.assign(u)
```

In each time step, we first read the current value for the velocity
field from the time series we have previously stored. We then solve
the nonlinear system, and assign the computed values to the left-hand
side values for the next time interval. When retrieving values from a
`TimeSeries`

, the values will by default be interpolated (linearly) to
the given time `t`

if the time does not exactly match a sample in the
series.

The solution at the final time is shown in Figure 11. We clearly see the advection of the species \( A \) and \( B \) and the formation of \( C \) along the center of the channel where \( A \) and \( B \) meet.

The complete code is presented below.

```
from fenics import *
T = 5.0 # final time
num_steps = 500 # number of time steps
dt = T / num_steps # time step size
eps = 0.01 # diffusion coefficient
K = 10.0 # reaction rate
# Read mesh from file
mesh = Mesh('navier_stokes_cylinder/cylinder.xml.gz')
# Define function space for velocity
W = VectorFunctionSpace(mesh, 'P', 2)
# Define function space for system of concentrations
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1, P1])
V = FunctionSpace(mesh, element)
# Define test functions
v_1, v_2, v_3 = TestFunctions(V)
# Define functions for velocity and concentrations
w = Function(W)
u = Function(V)
u_n = Function(V)
# Split system functions to access components
u_1, u_2, u_3 = split(u)
u_n1, u_n2, u_n3 = split(u_n)
# Define source terms
f_1 = Expression('pow(x[0]-0.1,2)+pow(x[1]-0.1,2)<0.05*0.05 ? 0.1 : 0',
degree=1)
f_2 = Expression('pow(x[0]-0.1,2)+pow(x[1]-0.3,2)<0.05*0.05 ? 0.1 : 0',
degree=1)
f_3 = Constant(0)
# Define expressions used in variational forms
k = Constant(dt)
K = Constant(K)
eps = Constant(eps)
# Define variational problem
F = ((u_1 - u_n1) / k)*v_1*dx + dot(w, grad(u_1))*v_1*dx \
+ eps*dot(grad(u_1), grad(v_1))*dx + K*u_1*u_2*v_1*dx \
+ ((u_2 - u_n2) / k)*v_2*dx + dot(w, grad(u_2))*v_2*dx \
+ eps*dot(grad(u_2), grad(v_2))*dx + K*u_1*u_2*v_2*dx \
+ ((u_3 - u_n3) / k)*v_3*dx + dot(w, grad(u_3))*v_3*dx \
+ eps*dot(grad(u_3), grad(v_3))*dx - K*u_1*u_2*v_3*dx + K*u_3*v_3*dx \
- f_1*v_1*dx - f_2*v_2*dx - f_3*v_3*dx
# Create time series for reading velocity data
timeseries_w = TimeSeries('navier_stokes_cylinder/velocity_series')
# Create VTK files for visualization output
vtkfile_u_1 = File('reaction_system/u_1.pvd')
vtkfile_u_2 = File('reaction_system/u_2.pvd')
vtkfile_u_3 = File('reaction_system/u_3.pvd')
# Create progress bar
progress = Progress('Time-stepping')
set_log_level(PROGRESS)
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Read velocity from file
timeseries_w.retrieve(w.vector(), t)
# Solve variational problem for time step
solve(F == 0, u)
# Save solution to file (VTK)
_u_1, _u_2, _u_3 = u.split()
vtkfile_u_1 << (_u_1, t)
vtkfile_u_2 << (_u_2, t)
vtkfile_u_3 << (_u_3, t)
# Update previous solution
u_n.assign(u)
# Update progress bar
progress.update(t / T)
# Hold plot
interactive()
```

This example program can be found in the file `ft09_reaction_system.py`.

Finally, we comment on three important techniques that are very useful when working with systems of PDEs: setting initial conditions, setting boundary conditions, and extracting components of the system for plotting or postprocessing.

In our example, we did not need to worry about setting an initial
condition, since we start with \( u_1 = u_2 = u_3 = 0 \). This happens
automatically in the code when we set `u_n = Function(V)`

. This
creates a `Function`

for the whole system and all degrees of freedom
are set to zero.

If we want to set initial conditions for the components of the system
separately, the easiest solution is to define the initial conditions
as a vector-valued `Expression`

and then project (or interpolate) this
to the `Function`

representing the whole system. For example,

```
u_0 = Expression(('sin(x[0])', 'cos(x[0]*x[1])', 'exp(x[1])'), degree=1)
u_n = project(u_0, V)
```

This defines \( u_1 \), \( u_2 \), and \( u_2 \) to be the projections of \( \sin x \), \( \cos (xy) \), and \( \exp(y) \), respectively.

In our example, we also did not need to worry about setting boundary
conditions since we used a natural Neumann condition. If we want to set
Dirichlet conditions for individual components of the system, this can
be done as usual by the class `DirichletBC`

, but we must specify for
which subsystem we set the boundary condition. For example, to specify
that \( u_2 \) should be equal to \( xy \) on the boundary defined by
`boundary`

, we do

```
u_D = Expression('x[0]*x[1]', degree=1)
bc = DirichletBC(V.sub(1), u_D, boundary)
```

The object `bc`

or a list of such objects containing different
boundary conditions, can then be passed to the `solve`

function as usual.
Note that numbering starts at \( 0 \) in FEniCS so the subspace
corresponding to \( u_2 \) is `V.sub(1)`

.

If `u`

is a `Function`

defined on a mixed function space in FEniCS,
there are several ways in which `u`

can be *split* into components.
Above we already saw an example of the first of these:

```
u_1, u_2, u_3 = split(u)
```

This extracts the components of `u`

as *symbols* that can be used in a
variational problem. The above statement is in fact equivalent to

```
u_1 = u[0]
u_2 = u[1]
u_3 = u[2]
```

Note that `u[0]`

is not really a `Function`

object, but merely a
symbolic expression, just like `grad(u)`

in FEniCS is a symbolic
expression and not a `Function`

representing the gradient. This means
that `u_1`

, `u_2`

, `u_3`

can be used in a variational problem, but
cannot be used for plotting or postprocessing.

To access the components of `u`

for plotting and saving the solution
to file, we need to use a different variant of the `split`

function:

```
u_1_, u_2_, u_3_ = u.split()
```

This returns three subfunctions as actual objects with access to the
common underlying data stored in `u`

, which makes plotting and saving
to file possible. Alternatively, we can do

```
u_1_, u_2_, u_3_ = u.split(deepcopy=True)
```

which will create `u_1_`

, `u_2_`

, and `u_3_`

as stand-alone `Function`

objects, each holding a copy of the subfunction data extracted from
`u`

. This is useful in many situations but is not necessary for
plotting and saving solutions to file.