This tutorial is part of multiphenicsx.

Tutorial 03: weak imposition of Dirichlet BCs by a Lagrange multiplier (nonlinear problem)

In this tutorial we solve the problem

$$\begin{align*} &\min_{u} \int_\Omega \left\{ (1 + u^2)\ |\nabla u|^2 - u \right\} dx,\\ &\text{s.t. } u = g\text{ on }\Gamma = \partial \Omega \end{align*}$$ where $\Omega$ is a ball in 2D.

The optimality conditions result in the following nonlinear problem

$$\begin{align*} &\int_\Omega (1+u^2)\ \nabla u \cdot \nabla v \ dx + \int_\Omega u \ |\nabla u|^2 v \ dx = \int_\Omega v \ dx\\ &\text{s.t. } u = g\text{ on }\Gamma = \partial \Omega \end{align*}$$

We compare the following two cases: the corresponding weak formulation is

• strong imposition of Dirichlet BCs:

$$ \text{find } u \in V_g \text{ s.t. } \int_\Omega (1+u^2)\ \nabla u \cdot \nabla v \ dx + \int_\Omega u \ |\nabla u|^2 v \ dx = \int_\Omega v \ dx, \quad \forall v \in V_0\\ $$ where $$ V_g = \{v \in H^1(\Omega): v|_\Gamma = g\},\\ V_0 = \{v \in H^1(\Omega): v|_\Gamma = 0\}.\\ $$

• weak imposition of Dirichlet BCs: this requires an introduction of a multiplier $\lambda$ which is restricted to $\Gamma$, and solves

$$ \text{find } u, \lambda \in V \times M \text{ s.t. }\\ \begin{align*} &\int_\Omega (1+u^2)\ \nabla u \cdot \nabla v \ dx + \int_\Omega u \ |\nabla u|^2 v \ dx & &+ \int_\Gamma \lambda v \ dx & &= \int_\Omega v \ dx, & \forall v \in V,\\ &\int_\Gamma u \mu \ ds & & & &= \int_\Gamma g \mu \ ds, & \forall \mu \in M \end{align*} $$

where $$ V = H^1(\Omega),\\ M = L^{2}(\Gamma).\\ $$

This example is a prototypical case of problems containing subdomain/boundary restricted variables (the Lagrange multiplier, in this case).

import typing

import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex

import multiphenicsx.fem
import multiphenicsx.fem.petsc

Geometrical parameters

r = 3
mesh_size = 1. / 4.

Mesh

gmsh.initialize()
gmsh.model.add("mesh")

p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, mesh_size)
p1 = gmsh.model.geo.addPoint(0.0, +r, 0.0, mesh_size)
p2 = gmsh.model.geo.addPoint(0.0, -r, 0.0, mesh_size)
c0 = gmsh.model.geo.addCircleArc(p1, p0, p2)
c1 = gmsh.model.geo.addCircleArc(p2, p0, p1)
boundary = gmsh.model.geo.addCurveLoop([c0, c1])
domain = gmsh.model.geo.addPlaneSurface([boundary])

gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [c0, c1], 1)
gmsh.model.addPhysicalGroup(2, [boundary], 0)
gmsh.model.mesh.generate(2)

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Circle)
Info    : [ 60%] Meshing curve 2 (Circle)
Info    : Done meshing 1D (Wall 0.000258373s, CPU 0.000397s)
Info    : Meshing 2D...
Info    : Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.0122594s, CPU 0.008594s)
Info    : 589 nodes 1177 elements

mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()

# Create connectivities required by the rest of the code
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)

facets_Gamma = boundaries.indices[boundaries.values == 1]

viskex.dolfinx.plot_mesh(mesh)

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Weak imposition of Dirichlet BCs

# Define a function space
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
M = V.clone()

Create element

# Define restrictions
dofs_V = np.arange(0, V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts)
dofs_M_Gamma = dolfinx.fem.locate_dofs_topological(M, boundaries.dim, facets_Gamma)
restriction_V = multiphenicsx.fem.DofMapRestriction(V.dofmap, dofs_V)
restriction_M_Gamma = multiphenicsx.fem.DofMapRestriction(M.dofmap, dofs_M_Gamma)
restriction = [restriction_V, restriction_M_Gamma]

# Define trial and test functions, as well as solution
(du, dl) = (ufl.TrialFunction(V), ufl.TrialFunction(M))
(u, l) = (dolfinx.fem.Function(V), dolfinx.fem.Function(M))
(v, m) = (ufl.TestFunction(V), ufl.TestFunction(M))

# Define problem block forms
g = dolfinx.fem.Function(V)
g.interpolate(lambda x: np.sin(3 * x[0] + 1) * np.sin(3 * x[1] + 1))
F = [(ufl.inner((1 + u**2) * ufl.grad(u), ufl.grad(v)) * ufl.dx
      + ufl.inner(ufl.dot(ufl.grad(u), ufl.grad(u)) * u, v) * ufl.dx
      + ufl.inner(l, v) * ufl.ds - ufl.inner(1, v) * ufl.dx),
     ufl.inner(u, m) * ufl.ds - ufl.inner(g, m) * ufl.ds]
J = [[ufl.derivative(F[0], u, du), ufl.derivative(F[0], l, dl)],
     [ufl.derivative(F[1], u, du), ufl.derivative(F[1], l, dl)]]

# Class for interfacing with the SNES
class NonlinearLagrangeMultplierBlockProblem:
    """Define a nonlinear problem, interfacing with SNES."""

    def __init__(  # type: ignore[no-any-unimported]
        self, F: list[ufl.Form], J: list[list[ufl.Form]],
        solutions: tuple[dolfinx.fem.Function, dolfinx.fem.Function],
        bcs: list[dolfinx.fem.DirichletBC],
        P: typing.Optional[list[list[ufl.Form]]] = None
    ) -> None:
        self._F = dolfinx.fem.form(F)
        self._J = dolfinx.fem.form(J)
        self._obj_vec = multiphenicsx.fem.petsc.create_vector_block(self._F, restriction)
        self._solutions = solutions
        self._bcs = bcs
        self._P = P

    def create_snes_solution(self) -> petsc4py.PETSc.Vec:  # type: ignore[no-any-unimported]
        """
        Create a petsc4py.PETSc.Vec to be passed to petsc4py.PETSc.SNES.solve.

        The returned vector will be initialized with the initial guesses provided in `self._solutions`,
        properly stacked together and restricted in a single block vector.
        """
        x = multiphenicsx.fem.petsc.create_vector_block(self._F, restriction=restriction)
        with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(x, [V.dofmap, M.dofmap], restriction) as x_wrapper:
            for x_wrapper_local, sub_solution in zip(x_wrapper, self._solutions):
                with sub_solution.x.petsc_vec.localForm() as sub_solution_local:
                    x_wrapper_local[:] = sub_solution_local
        return x

    def update_solutions(self, x: petsc4py.PETSc.Vec) -> None:  # type: ignore[no-any-unimported]
        """Update `self._solutions` with data in `x`."""
        x.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
        with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(x, [V.dofmap, M.dofmap], restriction) as x_wrapper:
            for x_wrapper_local, sub_solution in zip(x_wrapper, self._solutions):
                with sub_solution.x.petsc_vec.localForm() as sub_solution_local:
                    sub_solution_local[:] = x_wrapper_local

    def obj(  # type: ignore[no-any-unimported]
            self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec
    ) -> np.float64:
        """Compute the norm of the residual."""
        self.F(snes, x, self._obj_vec)
        return self._obj_vec.norm()  # type: ignore[no-any-return]

    def F(  # type: ignore[no-any-unimported]
        self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, F_vec: petsc4py.PETSc.Vec
    ) -> None:
        """Assemble the residual."""
        self.update_solutions(x)
        with F_vec.localForm() as F_vec_local:
            F_vec_local.set(0.0)
        multiphenicsx.fem.petsc.assemble_vector_block(  # type: ignore[misc]
            F_vec, self._F, self._J, self._bcs, x0=x, alpha=-1.0,
            restriction=restriction, restriction_x0=restriction)

    def J(  # type: ignore[no-any-unimported]
        self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, J_mat: petsc4py.PETSc.Mat,
        P_mat: petsc4py.PETSc.Mat
    ) -> None:
        """Assemble the jacobian."""
        J_mat.zeroEntries()
        multiphenicsx.fem.petsc.assemble_matrix_block(
            J_mat, self._J, self._bcs, diagonal=1.0,  # type: ignore[arg-type]
            restriction=(restriction, restriction))
        J_mat.assemble()
        if self._P is not None:
            P_mat.zeroEntries()
            multiphenicsx.fem.petsc.assemble_matrix_block(
                P_mat, self._P, self._bcs, diagonal=1.0,  # type: ignore[arg-type]
                restriction=(restriction, restriction))
            P_mat.assemble()

# Create problem
problem = NonlinearLagrangeMultplierBlockProblem(F, J, (u, l), [])
F_vec = multiphenicsx.fem.petsc.create_vector_block(problem._F, restriction=restriction)
J_mat = multiphenicsx.fem.petsc.create_matrix_block(problem._J, restriction=(restriction, restriction))

# Solve
snes = petsc4py.PETSc.SNES().create(mesh.comm)
snes.setTolerances(max_it=20)
snes.getKSP().setType("preonly")
snes.getKSP().getPC().setType("lu")
snes.getKSP().getPC().setFactorSolverType("mumps")
snes.setObjective(problem.obj)
snes.setFunction(problem.F, F_vec)
snes.setJacobian(problem.J, J=J_mat, P=None)
snes.setMonitor(lambda _, it, residual: print(it, residual))
solution = problem.create_snes_solution()
snes.solve(None, solution)
problem.update_solutions(solution)  # TODO can this be safely removed?
solution.destroy()
snes.destroy()

0 1.0636696615937857
1 0.5297763037133386
2 0.326660348283296
3 0.18048341314104063
4 9.069054492574919e-05
5 6.95366681680224e-10

<petsc4py.PETSc.SNES at 0x7f349bf1fa60>

viskex.dolfinx.plot_scalar_field(u, "u")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_scalar_field(l, "l")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Strong imposition of Dirichlet BCs

# Class for interfacing with the SNES
class NonlinearStrongBoundaryConditionsProblem:
    """Define a nonlinear problem, interfacing with SNES."""

    def __init__(  # type: ignore[no-any-unimported]
        self, F: ufl.Form, J: ufl.Form, solution: dolfinx.fem.Function,
        bcs: list[dolfinx.fem.DirichletBC], P: typing.Optional[ufl.Form] = None
    ) -> None:
        self._F = dolfinx.fem.form(F)
        self._J = dolfinx.fem.form(J)
        self._obj_vec = dolfinx.fem.petsc.create_vector(self._F)
        self._solution = solution
        self._bcs = bcs
        self._P = P

    def create_snes_solution(self) -> petsc4py.PETSc.Vec:  # type: ignore[no-any-unimported]
        """
        Create a petsc4py.PETSc.Vec to be passed to petsc4py.PETSc.SNES.solve.

        The returned vector will be initialized with the initial guess provided in `self._solution`.
        """
        x = self._solution.x.petsc_vec.copy()
        with x.localForm() as _x, self._solution.x.petsc_vec.localForm() as _solution:
            _x[:] = _solution
        return x

    def update_solution(self, x: petsc4py.PETSc.Vec) -> None:  # type: ignore[no-any-unimported]
        """Update `self._solution` with data in `x`."""
        x.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
        with x.localForm() as _x, self._solution.x.petsc_vec.localForm() as _solution:
            _solution[:] = _x

    def obj(  # type: ignore[no-any-unimported]
        self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec
    ) -> np.float64:
        """Compute the norm of the residual."""
        self.F(snes, x, self._obj_vec)
        return self._obj_vec.norm()  # type: ignore[no-any-return]

    def F(  # type: ignore[no-any-unimported]
        self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, F_vec: petsc4py.PETSc.Vec
    ) -> None:
        """Assemble the residual."""
        self.update_solution(x)
        with F_vec.localForm() as F_vec_local:
            F_vec_local.set(0.0)
        dolfinx.fem.petsc.assemble_vector(F_vec, self._F)
        dolfinx.fem.petsc.apply_lifting(F_vec, [self._J], [self._bcs], x0=[x], alpha=-1.0)
        F_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
        dolfinx.fem.petsc.set_bc(F_vec, self._bcs, x, -1.0)

    def J(  # type: ignore[no-any-unimported]
        self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, J_mat: petsc4py.PETSc.Mat,
        P_mat: petsc4py.PETSc.Mat
    ) -> None:
        """Assemble the jacobian."""
        J_mat.zeroEntries()
        dolfinx.fem.petsc.assemble_matrix(  # type: ignore[misc]
            J_mat, self._J, self._bcs, diagonal=1.0)  # type: ignore[arg-type]
        J_mat.assemble()
        if self._P is not None:
            P_mat.zeroEntries()
            dolfinx.fem.petsc.assemble_matrix(  # type: ignore[misc]
                P_mat, self._P, self._bcs, diagonal=1.0)  # type: ignore[arg-type]
            P_mat.assemble()

# Define problem block forms
u_ex = dolfinx.fem.Function(V)
F_ex = ufl.replace(F[0], {u: u_ex, l: 0})
J_ex = ufl.derivative(F_ex, u_ex, du)

# Define Dirichlet BC object on Gamma
dofs_V_Gamma = dolfinx.fem.locate_dofs_topological(V, boundaries.dim, facets_Gamma)
bc_ex = [dolfinx.fem.dirichletbc(g, dofs_V_Gamma)]

# Create problem
problem_ex = NonlinearStrongBoundaryConditionsProblem(F_ex, J_ex, u_ex, bc_ex)
F_ex_vec = dolfinx.fem.petsc.create_vector(problem_ex._F)
J_ex_mat = dolfinx.fem.petsc.create_matrix(problem_ex._J)

# Solve
snes = petsc4py.PETSc.SNES().create(mesh.comm)
snes.setTolerances(max_it=20)
snes.getKSP().setType("preonly")
snes.getKSP().getPC().setType("lu")
snes.getKSP().getPC().setFactorSolverType("mumps")
snes.setObjective(problem_ex.obj)
snes.setFunction(problem_ex.F, F_ex_vec)
snes.setJacobian(problem_ex.J, J=J_ex_mat, P=None)
snes.setMonitor(lambda _, it, residual: print(it, residual))
solution_ex = problem_ex.create_snes_solution()
snes.solve(None, solution_ex)
problem_ex.update_solution(solution_ex)  # TODO can this be safely removed?
solution_ex.destroy()
snes.destroy()

0 11.903671385189543
1 1.5766579840995987
2 0.3673953356946283
3 0.05286757521106234
4 0.0017909741046908975
5 2.2418618294007724e-06
6 3.348850757573569e-12

<petsc4py.PETSc.SNES at 0x7f3440307f10>

viskex.dolfinx.plot_scalar_field(u_ex, "u")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Comparison and error computation

u_ex_norm = np.sqrt(mesh.comm.allreduce(
    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u_ex), ufl.grad(u_ex)) * ufl.dx)),
    op=mpi4py.MPI.SUM))
err_norm = np.sqrt(mesh.comm.allreduce(
    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u_ex - u), ufl.grad(u_ex - u)) * ufl.dx)),
    op=mpi4py.MPI.SUM))
print("Relative error is equal to", err_norm / u_ex_norm)
assert np.isclose(err_norm / u_ex_norm, 0., atol=1.e-9)

Relative error is equal to 1.5375129401353526e-10