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:

\begin{array}{r} \begin{aligned} - η Δ u + α u (1 - u) = f & in Ω \\ u = g & on \partial Ω \end{aligned} \end{array}

with $Ω = [0, 1]^{2}$ . Being the problem non-linear, to find a solution (which esists unique for $Ω$ Lipschitz and $η > 0, α > 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 $L (u) = - η Δ u + α u (1 - u) - f$ , the Netwon scheme reads as follow:

\begin{array}{r} \begin{matrix} D L_{u^{k}} (h) = - L (u^{k}) \\ u^{k + 1} = u^{k} + h \end{matrix} \end{array}

being $D L_{u^{k}} (h)$ the Gateaux derivative of the operator $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):

D L_{u^{k}} (h) = - η Δ h + α (1 - 2 u^{k}) h

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

- η Δ h + α (1 - 2 u^{k}) h = f - (- η Δ u^{k} + α u^{k} (1 - u^{k})) .

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

- η Δ u^{k + 1} + α (1 - 2 u^{k}) u^{k + 1} = f - α (u^{k})^{2}

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

\int_{Ω} \nabla u^{k + 1} \cdot \nabla v + α (1 - u^{k}) u^{k + 1} v - \int_{Ω} α u^{k} u^{k + 1} v = \int_{Ω} f - α (u^{k})^{2} v \forall v \in V_{h}

which we iteratively solve for $u^{k + 1}$ until convergence.

Implementation¶

The first step in any finite element code is the definition of the problem geometry. As we are going to iterate several times over the mesh is convinient to compute the cells once and cache them to obtain a fast re-access to the mesh. This option is enabled by activating the cell caching.

Triangulation<2, 2> unit_square = Triangulation<2, 2>::UnitSquare(60, cache_cells);

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

FeSpace Vh(unit_square, P1<1>); // piecewise linear continuous scalar 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 (p)$ as $\sum_{i = 1}^{N} c_{i} ψ_{i} (p)$ given the expansion coefficient vector $c = [c_{1}, \dots, 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

- \int_{Ω} α (u^{k})^{2} v = \sum_{i = 1}^{N} \sum_{j = 1}^{N} [\int_{Ω} - α u^{k} ψ_{i} ψ_{j}] c_{i} = A_{2} c

being $A_{2}$ the discretization matrix of the bilinear form $b (u, v) = - \int_{Ω} α 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_{Ω} f v$ .

ScalarField<2, decltype([](const PointT& 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, QS2DP4)(f * v);

Observe that we explicitly require an higher order quadrature specifying the quadrature formula QS2DP4 for the exact integration of order 4 polynomials as second argument of the integral function. Finally, we define non-homegeneous Dirichlet boundary conditions $g (x) = 3 x^{2} + 2 y^{2}$ on all the boundary of the domain

ScalarField<2, decltype([](const PointT& p) { return 3 * p[0] * p[0] + 2 * p[1] * p[1]; })> g;
auto& 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 c = b$ deriving form the discretization of the variational problem, and that we must later enforce them on the pair $(A, 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 = Eigen::Matrix<double, Dynamic, 1>::Zero(Vh.n_dofs());   // initial guess u = 0
Eigen::SparseMatrix<double> A = a.assemble();
Eigen::Matrix<double, Dynamic, 1> v_ = F.assemble();
dof_handler.enforce_constraints(A, v_);
// linear system solve A*u_prev = v_ using Cholesky factorization
Eigen::SimplicialLLT<Eigen::SparseMatrix<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_{Ω} \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

\int_{Ω} \nabla u^{k + 1} \cdot \nabla v + α (1 - u^{k}) u^{k + 1} v - \int_{Ω} α u^{k} u^{k + 1} v = \int_{Ω} f - α (u^{k})^{2} v \forall v \in V_{h}

while (err > 1e-7) {
   Eigen::SparseMatrix<double> A1 = a.assemble();
   Eigen::SparseMatrix<double> A2 = b.assemble();
   Eigen::Matrix<double, Dynamic, 1> 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_{Ω} f - α (u^{k})^{2} v$ , enforces the Dirichlet constaints on the resulting linear system and solves the resulting linear system.

The complete script

#include <fdaPDE/finite_elements.h>
using namespace fdapde;

int main() {
   // useful typedef and constants definition
   constexpr int local_dim = 2;
   using PointT = Eigen::Matrix<double, local_dim, 1>;

   Triangulation<local_dim, local_dim> unit_square = Triangulation<2, 2>::UnitSquare(60, cache_cells);

   FeSpace Vh(unit_square, P1<1>);
   // 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 PointT& 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, QS2DP4)(f * v);

   // define dirichlet data
   ScalarField<2, decltype([](const PointT& p) { return 3 * p[0] * p[0] + 2 * p[1] * p[1]; })> g;
   auto& 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 = Eigen::Matrix<double, Dynamic, 1>::Zero(Vh.n_dofs());   // initial guess u = 0
   Eigen::SparseMatrix<double> A = a.assemble();   // this actually assembles dot(grad(u), grad(v)) + u * v
   Eigen::Matrix<double, Dynamic, 1> v_ = F.assemble();
   dof_handler.enforce_constraints(A, v_);
   // linear system solve A*u_prev = v_ using Cholesky factorization
   Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> lin_solver(A);
   u_prev = lin_solver.solve(v_);

   double err = std::numeric_limits<double>::max();
   while (err > 1e-7) {
       Eigen::SparseMatrix<double> A1 = a.assemble();
       Eigen::SparseMatrix<double> A2 = b.assemble();
       Eigen::Matrix<double, Dynamic, 1> 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;
}

[GT23]

Carlotta Gatti and Emanuele Tamburini. Solving nonlinear elliptic equations with femr. Advanced Programming for Scientific Computing project, 2023.