festim-dev
Penalty-approach for Sievert/Henry interfaces.
The current implementation of the Sievert/Henry interface condition with quadratic penalty assumes that we modify the boundary conditions, rather than adding a quadratic penalty. Following is an MWE, using DOLFINx nightly, that illustrates why we should use quadratic penalty as apposed to the modified boundary condition.
The manufactured solution of this problem is: Given a unit square, divided into two materials $\Omega_A=[0,d]\times[0,1]$ and $\Omega_B=[d, 1]x[0,1]$ with
$$ \begin{align*} -\nabla \cdot (D_A\nabla u_A)&=f_A&&\text{in } \Omega_A,\ -\nabla \cdot (D_B\nabla u_B)&=f_B&&\text{in } \Omega_B,\ u_A &= K_A\left(\frac{u_B}{K_B}\right)^n&&\text{on } \Gamma,\ D_A\nabla u_A \cdot n_A &= D_B\nabla u_B\cdot n_A && \text{on }\Gamma \end{align*} $$
where $n=\frac{1}{2}$ or $1$, where 0.5 correponds to $\Omega_A$ is Sievert, $\Omega_B$ Henry.
$$ \begin{align*} u_A&=x\ u_B&=\frac{n}{d^{\frac{1}{n}-1}}x^{\frac{1}{x^n}}\ K_A&= 1\ K_B&= n d^{\frac{n-1}{n}}\ D_A&=D_B=1 \end{align*} $$
from mpi4py import MPI
import dolfinx.fem.petsc
import ufl
import numpy as np
import numpy.typing as npt
from enum import Enum
import argparse
class PenaltyMode(Enum):
BC = 1
QUADRATIC = 2
def tr_a(v):
return v("+")
def tr_b(v):
return v("-")
def left_dofs(x):
return np.isclose(x[0], 0.0)
def right_dofs(x):
return np.isclose(x[0], 1.0)
def compute_interface_data(
cell_tags: dolfinx.mesh.MeshTags, facet_indices: npt.NDArray[np.int32]
) -> npt.NDArray[np.int32]:
"""
Compute interior facet integrals that are consistently ordered according to the `cell_tags`,
such that the data `(cell0, facet_idx0, cell1, facet_idx1)` is ordered such that
`cell_tags[cell0]`<`cell_tags[cell1]`, i.e the cell with the lowest cell marker is considered the
"+" restriction".
Args:
cell_tags: MeshTags that must contain an integer marker for all cells adjacent to the `facet_indices`
facet_indices: List of facets (local index) that are on the interface.
Returns:
The integration data.
"""
cell_tags.topology.create_connectivity(cell_tags.dim - 1, cell_tags.dim)
idata = dolfinx.cpp.fem.compute_integration_domains(
dolfinx.fem.IntegralType.interior_facet,
cell_tags.topology,
facet_indices,
)
ordered_idata = idata.reshape(-1, 4).copy()
switch = (
cell_tags.values[ordered_idata[:, 0]] > cell_tags.values[ordered_idata[:, 2]]
)
if True in switch:
ordered_idata[switch, :] = ordered_idata[switch][:, [2, 3, 0, 1]]
return ordered_idata
def solve_penalty_problem(
penalty_mode: PenaltyMode,
d: float = 0.5,
Nx: int = 10,
Ny: int = 10,
alpha_c: float = 1000.0,
):
def split(x, tol=1e2 * np.finfo(dolfinx.default_scalar_type).eps):
"""Locator for right subdomain"""
return x[0] <= d + tol
mesh = dolfinx.mesh.create_unit_square(
MPI.COMM_WORLD, Nx, Ny, ghost_mode=dolfinx.mesh.GhostMode.shared_facet
)
tdim = mesh.topology.dim
cell_map = mesh.topology.index_map(tdim)
num_cells = cell_map.size_local + cell_map.num_ghosts
values = np.full(num_cells, 2, dtype=np.int32)
values[dolfinx.mesh.locate_entities(mesh, tdim, split)] = 1
ct = dolfinx.mesh.meshtags(mesh, tdim, np.arange(num_cells, dtype=np.int32), values)
omega_A, map_A, v_map_A, n_map_A = dolfinx.mesh.create_submesh(
mesh, tdim, ct.find(1)
)
omega_B, map_B, v_map_B, n_map_B = dolfinx.mesh.create_submesh(
mesh, tdim, ct.find(2)
)
V_A = dolfinx.fem.functionspace(omega_A, ("Lagrange", 1))
V_B = dolfinx.fem.functionspace(omega_B, ("Lagrange", 1))
W = ufl.MixedFunctionSpace(*(V_A, V_B))
u_a = dolfinx.fem.Function(V_A)
u_b = dolfinx.fem.Function(V_B)
v_a, v_b = ufl.TestFunctions(W)
facets_A = dolfinx.mesh.compute_incident_entities(
mesh.topology, ct.find(1), tdim, tdim - 1
)
facets_B = dolfinx.mesh.compute_incident_entities(
mesh.topology, ct.find(2), tdim, tdim - 1
)
gamma = np.intersect1d(facets_A, facets_B)
dfx = dolfinx.default_scalar_type
n_sorption = dolfinx.fem.Constant(mesh, dfx(1.0 / 2))
bc_A = dolfinx.fem.dirichletbc(
dolfinx.fem.Constant(omega_A, dfx(0.0)),
dolfinx.fem.locate_dofs_geometrical(V_A, left_dofs),
V_A,
)
bc_B = dolfinx.fem.dirichletbc(
dolfinx.fem.Constant(omega_B, dfx(n_sorption / d ** (1 / n_sorption - 1))),
dolfinx.fem.locate_dofs_geometrical(V_B, right_dofs),
V_B,
)
bcs = [bc_A, bc_B]
idata = compute_interface_data(ct, gamma)
dGamma = ufl.Measure(
"dS", domain=mesh, subdomain_data=[(1, idata.flatten())], subdomain_id=1
)
dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)
dOmegaA = dx(1)
dOmegaB = dx(2)
x = ufl.SpatialCoordinate(mesh)
f_A = -ufl.div(ufl.grad(x[0]))
f_B = -ufl.div(
ufl.grad(n_sorption / (d ** (1 / n_sorption - 1)) * x[0] ** (1 / n_sorption))
)
F = ufl.inner(ufl.grad(u_a), ufl.grad(v_a)) * dOmegaA
F += ufl.inner(ufl.grad(u_b), ufl.grad(v_b)) * dOmegaB
F -= ufl.inner(f_A, v_a) * dOmegaA
F -= ufl.inner(f_B, v_b) * dOmegaB
alpha = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(alpha_c))
K_a = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(1.0))
K_b = dolfinx.fem.Constant(
mesh,
dolfinx.default_scalar_type(n_sorption * d ** ((n_sorption - 1) / n_sorption)),
)
# Add penalty term for alpha/2 * equation**2 dGamma
tol = dolfinx.fem.Constant(mesh, 1e2 * np.finfo(dolfinx.default_scalar_type).eps)
cond_b = ufl.gt(abs(tr_b(u_b)), tol)
u_b_padded = ufl.conditional(cond_b, tr_b(u_b), tol)
equation = tr_a(u_a) - K_a * (u_b_padded / K_b) ** n_sorption
deq_du_a = ufl.diff(equation, u_a) * tr_a(v_a)
deq_du_b = ufl.diff(equation, u_b) * tr_b(v_b)
if penalty_mode == PenaltyMode.QUADRATIC:
F += alpha * equation * (deq_du_a + deq_du_b) * dGamma
elif penalty_mode == PenaltyMode.BC:
F += alpha * equation * (tr_a(v_a) - tr_b(v_b)) * dGamma
else:
raise ValueError(f"Unknown penalty mode: {penalty_mode}")
F_blocked = ufl.extract_blocks(F)
petsc_options = {
"snes_type": "newtonls",
"snes_linesearch_type": "none",
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"snes_monitor": None,
"snes_error_if_not_converged": True,
"ksp_error_if_not_converged": True,
"ksp_monitor": None,
"snes_atol": 1e-10,
"snes_rtol": 1e-10,
}
solver = dolfinx.fem.petsc.NonlinearProblem(
F_blocked,
[u_a, u_b],
bcs=bcs,
petsc_options=petsc_options,
petsc_options_prefix="penalty_solver_",
entity_maps=[map_A, map_B],
kind="mpi",
)
solver.solve()
with dolfinx.io.VTXWriter(omega_A.comm, "u_A.bp", [u_a]) as bp:
bp.write(0.0)
with dolfinx.io.VTXWriter(omega_B.comm, "u_B.bp", [u_b]) as bp:
bp.write(0.0)
ub_ex = n_sorption * x[0] ** (1 / n_sorption) / d ** (1 / n_sorption - 1)
ua_ex = x[0]
u_b_error = ufl.inner(u_b - ub_ex, u_b - ub_ex) * dOmegaB
L2_b = dolfinx.fem.form(u_b_error, entity_maps=[map_B])
print("L2 error u_b:", np.sqrt(dolfinx.fem.assemble_scalar(L2_b)))
u_a_error = ufl.inner(u_a - ua_ex, u_a - ua_ex) * dOmegaA
L2_a = dolfinx.fem.form(u_a_error, entity_maps=[map_A])
print("L2 error u_a:", np.sqrt(dolfinx.fem.assemble_scalar(L2_a)))
if __name__ == "__main__":
parser = argparse.ArgumentParser(
description="Solve penalty problem.",
formatter_class=argparse.ArgumentDefaultsHelpFormatter,
)
parser.add_argument(
"--penalty_mode",
type=str,
choices=[e.name for e in PenaltyMode],
default=PenaltyMode.QUADRATIC.name,
help="Penalty mode to use: BC or QUADRATIC.",
)
parser.add_argument(
"--Nx", type=int, default=10, help="Number of elements in the x direction."
)
parser.add_argument(
"--Ny", type=int, default=10, help="Number of elements in the y direction."
)
parser.add_argument(
"--d",
type=float,
default=0.1,
help="Parameter d for the penalty problem.",
)
parser.add_argument(
"--alpha", type=float, default=1000.0, help="Penalty parameter alpha."
)
args = parser.parse_args()
penalty_mode = PenaltyMode[args.penalty_mode]
solve_penalty_problem(
penalty_mode=penalty_mode, d=args.d, Nx=args.Nx, Ny=args.Ny, alpha_c=args.alpha
)
Running this with
python3 mwe_penalty.py --penalty_mode=QUADRATIC --d 0.1 --Nx=100 --Ny=100 --alpha=10000
yields
0 SNES Function norm 7.079204673121e+01
Residual norms for penalty_solver_ solve.
0 KSP Residual norm 7.079204673121e+01
1 KSP Residual norm 1.538620576106e-13
1 SNES Function norm 4.460940922661e+02
Residual norms for penalty_solver_ solve.
0 KSP Residual norm 4.460940922661e+02
1 KSP Residual norm 1.349764383516e-13
2 SNES Function norm 9.981139751230e-02
Residual norms for penalty_solver_ solve.
0 KSP Residual norm 9.981139751230e-02
1 KSP Residual norm 4.811261996135e-14
3 SNES Function norm 1.805742841077e+02
Residual norms for penalty_solver_ solve.
0 KSP Residual norm 1.805742841077e+02
1 KSP Residual norm 3.182108548435e-14
4 SNES Function norm 8.203457721008e-05
Residual norms for penalty_solver_ solve.
0 KSP Residual norm 8.203457721008e-05
1 KSP Residual norm 5.142089122414e-17
5 SNES Function norm 3.109902078367e-04
Residual norms for penalty_solver_ solve.
0 KSP Residual norm 3.109902078367e-04
1 KSP Residual norm 5.577500438612e-20
6 SNES Function norm 7.057244655731e-11
L2 error u_b: 0.0005641090470741396
L2 error u_a: 0.00014560136774224966
while using the BC modification (see festim 2 overleaf for full description of the variational form), which is the one that is implemented in FESTIM (and Moose):
python3 mwe_penalty.py --penalty_mode=BC --d 0.1 --Nx=100 --Ny=100 --alpha=10000
0 SNES Function norm 7.079204673592e+01
Residual norms for penalty_solver_ solve.
0 KSP Residual norm 7.079204673592e+01
1 KSP Residual norm 1.503161040347e-13
1 SNES Function norm 6.148980776698e+02
Residual norms for penalty_solver_ solve.
0 KSP Residual norm 6.148980776698e+02
1 KSP Residual norm 1.049846983884e-13
Traceback (most recent call last):
File "/root/shared/mwe_penalty.py", line 234, in <module>
solve_penalty_problem(
File "/root/shared/mwe_penalty.py", line 185, in solve_penalty_problem
solver.solve()
File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/petsc.py", line 1399, in solve
self.solver.solve(None, self.x)
File "petsc4py/PETSc/SNES.pyx", line 1738, in petsc4py.PETSc.SNES.solve
petsc4py.PETSc.Error: error code 91
[0] SNESSolve() at /usr/local/petsc/src/snes/interface/snes.c:4846
[0] SNESSolve_NEWTONLS() at /usr/local/petsc/src/snes/impls/ls/ls.c:240
[0] SNESSolve has not converged due to Nan or Inf norm
Thanks for this @jorgensd !
Is it possible to use this method for the interface condition with the current implementation of the HydrogenTransportProblemDiscontinuous class, which involves having formulations per subdomain? Can we evaluate formulations per subdomain, assemble the forms, but also add the interface terms (F += alpha * equation * (deq_du_a + deq_du_b) * dGamma) to the global formulation?
Thanks for this @jorgensd !
Is it possible to use this method for the interface condition with the current implementation of the
HydrogenTransportProblemDiscontinuousclass, which involves having formulations per subdomain? Can we evaluate formulations per subdomain, assemble the forms, but also add the interface terms (F += alpha * equation * (deq_du_a + deq_du_b) * dGamma) to the global formulation?
@jhdark, yes it is very much doable.
Note that the current implementation is implemented as: https://github.com/festim-dev/FESTIM/blob/fenicsx/src/festim/hydrogen_transport_problem.py#L1519-L1531 which corresponds to a version of:
Note that for Sievert-Sievert or Henry-Henry interfaces, where you have u_A/K_A = u_B/K_B, the approaches are the same, as $n=1$
If we want to implement the latter properly in FESTIm, we need to ensure that we write the penalty such that the equality is written as $u_A - K_A\left(\frac{u_B}{K_B}\right)^n$, where $n=1$ or $\frac{1}{2}$, meaning that the Sievert material should always be $u_A$ (correct me if i mess up which material is which. I am basing it on
@jorgensd we realised with @hhy2022 that we really need to implement this as the Sivert/henry interface simply don't work. I'll start working on it
Now @RemDelaporteMathurin has a prototype with the new method, after an incredibly efficient live coding session.
It's happening there: #1030
@jorgensd quick question: will this work in cylindrical coordinates?