This tutorial is part of multiphenicsx.

Tutorial 01, case 2: 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:

• strong imposition of Dirichlet BCs: the corresponding weak formulation is $$ \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 dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
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.00028805s, CPU 0.00061s)
Info    : Meshing 2D...
Info    : Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.0120053s, CPU 0.012289s)
Info    : 589 nodes 1177 elements

mesh, subdomains, boundaries, *other_tags = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize
assert subdomains is not None
assert boundaries is not None

# 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)

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

Weak imposition of Dirichlet BCs

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

# 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 test functions, as well as solutions
(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]

# Solve
petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "ksp_error_if_not_converged": True,
    "snes_monitor": None,
    "snes_error_if_not_converged": True
}
problem = multiphenicsx.fem.petsc.NonlinearProblem(
    F, (u, l),
    petsc_options_prefix="tutorial_1_lagrange_multipliers_nonlinear_", petsc_options=petsc_options,
    kind="mpi", restriction=restriction
)
problem.solve()
del problem

  0 SNES Function norm 1.063669661594e+00
  1 SNES Function norm 8.244952454609e-01
  2 SNES Function norm 5.669031385461e-01
  3 SNES Function norm 3.552462742210e-01
  4 SNES Function norm 2.757370255700e-01
  5 SNES Function norm 4.004704369053e-04
  6 SNES Function norm 3.997907337463e-08
  7 SNES Function norm 2.701705814475e-14

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

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

Strong imposition of Dirichlet BCs

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

# 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)]

# Solve
problem_ex = dolfinx.fem.petsc.NonlinearProblem(
    F_ex, u_ex, bcs=bc_ex,
    petsc_options_prefix="tutorial_1_lagrange_multipliers_nonlinear_ex_", petsc_options=petsc_options
)
problem_ex.solve()
del problem_ex

  0 SNES Function norm 1.190367138519e+01
  1 SNES Function norm 1.576657984100e+00
  2 SNES Function norm 3.673953356946e-01
  3 SNES Function norm 5.286757521106e-02
  4 SNES Function norm 1.790974104689e-03
  5 SNES Function norm 2.241861828302e-06
  6 SNES Function norm 3.348422405820e-12

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

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 9.157442065436912e-13