Nonlinear Poisson equation

In this example, we solve a nonlinear elliptic equation with Newton's method,

−Δu+u4=f, in Ω,u=0, on ∂Ω, -\Delta u + u^4 = f,\qquad\text{ in }\Omega, \\ u = 0,\qquad\text{ on }\partial\Omega,

in the unit square Ω=(0,1)2\Omega=(0,1)^2 with the right-hand side f(x)=−(400∥x∥2−20d) e−10∥x∥2f(\mathbf{x}) = -(400\|\mathbf{x}\|^2 - 20d)\,e^{-10\|\mathbf{x}\|^2} from the Poisson example.

The Poisson and heat examples set up one linear system and solve it. Here, the solution is found by a sequence of linear systems. Each linearized problem is formulated in terms of operators whose coefficients depend on the current iterate. NewtonIteration then re-assembles and solves the system in each step.

Discretization

Newton's method

We write the weak form of the equation as F(u)=0F(u) = 0 with

⟨F(u),ϑ⟩:=∫Ω∇u⋅∇ϑ+u4ϑ−fϑ dx,∀ϑ∈H01(Ω). \langle F(u), \vartheta\rangle := \int_\Omega \nabla u\cdot\nabla\vartheta + u^4\vartheta - f\vartheta\,\text{d}\mathbf{x},\qquad\forall\vartheta\in H^1_0(\Omega).

Starting from an initial guess u0u_0, Newton's method computes a sequence of iterates uk+1=uk+δuu_{k+1} = u_k + \delta u, where the update δu\delta u solves the linear problem DF(uk)[δu]=−F(uk)DF(u_k)[\delta u] = -F(u_k), i.e.,

∫Ω∇δu⋅∇ϑ+4uk3 δu ϑ dx=−∫Ω∇uk⋅∇ϑ dx+∫Ω(f−uk4) ϑ dx,∀ϑ∈H01(Ω), \int_\Omega \nabla\delta u\cdot\nabla\vartheta + 4u_k^3\,\delta u\,\vartheta\,\text{d}\mathbf{x} = -\int_\Omega \nabla u_k\cdot\nabla\vartheta\,\text{d}\mathbf{x} + \int_\Omega \big(f - u_k^4\big)\,\vartheta\,\text{d}\mathbf{x},\qquad\forall\vartheta\in H^1_0(\Omega),

with boundary condition δu=−uk\delta u = -u_k on ∂Ω\partial\Omega, such that uk+1=0u_{k+1} = 0 on the boundary. The iteration stops when the norm of the update ∥δu∥\|\delta u\| is smaller than a given tolerance.

Implementation

Problem definition

We use a YaspGrid created from the initfile parameters with the prefix mesh, and second-order Lagrange elements. The problem is constructed directly from the grid and a dune-functions basis factory:

#include <amdis/AMDiS.hpp>
#include <amdis/AdaptStationary.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/nonlin/NewtonIteration.hpp>

#include <dune/grid/yaspgrid.hh>

using namespace AMDiS;
using namespace Dune::Functions::BasisFactory;

// ...
using Grid = Dune::YaspGrid<GRIDDIM>;
auto gridPtr = MeshCreator<Grid>{"mesh"}.create();

ProblemStat prob("poisson", *gridPtr, lagrange<2>());
prob.initialize(INIT_ALL);

In the Newton iteration, the solution vector of the problem prob.solution() holds the current iterate uku_k. The system matrix holds the Jacobian DF(uk)DF(u_k), and the right-hand side vector holds −F(uk)-F(u_k).

The Newton iteration

The class NewtonIteration wraps a ProblemStat. It implements one Newton step: it assembles the Jacobian and the residual, solves for the update δu\delta u (stored in an extra vector, accessible via newton.stepSolution()), and adds it to the solution vector. The norm of the update is reported as error estimate to the AdaptInfo.

NewtonIteration newton{"newton", prob};

Two parameters can be set in the initfile:

  • newton->build cycle: Re-assemble the Jacobian only in every ii-th iteration, and use the old matrix in between. With 0, the Jacobian is assembled only once. Default: 1.
  • newton->norm: The norm to measure the update, 1 for the L2L^2 norm and any value larger than 1 for the H1H^1 norm. Default: 1.

Jacobian and residual

The bilinear form DF(uk)DF(u_k) consists of the Laplace term and a mass term with a coefficient that depends on the current iterate. The expression 4*pow<3>(prob.solution()) is a grid function that refers to the solution vector. So it is evaluated with the current iterate uku_k each time the operator is assembled:

// <grad(du), grad(theta)>
prob.addMatrixOperator(sot(1.0));
// <4*u_k^3 du, theta>
prob.addMatrixOperator(zot(4*pow<3>(prob.solution())));

The negative residual −F(uk)-F(u_k) is added as vector operators. The term −⟨∇uk,∇ϑ⟩-\langle\nabla u_k, \nabla\vartheta\rangle uses the operator tag::gradtest, which tests a vector-valued coefficient against the gradient of the test functions. gradientOf(...) gives the gradient of a discrete function. The last argument of zot(f, 6) is the quadrature order used for the non-polynomial function f:

// -<grad(u_k), grad(theta)>
prob.addVectorOperator(makeOperator(tag::gradtest{}, -gradientOf(prob.solution())));
// <-u_k^4, theta>
prob.addVectorOperator(zot(-pow<4>(prob.solution())));
// <f, theta>
prob.addVectorOperator(zot(f, 6));

Boundary condition

The Dirichlet value for the update is again an expression of the current iterate. The predicate selects the whole boundary:

prob.addDirichletBC([](auto const&) { return true; }, -prob.solution());

Running the iteration

NewtonIteration implements the same ProblemIterationInterface as ProblemStat. So it can be passed to AdaptStationary, which repeats the iteration until the estimated error is below the tolerance, or until the maximal number of iterations is reached:

AdaptInfo adaptInfo{"adapt"};
AdaptStationary adaptStat{"adapt", newton, adaptInfo};
adaptStat.adapt();

The parameters are read from the initfile. The tolerance refers to the error estimate set by NewtonIteration, i.e., the norm of the Newton update:

adapt->max iteration: 20
adapt[]->tolerance: 1.e-10

newton->build cycle: 2
newton->norm: 2

With build cycle: 2, the Jacobian is re-assembled only in every second step. This makes the iteration cheaper per step, but it may need more steps to converge.

Complete Example

The full source code can be found in the repository at examples/nonlin-poisson.cc.

Compile with 

cmake --build [build-dir] --target nonlin-poisson.2d

and run in the examples build directory with 

./nonlin-poisson.2d [amdis-dir]/examples/init/nonlin-poisson.dat.2d