I am solving a Poisson problem using hybridized discontinuous Galerkin and want to perform static condensation of the local solves DoFs.
The main problem is
-\nabla \cdot \nabla u =0
with u_{\Gamma_{left}}=0, u_{\Gamma_{right}}=1, and \nabla u \cdot \pmb{n}|_{\partial \Omega \setminus (\Gamma_{left}\bigcup\Gamma_{right}})=0.
Where \pmb{x} are the local solutions and \pmb{y} is the global solution, I’m solving:
\begin{bmatrix}
\pmb{A} & \pmb{B}\\
\pmb{C} & \pmb{D}
\end{bmatrix}
\begin{pmatrix}
\pmb{x}\\
\pmb{y}
\end{pmatrix}
=
\begin{pmatrix}
\pmb{0}\\
\pmb{f_y}
\end{pmatrix}
and want to have
\pmb{K_y} = \pmb{D} - \pmb{C}\pmb{A}^{-1}\pmb{B}
\pmb{b_y} = \pmb{f_y}
so that I can solve \pmb{K_y}\pmb{y}=\pmb{b_y}.
In the below MWE, I’m unclear what goes into dolfinx.fem.form_cpp_class especially since I’m working with submeshes. As a result, I’m getting an error that
Traceback (most recent call last):
File "/home/lmolel/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/numba/np/linalg.py", line 841, in _check_finite_matrix
raise np.linalg.LinAlgError(
numpy.linalg.LinAlgError: Array must not contain infs or NaNs.
inf
Here’s the MWE
#!/usr/bin/env python3
import argparse
import json
import os
import shutil
import sys
import time
sys.path.append("../")
import gmsh
import matplotlib.pyplot as plt
import matspy
import numpy as np
import scifem
import scipy
import ufl
from dolfinx import cpp, fem, io, log, mesh
from dolfinx import default_real_type, default_scalar_type
from dolfinx.graph import partitioner_scotch
from dolfinx.jit import ffcx_jit
import dolfinx.fem.petsc as petsc
import numba.core.typing.cffi_utils as cffi_support
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import div, dot, grad, inner
import cffi
import numba
from ffcx.codegeneration.utils import empty_void_pointer
from ffcx.codegeneration.utils import numba_ufcx_kernel_signature as ufcx_signature
Print = PETSc.Sys.Print
dtype = default_scalar_type
rtype = default_real_type
_type_to_offset_index = {fem.IntegralType.cell: 0, fem.IntegralType.exterior_facet: 1}
ffi = cffi.FFI()
left = 1
right = 2
def compute_cell_boundary_facets(msh):
"""Compute the integration entities for integrals around the
boundaries of all cells in msh.
Parameters:
msh: The mesh.
Returns:
Facets to integrate over, identified by ``(cell, local facet
index)`` pairs.
"""
tdim = msh.topology.dim
fdim = tdim - 1
n_f = cpp.mesh.cell_num_entities(msh.topology.cell_type, fdim)
n_c = msh.topology.index_map(tdim).size_local
return np.vstack((np.repeat(np.arange(n_c), n_f), np.tile(np.arange(n_f), n_c))).T.flatten()
def create_mesh(LX, LY):
gmsh.initialize()
gmsh.model.add('ssb')
# gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.01)
points = [
(0, 0, 0),
(LX, 0, 0),
(LX, LY, 0),
(0, LY, 0),
]
gpoints = [gmsh.model.occ.addPoint(*p) for p in points]
lines = [gmsh.model.occ.addLine(gpoints[i], gpoints[i+1]) for i in range(-1, len(gpoints)-1)]
curve_loop = gmsh.model.occ.addCurveLoop(lines)
surf = gmsh.model.occ.addPlaneSurface([curve_loop])
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(1, [1], left, "Left")
gmsh.model.addPhysicalGroup(1, [4], right, "Right")
gmsh.model.addPhysicalGroup(2, [surf], 1, "domain")
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
gmsh.finalize()
if __name__ == '__main__':
comm = MPI.COMM_WORLD
if comm.rank == 0:
create_mesh(1, 1)
partitioner = mesh.create_cell_partitioner(partitioner_scotch(), mesh.GhostMode.shared_facet)
domain, ct, ft = io.gmsh.read_from_msh("mesh.msh", comm, partitioner=partitioner)[:3]
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_entities(fdim)
domain.topology.create_connectivity(tdim, fdim)
domain.topology.create_connectivity(fdim, tdim)
domain.topology.create_connectivity(fdim, fdim)
ct_imap = domain.topology.index_map(tdim)
num_entities_local = ct_imap.size_local + ct_imap.num_ghosts
# tag internal facets as 0
ft_imap = domain.topology.index_map(fdim)
num_facets_local = ft_imap.size_local + ft_imap.num_ghosts
# facets
num_facets_local = (
domain.topology.index_map(fdim).size_local + domain.topology.index_map(fdim).num_ghosts
)
facets = np.arange(num_facets_local, dtype=np.int32)
values = np.full_like(facets, 0, dtype=np.int32)
values[ft.find(left)] = left
values[ft.find(right)] = right
ft = mesh.meshtags(domain, fdim, facets, values)
domain.topology.create_connectivity(fdim, tdim)
f_to_c = domain.topology.connectivity(fdim, tdim)
x = ufl.SpatialCoordinate(domain)
dx = ufl.Measure('dx', domain=domain, subdomain_data=ct)
# facet mesh
facet_mesh, facet_mesh_emap = mesh.create_submesh(domain, fdim, ft.indices)[:2]
entity_maps = [facet_mesh_emap]
k = 1
V = fem.functionspace(domain, ("DG", k))
Vbar = fem.functionspace(facet_mesh, ("DG", k))
W = ufl.MixedFunctionSpace(V, Vbar)
u, ubar = ufl.TrialFunctions(W)
v, vbar = ufl.TestFunctions(W)
h = ufl.CellDiameter(domain)
n = ufl.FacetNormal(domain)
gamma = 16 * k ** 2 / h
cell_boundaries = 1 # A tag
cell_boundary_facets = compute_cell_boundary_facets(domain)
# Create the measure
ds_c = ufl.Measure("ds", subdomain_data=[(cell_boundaries, cell_boundary_facets)], domain=domain)
f_to_c = domain.topology.connectivity(fdim, tdim)
c_to_f = domain.topology.connectivity(tdim, fdim)
a = (
ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
- ufl.inner((u - ubar), ufl.dot(ufl.grad(v), n)) * ds_c(cell_boundaries)
- ufl.inner(ufl.dot(ufl.grad(u), n), (v - vbar)) * ds_c(cell_boundaries)
+ gamma * ufl.inner((u - ubar), v - vbar) * ds_c(cell_boundaries)
)
a_blocks = ufl.extract_blocks(a)
a00 = a_blocks[0][0]
a01 = a_blocks[0][1]
a10 = a_blocks[1][0]
a11 = a_blocks[1][1]
facet_mesh.topology.create_connectivity(fdim, fdim)
left_boundary_facets = facet_mesh_emap.sub_topology_to_topology(
ft.find(left), inverse=True
)
right_boundary_facets = facet_mesh_emap.sub_topology_to_topology(
ft.find(right), inverse=True
)
# Get the dofs and apply the boundary condition
facet_mesh.topology.create_connectivity(fdim, fdim)
left_dofs = fem.locate_dofs_topological(Vbar, fdim, left_boundary_facets)
right_dofs = fem.locate_dofs_topological(Vbar, fdim, right_boundary_facets)
left_bc = fem.dirichletbc(dtype(0.0), left_dofs, Vbar)
right_bc = fem.dirichletbc(dtype(1.0), right_dofs, Vbar)
bcs = [right_bc]
ufcx00, _, _ = ffcx_jit(comm, a00, form_compiler_options={"scalar_type": dtype})
kernel00 = getattr(ufcx00.form_integrals[0], f"tabulate_tensor_{np.dtype(dtype).name}")
ufcx01, _, _ = ffcx_jit(comm, a01, form_compiler_options={"scalar_type": dtype})
kernel01 = getattr(ufcx01.form_integrals[0], f"tabulate_tensor_{np.dtype(dtype).name}")
ufcx10, _, _ = ffcx_jit(comm, a10, form_compiler_options={"scalar_type": dtype})
kernel10 = getattr(ufcx10.form_integrals[0], f"tabulate_tensor_{np.dtype(dtype).name}")
ufcx11, _, _ = ffcx_jit(comm, a11, form_compiler_options={"scalar_type": dtype})
kernel11 = getattr(ufcx11.form_integrals[0], f"tabulate_tensor_{np.dtype(dtype).name}")
Vsize = V.element.space_dimension
Vbarsize = Vbar.element.space_dimension
@numba.cfunc(ufcx_signature(dtype, rtype), nopython=True)
def tabulate_A(A_, w_, c_, coords_, entity_local_index, permutation=ffi.NULL, custom_data=None):
"""Element kernel that applies static condensation."""
# Prepare target condensed local element tensor
A = numba.carray(A_, (Vbarsize, Vbarsize), dtype=dtype)
# Tabulate all sub blocks locally
A00 = np.zeros((Vsize, Vsize), dtype=dtype)
kernel00(
ffi.from_buffer(A00),
w_,
c_,
coords_,
entity_local_index,
permutation,
empty_void_pointer(),
)
A01 = np.zeros((Vsize, Vbarsize), dtype=dtype)
kernel01(
ffi.from_buffer(A01),
w_,
c_,
coords_,
entity_local_index,
permutation,
empty_void_pointer(),
)
A10 = np.zeros((Vbarsize, Vsize), dtype=dtype)
kernel10(
ffi.from_buffer(A10),
w_,
c_,
coords_,
entity_local_index,
permutation,
empty_void_pointer(),
)
A11 = np.zeros((Vbarsize, Vbarsize), dtype=dtype)
kernel11(
ffi.from_buffer(A11),
w_,
c_,
coords_,
entity_local_index,
permutation,
empty_void_pointer(),
)
# A = A11 - A10 * A00^{-1} * A01
A[:, :] = A11 - A10 @ np.linalg.solve(A00, A01)
formtype = fem.form_cpp_class(dtype)
facets = np.arange(domain.topology.index_map(fdim).size_local)
cells = np.arange(domain.topology.index_map(domain.topology.dim).size_local)
integrals = {fem.IntegralType.exterior_facet: [(0, tabulate_A.address, facets, np.array([], dtype=np.int8))]}
a_cond = fem.Form(
formtype([Vbar._cpp_object, Vbar._cpp_object], integrals, [], [], False, [], mesh=facet_mesh._cpp_object)
)
A_cond = fem.petsc.assemble_matrix(a_cond, bcs=bcs)
A_cond.assemble()
Print(A_cond)
b = fem.petsc.create_vector([Vbar])
fem.petsc.apply_lifting(b, [a_cond], bcs=[[right_bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
right_bc.set(b)
uc = fem.Function(Vbar, name="from_condensation")
solver = PETSc.KSP().create(A_cond.getComm())
solver.setType("gmres")
solver.getPC().setType("lu")
solver.setOperators(A_cond)
solver.solve(b, uc.x.petsc_vec)
Print(uc.x.petsc_vec.norm(0))
solver.destroy()