Example 1: A non-linear diffusion reaction problem¶

In this example we solve a non-linear diffusion reaction problem, corresponding to the stationary case of the Fisher KPP equation. Specifically, the problem is described by the following equation:

with \(\Omega = [0,1]^2\). Being the problem non-linear, to find a solution (which esists unique for \(\Omega\) Lipschitz and \(\eta > 0, \alpha > 0\)) we need some iterative scheme which converges, with some reasonable degree of accuracy, to the solution. We here develop a Newton scheme following the work done in [GT23]. Let \(\mathcal{L}(u) = -\eta \Delta u + \alpha u(1-u) - f\), the Netwon scheme reads as follow:

being \(D \mathcal{L}_{u^k}(h)\) the Gateaux derivative of the operator \(\mathcal{L}\) in the direction \(h\). Computations leads to the following expression for the Gateaux derivative (it is not the point of the following example how to compute such functional derivative):

Substituting the above expression in the Newton formula we obtain the following linearized problem to be solved for \(h\)

Finally, observing that by definition of Newton update, \(h = u^{k+1} - u^k\), we get the more computationally efficient iteration

As always, fdaPDE requires the weak formulation of the differential problem above. Following the Galerkin approach, let \(V_h\) a proper finite dimensional subspace of the space where we seek the solution, and let \(u^{k+1}, u^k, v \in V_h\). Multiplying both sides of the last equation by the test function \(v\) and integrating we recover the following variational identity

which we iteratively solve for \(u^{k+1}\) until convergence. Follows a step by step description of the program which enables us to find a solution to the considered problem in less than 50 lines of code.

First we load the geometry (currently generated from some third party triangulator), enabling the cell caching. As we are going to iterate several times over the mesh is convinient to compute the cells once to obtain a fast re-access to the mesh.

Triangulation<2, 2> unit_square = read_mesh<2, 2>("unit_square", cache_cells);

Once we have a discretization of the domain \(\Omega\), we can instatiate a (linear) finite element space on it togheter with trial and test functions.

FiniteElementSpace Vh(unit_square, P1); // P1 denotes the space of linear finite elements
TrialFunction u(Vh);
TestFunction  v(Vh);

Because we are dealing with a non-linear problem, we need an object describing the current and previous solution \(u^{k+1}\) and \(u^k\). This is achieved using the FeFunction type, which represents a scalar field \(u\) written with respect to a finite element basis system, e.g. it represents \(u(\boldsymbol{p})\) as \(\sum_{i=1}^N c_i \psi_i(\boldsymbol{p})\) given the expansion coefficient vector \(\boldsymbol{c} = [c_1, \ldots, c_N]\).

FeFunction u_prev(Vh), solution(Vh);

We now have all the ingredients to assemble our bilinear forms. First we notice that we can save some computational time observing that

being \(A_2\) the discretization matrix of the bilinear form \(b(u, v) = - \int_{\Omega} \alpha u^k u v\), which can therefore be computed once for the lhs and rhs. For this reason, we defined two bilinear forms as follow (look how the writing is almost equivalent to the mathematical one):

auto a = integral(unit_square)(dot(grad(u), grad(v)) + (1 - u_prev) * u * v);
auto b = integral(unit_square)(-u_prev * u * v);

We then define the forcing term as a plain ScalarField togheter with the forcing functional \(F(v) = \int_{\Omega} f v\).

ScalarField<2, decltype([](const SVector<2>& p) {
    return -9*std::pow(p[0], 4) - 12*p[0]*p[0]*p[1]*p[1] + 3*p[0]*p[0] +
            2*p[1]*p[1] - 4*std::pow(p[1], 4) - 10;
})> f;
auto F = integral(unit_square, QS2D6P)(f * v);

Observe that we explicitly require an higher order quadrature specifying the 6 points quadrature formula QS2D6P as second argument of the integral function. Finally, we define non-homegeneous Dirichlet boundary conditions \(g(\boldsymbol{x}) = 3x^2 + 2y^2\) on all the boundary of the domain

ScalarField<2, decltype([](const SVector<2>& p) { return 3 * p[0] * p[0] + 2 * p[1] * p[1]; })> g;
DofHandler<2, 2>& dof_handler = Vh.dof_handler();
dof_handler.set_dirichlet_constraint(/* on = */ BoundaryAll, /* data = */ g);

Recall that Dirichlet boundary conditions are implemented as constraints on the degrees of freedom of the linear system \(A \boldsymbol{c} = \boldsymbol{b}\) deriving form the discretization of the variational problem, and that we must later enforce them on the pair \((A, \boldsymbol{b})\) before solving the linear system, using the enforce_constraints method.

We can now find an initial point for the Newton scheme. To this end, we solve the linerized problem with initial guess \(u = 0\).

u_prev = DVector<double>::Zero(Vh.n_dofs());   // initial guess u = 0
SpMatrix<double> A = a.assemble();
DVector<double> v_ = F.assemble();
dof_handler.enforce_constraints(A, v_);
// linear system solve A*u_prev = v_ using Cholesky factorization
Eigen::SimplicialLLT<SpMatrix<double>> lin_solver(A);
u_prev = lin_solver.solve(v_);

The code fragment above effectivelly assemble the discretization matrix A for the bilinear form \(\int_{\Omega} \nabla u^0 \cdot \nabla v + u^0 v\) togheter with the discretizing vector v_ of the forcing functional \(F\). Then, it sets the Dirichlet conditions at the boundary via the enforce_constaints method of the dof_handler object. Finally, observing that the bilinear form is SPD, solves the FEM linear system using a Cholesky factorization and sets \(u^0\) to the solution of this linear system.

We can finally start looping until convergence, iteratively solving the recurrence

while (err > 1e-7) {
   SpMatrix<double> A1 = a.assemble();
   SpMatrix<double> A2 = b.assemble();
   DVector<double> v = v_ + A2 * u_prev.coeff();
   A = A1 + A2;
   dof_handler.enforce_constraints(A, v);
   lin_solver.compute(A);
   solution = lin_solver.solve(v);

   // update error
   err = (u_prev.coeff() - solution.coeff()).norm();
   u_prev = solution;
}

The code just assembles A1 and A2, updates the right hand side \(\int_{\Omega} f - \alpha (u^k)^2 v\), enforces the Dirichlet constaints on the resulting linear system and solves the resulting linear system.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

#include <fdaPDE/fields.h>
#include <fdaPDE/geometry.h>
#include <fdaPDE/finite_elements.h>

using namespace fdapde;

int main() {
   // import mesh, enable cell caching for fast re-cycling
   Triangulation<2, 2> unit_square = read_mesh<2, 2>("unit_square", cache_cells);

   FiniteElementSpace Vh(unit_square, P1);
   // create trial and test functions
   TrialFunction u(Vh);
   TestFunction  v(Vh);
   // current solution and solution at previous step
   FeFunction u_prev(Vh), solution(Vh);

   // define bilinear forms
   auto a = integral(unit_square)(dot(grad(u), grad(v)) + (1 - u_prev) * u * v);
   auto b = integral(unit_square)(-u_prev * u * v);

   // define forcing functional
   ScalarField<2, decltype([](const SVector<2>& p) {
       return -9*std::pow(p[0], 4) - 12*p[0]*p[0]*p[1]*p[1] + 3*p[0]*p[0] + 2*p[1]*p[1] - 4*std::pow(p[1], 4) - 10;
   })> f;
   auto F = integral(unit_square, QS2D6P)(f * v);

   // define dirichlet data
   ScalarField<2, decltype([](const SVector<2>& p) { return 3 * p[0] * p[0] + 2 * p[1] * p[1]; })> g;
   DofHandler<2, 2>& dof_handler = Vh.dof_handler();
   dof_handler.set_dirichlet_constraint(/* on = */ BoundaryAll, /* data = */ g);

   // Newton scheme initialization (solve linearized problem with initial guess u = 0)
   u_prev = DVector<double>::Zero(Vh.n_dofs());   // initial guess u = 0
   SpMatrix<double> A = a.assemble();   // this actually assembles dot(grad(u), grad(v)) + u * v
   DVector<double> v_ = F.assemble();
   dof_handler.enforce_constraints(A, v_);
   // linear system solve A*u_prev = v_ using Cholesky factorization
   Eigen::SimplicialLLT<SpMatrix<double>> lin_solver(A);
   u_prev = lin_solver.solve(v_);

   double err = std::numeric_limits<double>::max();
   while (err > 1e-7) {
       SpMatrix<double> A1 = a.assemble();
       SpMatrix<double> A2 = b.assemble();
       DVector<double> v = v_ + A2 * u_prev.coeff();    // update rhs
       A = A1 + A2;
       dof_handler.enforce_constraints(A, v);
       lin_solver.compute(A);
       solution = lin_solver.solve(v);

       // update error
       err = (u_prev.coeff() - solution.coeff()).norm();
       u_prev = solution;
   }
   return 0;
}