$$\newcommand{\dt}{\Delta t} \newcommand{\tp}{\thinspace .} \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\x}{\boldsymbol{x}} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\uI}{u_{_0}} \newcommand{\ub}{u_{_\mathrm{D}}} \newcommand{\GD}{\Gamma_{_\mathrm{D}}} \newcommand{\GN}{\Gamma_{_\mathrm{N}}} \newcommand{\GR}{\Gamma_{_\mathrm{R}}} \newcommand{\inner}{\langle #1, #2 \rangle}$$

# A nonlinear Poisson equation

We shall now address how to solve nonlinear PDEs. We will see that nonlinear problems can be solved just as easily as linear problems in FEniCS, by simply defining a nonlinear variational problem and calling the solve function. When doing so, we will encounter a subtle difference in how the variational problem is defined.

## PDE problem

As a model problem for the solution of nonlinear PDEs, we take the following nonlinear Poisson equation: $$\begin{equation} -\nabla\cdot\left(q(u)\nabla u\right) = f, \tag{3.16} \end{equation}$$ in $$\Omega$$, with $$u=\ub$$ on the boundary $$\partial\Omega$$. The coefficient $$q = q(u)$$ makes the equation nonlinear (unless $$q(u)$$ is constant in $$u$$).

## Variational formulation

As usual, we multiply our PDE by a test function $$v\in\hat V$$, integrate over the domain, and integrate the second-order derivatives by parts. The boundary integral arising from integration by parts vanishes wherever we employ Dirichlet conditions. The resulting variational formulation of our model problem becomes: find $$u \in V$$ such that $$\begin{equation} F(u; v) = 0 \quad \forall v \in \hat{V}, \tag{3.17} \end{equation}$$ where $$\begin{equation} F(u; v) = \int_\Omega (q(u)\nabla u\cdot \nabla v - fv) \dx, \tag{3.18} \end{equation}$$ and \begin{align*} V &= \{v \in H^1(\Omega) : v = \ub \mbox{ on } \partial\Omega\},\\ \hat{V} &= \{v \in H^1(\Omega) : v = 0 \mbox{ on } \partial\Omega\}\tp \end{align*}

The discrete problem arises as usual by restricting $$V$$ and $$\hat V$$ to a pair of discrete spaces. As before, we omit any subscript on the discrete spaces and discrete solution. The discrete nonlinear problem is then written as: find $$u\in V$$ such that $$\begin{equation} F(u; v) = 0 \quad \forall v \in \hat{V}, \tag{3.19} \end{equation}$$ with $$u = \sum_{j=1}^N U_j \phi_j$$. Since $$F$$ is nonlinear in $$u$$, the variational statement gives rise to a system of nonlinear algebraic equations in the unknowns $$U_1,\ldots,U_N$$.

## FEniCS implementation

### Test problem

To solve a test problem, we need to choose the right-hand side $$f$$, the coefficient $$q(u)$$ and the boundary value $$\ub$$. Previously, we have worked with manufactured solutions that can be reproduced without approximation errors. This is more difficult in nonlinear problems, and the algebra is more tedious. However, we may utilize SymPy for symbolic computing and integrate such computations in the FEniCS solver. This allows us to easily experiment with different manufactured solutions. The forthcoming code with SymPy requires some basic familiarity with this package. In particular, we will use the SymPy functions diff for symbolic differentiation and ccode for C/C++ code generation.

We take $$q(u) = 1 + u^2$$ and define a two-dimensional manufactured solution that is linear in $$x$$ and $$y$$:

# Warning: from fenics import * will import both sym and
# q from FEniCS. We therefore import FEniCS first and then
# overwrite these objects.
from fenics import *

def q(u):
"Return nonlinear coefficient"
return 1 + u**2

# Use SymPy to compute f from the manufactured solution u
import sympy as sym
x, y = sym.symbols('x, x')
u = 1 + x + 2*y
f = - sym.diff(q(u)*sym.diff(u, x), x) - sym.diff(q(u)*sym.diff(u, y), y)
f = sym.simplify(f)
u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)
print('u =', u_code)
print('f =', f_code)


Define symbolic coordinates as required in Expression objects. Note that we would normally write x, y = sym.symbols('x, y'), but if we want the resulting expressions to have valid syntax for FEniCS Expression objects, we must use x and x. This is easily accomplished with sympy by defining the names of x and y as x and x: x, y = sym.symbols('x, x').

Turning the expressions for u and f into C or C++ syntax for FEniCS Expression objects needs two steps. First, we ask for the C code of the expressions:

u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)


In some cases, one will need to edit the result to match the required syntax of Expression objects, but not in this case. (The primary example is that M_PI for $$\pi$$ in C/C++ must be replaced by pi for Expression objects.) In the present case, the output of u_code and f_code is

x + 2*x + 1
-10*x - 20*x - 10


After having defined the mesh, the function space, and the boundary, we define the boundary value u_D as

u_D = Expression(u_code, degree=1)


Similarly, we define the right-hand side function as

f = Expression(f_code, degree=1)


Name clash between FEniCS and program variables. In a program like the one above, strange errors may occur due to name clashes. If you define sym and q prior to doing from fenics import *, the latter statement will also import variables with the names sym and q, overwriting the objects you have previously defined! This may lead to strange errors. The safest solution is to do import fenics instead of from fenics import * and then prefix all FEniCS object names by fenics. The next best solution is to do from fenics import * first and then define your own variables that overwrite those imported from fenics. This is acceptable if we do not need sym and q from fenics.

### FEniCS implementation

A solver for the nonlinear Poisson equation is as easy to implement as a solver for the linear Poisson equation. All we need to do is to state the formula for $$F$$ and call solve(F == 0, u, bc) instead of solve(a == L, u, bc) as we did in the linear case. Here is a minimalistic code:

from fenics import *

def q(u):
return 1 + u**2

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)
u_D = Expression(u_code, degree=1)

def boundary(x, on_boundary):
return on_boundary

bc = DirichletBC(V, u_D, boundary)

u = Function(V)
v = TestFunction(V)
f = Expression(f_code, degree=1)

The major difference from a linear problem is that the unknown function u in the variational form in the nonlinear case must be defined as a Function, not as a TrialFunction. In some sense this is a simplification from the linear case where we must define u first as a TrialFunction and then as a Function.
The solve function takes the nonlinear equations, derives symbolically the Jacobian matrix, and runs a Newton method to compute the solution.
When we run the code, FEniCS reports on the progress of the Newton iterations. With $$2\cdot(8\times 8)$$ cells, we reach convergence in eight iterations with a tolerance of $$10^{-9}$$, and the error in the numerical solution is about $$10^{-16}$$. These results bring evidence for a correct implementation. Thinking in terms of finite differences on a uniform mesh, $$\mathsf{P}_1$$ elements mimic standard second-order differences, which compute the derivative of a linear or quadratic function exactly. Here, $$\nabla u$$ is a constant vector, but then multiplied by $$(1+u^2)$$, which is a second-order polynomial in $$x$$ and $$y$$, which the divergence "difference operator" should compute exactly. We can therefore, even with $$\mathcal{P}_1$$ elements, expect the manufactured $$u$$ to be reproduced by the numerical method. With a nonlinearity like $$1+u^4$$, this will not be the case, and we would need to verify convergence rates instead.