Real function spaces

Real function spaces#

Author: Jørgen S. Dokken License: MIT In this example we will show how to use the “real” function space to solve a singular Poisson problem. The problem at hand is: Find

u \in H^{1} (Ω)

such that

(1)#

\begin{aligned} - Δ u & = f in Ω, \\ \frac{\partial u}{\partial n} & = g on \partial Ω, \\ \int_{Ω} u & = h . \end{aligned}

Lagrange multiplier#

We start by considering the equivalent optimization problem: Find $u \in H^{1} (Ω)$ such that

(2)#

\begin{array}{r} min_{u \in H^{1} (Ω)} J (u) = min_{u \in H^{1} (Ω)} \frac{1}{2} \int_{Ω} | \nabla u \cdot \nabla u | d x - \int_{Ω} f u d x - \int_{\partial Ω} g u d s, \end{array}

such that

(3)#

\begin{array}{r} \int_{Ω} u = h . \end{array}

We introduce a Lagrange multiplier $λ$ to enforce the constraint:

(4)#

\begin{array}{r} min_{u \in H^{1} (Ω), λ \in R} L (u, λ) = min_{u \in H^{1} (Ω), λ \in R} J (u) + λ (\int_{Ω} u d x - h) . \end{array}

We then compute the optimality conditions for the problem above

(5)#

\begin{aligned} \frac{\partial L}{\partial u} [δ u] & = \int_{Ω} \nabla u \cdot \nabla δ u d x + λ \int δ u d x - \int_{Ω} f δ u d x - \int_{\partial Ω} g δ u d s = 0, \\ \frac{\partial L}{\partial λ} [δ λ] & = δ λ (\int_{Ω} u d x - h) = 0. \end{aligned}

We write the weak formulation:

\begin{array}{r} \begin{matrix} \int_{Ω} \nabla u \cdot \nabla δ u d x + \int_{Ω} λ δ u d x = \int_{Ω} f δ u d x + \int_{\partial Ω} g v d s \\ \int_{Ω} u δ λ d x = h δ λ . \end{matrix} \end{array}

where we have moved $δ λ$ into the integral as it is a spatial constant.

Implementation#

We start by creating the domain and derive the source terms $f$ , $g$ and $h$ from our manufactured solution For this example we will use the following exact solution

(6)#

\begin{array}{r} u_{e x a c t} (x, y) = 0.3 y^{2} + \sin (2 π x) . \end{array}

from packaging.version import Version
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.la.petsc import scatter_local_vectors, get_local_vectors
import dolfinx.fem.petsc

import numpy as np
from scifem import create_real_functionspace, assemble_scalar
import ufl

M = 20
mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, M, M, dolfinx.mesh.CellType.triangle, dtype=np.float64
)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

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

x = ufl.SpatialCoordinate(mesh)
n = ufl.FacetNormal(mesh)
g = ufl.dot(ufl.grad(u_exact(x)), n)
f = -ufl.div(ufl.grad(u_exact(x)))
h = assemble_scalar(u_exact(x) * ufl.dx)

We then create the Lagrange multiplier space

R = create_real_functionspace(mesh)

Next, we can create a mixed-function space for our problem

if dolfinx.__version__ == "0.8.0":
    u = ufl.TrialFunction(V)
    lmbda = ufl.TrialFunction(R)
    du = ufl.TestFunction(V)
    dl = ufl.TestFunction(R)
elif Version(dolfinx.__version__) >= Version("0.9.0.0"):
    W = ufl.MixedFunctionSpace(V, R)
    u, lmbda = ufl.TrialFunctions(W)
    du, dl = ufl.TestFunctions(W)
else:
    raise RuntimeError("Unsupported version of dolfinx")

We can now define the variational problem

zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0))

a00 = ufl.inner(ufl.grad(u), ufl.grad(du)) * ufl.dx
a01 = ufl.inner(lmbda, du) * ufl.dx
a10 = ufl.inner(u, dl) * ufl.dx
L0 = ufl.inner(f, du) * ufl.dx + ufl.inner(g, du) * ufl.ds
L1 = ufl.inner(zero, dl) * ufl.dx

a = dolfinx.fem.form([[a00, a01], [a10, None]])
L = dolfinx.fem.form([L0, L1])

Note that we have defined the variational form in a block form, and that we have not included $h$ in the variational form. We will enforce this once we have assembled the right hand side vector.

We can now assemble the matrix and vector

A = dolfinx.fem.petsc.assemble_matrix_block(a)
A.assemble()
b = dolfinx.fem.petsc.assemble_vector_block(L, a, bcs=[])

Next, we modify the second part of the block to contain h We start by enforcing the multiplier constraint $h$ by modifying the right hand side vector

if dolfinx.__version__ == "0.8.0":
    maps = [(V.dofmap.index_map, V.dofmap.index_map_bs), (R.dofmap.index_map, R.dofmap.index_map_bs)]
elif Version(dolfinx.__version__) >= Version("0.9.0.0"):
    maps = [(Wi.dofmap.index_map, Wi.dofmap.index_map_bs) for Wi in W.ufl_sub_spaces()]

b_local = get_local_vectors(b, maps)
b_local[1][:] = h
scatter_local_vectors(
        b,
        b_local,
        maps,
    )
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

We can now solve the linear system

ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")

xh = dolfinx.fem.petsc.create_vector_block(L)
ksp.solve(b, xh)
xh.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

Finally, we extract the solution u from the blocked system and compute the error

uh = dolfinx.fem.Function(V, name="u")
x_local = get_local_vectors(xh, maps)
uh.x.array[: len(x_local[0])] = x_local[0]
uh.x.scatter_forward()

diff = uh - u_exact(x)
error = dolfinx.fem.form(ufl.inner(diff, diff) * ufl.dx)

print(f"L2 error: {np.sqrt(assemble_scalar(error)):.2e}")

L2 error: 6.73e-03

We can now plot the solution

vtk_mesh = dolfinx.plot.vtk_mesh(V)

import pyvista

pyvista.start_xvfb()
grid = pyvista.UnstructuredGrid(*vtk_mesh)
grid.point_data["u"] = uh.x.array.real

warped = grid.warp_by_scalar("u", factor=1)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, style="wireframe")
plotter.add_mesh(warped)
if not pyvista.OFF_SCREEN:
    plotter.show()

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

Real function spaces

Contents

Real function spaces#

Lagrange multiplier#

Implementation#