In this tutorial we compute the inf-sup constant of the saddle point problem resulting from the discretization of the following Laplace problem $$\begin{cases} -\Delta u = f, & \text{in } \Omega,\\ u = 0, & \text{on } \Gamma = \partial\Omega, \end{cases}$$
where $\Omega$ is a ball in 2D, and for which the non-homogeneous Dirichlet boundary conditions are imposed by a Lagrange multiplier.
The resulting eigenvalue problem is \begin{align*} &\text{find } \eta, u, \lambda \in \mathbb{R} \times V \times M \text{ s.t. }\\ &\begin{cases} \int_\Omega \nabla u \cdot \nabla v + \int_\Gamma \lambda v = 0, & \forall v \in V,\\ \int_\Gamma u \mu = \eta \int_\Gamma \lambda \mu, & \forall \mu \in M \end{cases} \end{align*} where $$ V = H^1(\Omega),\\ M = L^{2}(\Gamma).\\ $$
import dolfinx.fem
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import slepc4py.SLEPc
import ufl
import viskex
import multiphenicsx.fem
import multiphenicsx.fem.petsc
r = 3
mesh_size = 1. / 4.
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.000297484s, CPU 0.000469s) Info : Meshing 2D... Info : Meshing surface 1 (Plane, Frontal-Delaunay) Info : Done meshing 2D (Wall 0.0125968s, CPU 0.01273s) 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
# Define a function space
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
# Define restrictions.
dofs_V = np.arange(0, V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts)
dofs_V_Gamma = dolfinx.fem.locate_dofs_topological(V, boundaries.dim, facets_Gamma)
restriction_V = multiphenicsx.fem.DofMapRestriction(V.dofmap, dofs_V)
restriction_V_Gamma = multiphenicsx.fem.DofMapRestriction(V.dofmap, dofs_V_Gamma)
restriction = [restriction_V, restriction_V_Gamma]
# Define trial and test functions
(u, l) = (ufl.TrialFunction(V), ufl.TrialFunction(V))
(v, m) = (ufl.TestFunction(V), ufl.TestFunction(V))
# Define problem block forms
a = [[ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx, - ufl.inner(l, v) * ufl.ds],
[- ufl.inner(u, m) * ufl.ds, None]]
b = [[None, None],
[None, - ufl.inner(l, m) * ufl.ds]]
b[0][0] = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)) * ufl.inner(u, v) * ufl.dx
# Assemble lhs and rhs matrices
A = multiphenicsx.fem.petsc.assemble_matrix_block(
dolfinx.fem.form(a), bcs=[], restriction=(restriction, restriction))
A.assemble()
B = multiphenicsx.fem.petsc.assemble_matrix_block(
dolfinx.fem.form(b), bcs=[], restriction=(restriction, restriction))
B.assemble()
# Solve
eps = slepc4py.SLEPc.EPS().create(mesh.comm)
eps.setOperators(A, B)
eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GNHEP)
eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)
eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.TARGET_REAL)
eps.setTarget(1.e-5)
eps.getST().setType(slepc4py.SLEPc.ST.Type.SINVERT)
eps.getST().getKSP().setType("preonly")
eps.getST().getKSP().getPC().setType("lu")
eps.getST().getKSP().getPC().setFactorSolverType("mumps")
eps.solve()
assert eps.getConverged() >= 1
# Extract leading eigenvalue
eigv = eps.getEigenvalue(0)
r, i = eigv.real, eigv.imag
assert abs(i) < 1.e-10
assert r > 0., "r = " + str(r) + " is not positive"
print("Inf-sup constant: ", np.sqrt(r))
assert np.isclose(np.sqrt(r), 0.125429, rtol=np.inf)
Inf-sup constant: 0.12542987783418216
eps.destroy()
<slepc4py.SLEPc.EPS at 0x7f156026c3b0>