The standard way of compiling code with DOLFINx#

In this section, we will focus on the approach most users use to interact with UFL, FFCx and basix. Here we will start by creating the domain we want to solve a problem on. In this case, we will use a unit square

from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
import scipy

N = 10
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N)
tdim = mesh.topology.dim

Problem specification: Non-linear Poisson.#

Next, let’s consider the problem

\[\begin{split} \begin{align} -\nabla \cdot p(u) \nabla u &= f \quad \text{in } \Omega, \\ u &= g \quad \text{on } \partial \Omega \end{align} \end{split}\]

where \(p(u)=1+u^2\). We choose to use a first order Lagrange space for the unknown.

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

We do as in the previous section, and define a manufactured solution

def u_exact(module, x):
    return module.sin(2 * module.pi * x[1]) + x[0] * x[1] ** 2


x = ufl.SpatialCoordinate(mesh)
u_ex = u_exact(ufl, x)

As the problem is non-linear we need to prepare for solving this with a Newton-solver.

uh = dolfinx.fem.Function(V)

We define the residual

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


f = -ufl.div(p(u_ex) * ufl.grad(u_ex))
v = ufl.TestFunction(V)
F = ufl.inner(p(uh) * ufl.grad(uh), ufl.grad(v)) * ufl.dx - ufl.inner(f, v) * ufl.dx

gh = dolfinx.fem.Function(V)
gh.interpolate(lambda x: u_exact(np, x))

mesh.topology.create_entities(tdim - 1)
mesh.topology.create_connectivity(tdim - 1, tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, tdim - 1, boundary_facets)
bc = dolfinx.fem.dirichletbc(gh, boundary_dofs)
bcs = [bc]

We compute the Jacobian of the system using UFL.

J = ufl.derivative(F, uh)

Now that we have associated the forms with the discrete problem already, we use dolfinx.fem.form

residual = dolfinx.fem.form(F)
jacobian = dolfinx.fem.form(J)
A = dolfinx.fem.create_matrix(jacobian)
b = dolfinx.fem.create_vector(residual)

We are now ready to use these forms to solve the variational problem.

\[\begin{split} \begin{align} \frac{\partial F}{\partial u}[\delta u] &= F(u_k) \\ u_{k+1} = u_k - \delta u \end{align} \end{split}\]

As we want \(u_{k+1}=g\), but we do not know if \(u_k=g\), we need to take this into account when assembling the right hand side for the Jacobian equation.

We will go through how we apply this boundary condition using lifting

Lifting#

Let us split the degrees of freedom into two disjoint sets, \(u_d\), and \(u_{bc}\), and set up the corresponding linear system

(3)#\[\begin{split} \begin{align} \begin{pmatrix} A_{d,d} & A_{d, bc} \\ A_{bc,d} & A_{bc, bc} \end{pmatrix} \begin{pmatrix} u_d \\ u_{bc} \end{pmatrix} &= \begin{pmatrix} b_d \\ b_{bc} \end{pmatrix} \end{align} \end{split}\]

In the identity row approach, we set the rows corresponding to the Dirichlet conditions to the identity row and set the appropriate dofs on the right hand side to contain the Dirichlet values:

\[\begin{split} \begin{align} \begin{pmatrix} A_{d,d} & A_{d, bc} \\ 0 & I \end{pmatrix} \begin{pmatrix} u_d \\ u_{bc} \end{pmatrix} &= \begin{pmatrix} b_d \\ g \end{pmatrix} \end{align} \end{split}\]

where \(g\) is the vector satisfying the various Dirichlet conditions. We can now reduce this to a smaller system

\[ \begin{align} A_{d,d} u_d &= b_d - A_{d, bc}g \end{align} \]

which is symmetric if \(A_{d,d}\) is symmetric. However, we do not want to remove the degrees of freedom from the sparse matrix, as it makes the matrix less adaptable for varying boundary conditions, so we set up the system

\[\begin{split} \begin{align} \begin{pmatrix} A_{d,d} & 0 \\ 0 & I \end{pmatrix} \begin{pmatrix} u_d \\ u_{bc} \end{pmatrix} &= \begin{pmatrix} b_d - A_{d,bc}g