A better way to restrict this is to follow the suggestion of @MiroK from:
which ensures a consistent orientation. I.e. mark all cells on the side you would like to have a “+” restriction on with a larger value than on the other side of the interface.
Instead of defining two different measures, I would suggest using a single measure, like
dx = Measure("dx", domain=mesh, subdomain_data=sf)
dx_shape = dx(shape_id)
dx_square = dx(square_id)
So let’s go through this step by step.
If we define:
dx = Measure("dx", domain=mesh, subdomain_data=sf)
dx_shape = dx(shape_id)
dx_square = dx(square_id)
dS = Measure("dS", domain=mesh, subdomain_data=mf) # all interior facets
ds_outer = Measure("ds", domain=mesh) # outer boundary
dS_shape = dS(INTERFACE_ID)
dS_I_shape = dS(INT_SHAPE_ID)
dS_I_square = dS(INT_SQUARE_ID)
if shape_id > square_id:
shape_label = "+"
square_label = "-"
else:
shape_label = "-"
square_label = "+"
and then compare the first part of your forms, i.e.
F_0 = (
u * u.dx(i) * nu_u.dx(i) + f * nu_u
) * dx_shape + (
u.dx(i) * nu_u.dx(i) + f * nu_u
) * dx_square
a0 = assemble(F_0)
# option 2 (working)
F_01 = ufl_conditional_form(mesh,
sf,
u * u.dx(i) * nu_u.dx(i) + f * nu_u,
u.dx(i) * nu_u.dx(i) + f * nu_u,
shape_id,
square_id
) * dx
a0_1 = assemble(F_01)
np.testing.assert_allclose(a0.get_local(), a0_1.get_local())
I’ve also tried to create the full forms:
F_0 = (
u * u.dx(i) * nu_u.dx(i) + f * nu_u
) * dx_shape + (
u.dx(i) * nu_u.dx(i) + f * nu_u
) * dx_square \
- n[i] * u.dx(i) * nu_u * ds_outer \
- n(square_label)[i] * (u(square_label)).dx(i) * nu_u(square_label) * dS_shape \
- n(shape_label)[i] * u(shape_label) * (u(shape_label)).dx(i) * nu_u(shape_label) * dS_shape
# option 2 (working)
F_01 = ufl_conditional_form(mesh,
sf,
u * u.dx(i) * nu_u.dx(i) + f * nu_u,
u.dx(i) * nu_u.dx(i) + f * nu_u,
shape_id,
square_id
) * dx \
- n[i] * u.dx(i) * nu_u * ds_outer \
- n(square_label)[i] * (u(square_label)).dx(i) * nu_u(square_label) * dS_shape \
- n(shape_label)[i] * u(shape_label) * (u(shape_label)).dx(i) * nu_u(shape_label) * dS_shape
and then solve the problem twice to see if I can spot the issue:
a = assemble(F)
a2 = assemble(F2)
np.testing.assert_allclose(a.get_local(), a2.get_local())
J = derivative(F, u)
A = assemble(J)
A2 = assemble(derivative(F2, u))
np.testing.assert_allclose(A.array(), A2.array(), atol=1e-13)
solve(F2 == 0, u, solver_parameters={"newton_solver": {"linear_solver": "mumps"}})
d = u.copy(deepcopy=True)
u.vector()[:] = 1.0
solve(F == 0, u, solver_parameters={"newton_solver": {"linear_solver": "mumps"}})
and then I get that the jacobians are the same at the first iteration, but they seem to deviate afterwards, as I get the two solver outputs:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.073e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 3.098e-02 (tol = 1.000e-10) r (rel) = 2.886e-05 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 2.035e-03 (tol = 1.000e-10) r (rel) = 1.895e-06 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 8.184e-06 (tol = 1.000e-10) r (rel) = 7.625e-09 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 1.631e-10 (tol = 1.000e-10) r (rel) = 1.520e-13 (tol = 1.000e-09)
Newton solver finished in 4 iterations and 4 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.073e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 6.052e-01 (tol = 1.000e-10) r (rel) = 5.639e-04 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 3.228e-03 (tol = 1.000e-10) r (rel) = 3.007e-06 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 8.185e-06 (tol = 1.000e-10) r (rel) = 7.625e-09 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 1.631e-10 (tol = 1.000e-10) r (rel) = 1.520e-13 (tol = 1.000e-09)
Newton solver finished in 4 iterations and 4 linear solver iterations.
where we see that the initial residual is the same.
However, if you lessen the restriction on the jacobians, you will get that:
Traceback (most recent call last):
File "/root/shared/mwe_mekong.py", line 286, in <module>
np.testing.assert_allclose(A.array(), A2.array())
File "/usr/lib/python3/dist-packages/numpy/testing/_private/utils.py", line 1530, in assert_allclose
assert_array_compare(compare, actual, desired, err_msg=str(err_msg),
File "/usr/lib/python3/dist-packages/numpy/testing/_private/utils.py", line 844, in assert_array_compare
raise AssertionError(msg)
AssertionError:
Not equal to tolerance rtol=1e-07, atol=0
Mismatched elements: 2753 / 34152336 (0.00806%)
Max absolute difference: 1.52100554e-14
Max relative difference: 82.
x: array([[37.001268, -4.704807, -4.403054, ..., 0. , 0. ,
0. ],
[-4.704807, 35.352108, -3.959018, ..., 0. , 0. ,...
y: array([[37.001268, -4.704807, -4.403054, ..., 0. , 0. ,
0. ],
[-4.704807, 35.352108, -3.959018, ..., 0. , 0. ,...
which is floating error accumulated in the conditional logic.
However, I think there might be a bug in FFC here, and I have no further time to play around with it. As I cannot reproduce the descrepancy in DOLFINx, this is likely something in legacy FFC:
from mpi4py import MPI
from dolfinx import default_scalar_type
from dolfinx.fem import functionspace, Function
from ufl import Measure, FacetNormal, CellDiameter, conditional, TestFunction, indices, eq, derivative
from dolfinx.io.gmsh import read_from_msh
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.fem import assemble_scalar, IntegralType, compute_integration_domains, Constant
from dolfinx.mesh import compute_incident_entities, meshtags, GhostMode
from scifem import compute_interface_data
import gmsh
import numpy as np
from dolfinx.cpp.mesh import create_cell_partitioner
# =========================================================
# 1. MESH GENERATION (Gmsh, with embedded interface)
# =========================================================
gmsh.initialize()
gmsh.model.add("model")
L, H = 1.0, 1.0
resolution = 0.05
shape_id = 1 # inner square (shape)
square_id = 2 # outer region
# outer boundary
p_out = [gmsh.model.geo.addPoint(x, y, 0)
for x, y in [(0,0),(L,0),(L,H),(0,H)]]
l_out = [gmsh.model.geo.addLine(p_out[i], p_out[(i+1)%4]) for i in range(4)]
outer_loop = gmsh.model.geo.addCurveLoop(l_out)
# inner shape (square)
shape_pts = [[0.3,0.3],[0.7,0.3],[0.7,0.7],[0.3,0.7]]
p_in = [gmsh.model.geo.addPoint(x, y, 0) for x, y in shape_pts]
l_in = [gmsh.model.geo.addLine(p_in[i], p_in[(i+1)%4]) for i in range(4)]
inner_loop = gmsh.model.geo.addCurveLoop(l_in)
surf_square = gmsh.model.geo.addPlaneSurface([outer_loop, inner_loop])
surf_shape = gmsh.model.geo.addPlaneSurface([inner_loop])
gmsh.model.geo.synchronize()
gmsh.model.mesh.embed(1, l_in, 2, surf_square) # embed interface in outer region
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(2, [surf_shape], shape_id)
gmsh.model.addPhysicalGroup(2, [surf_square], square_id)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", resolution)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", resolution)
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
gmsh.finalize()
gmode = GhostMode.shared_facet
partitioner = create_cell_partitioner(gmode, max_facet_to_cell_links=2)
mesh_data = read_from_msh("mesh.msh",MPI.COMM_WORLD, 0, gdim=2, partitioner=partitioner)
mesh = mesh_data.mesh
sf = mesh_data.cell_tags
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
num_facets_on_proc = mesh.topology.index_map(mesh.topology.dim-1).size_local + mesh.topology.index_map(mesh.topology.dim-1).num_ghosts
INT_SHAPE_ID = 1
INT_SQUARE_ID = 2
INT_INTERFACE_ID = INT_SHAPE_ID + INT_SQUARE_ID
facet_markers = np.zeros(num_facets_on_proc, dtype=np.int32)
facet_markers[compute_incident_entities(mesh.topology, sf.find(square_id), mesh.topology.dim, mesh.topology.dim-1)] += INT_SQUARE_ID
facet_markers[compute_incident_entities(mesh.topology, sf.find(shape_id), mesh.topology.dim, mesh.topology.dim-1)] += INT_SHAPE_ID
ft = meshtags(mesh, mesh.topology.dim-1, np.arange(num_facets_on_proc, dtype=np.int32), facet_markers)
interface_entities = compute_interface_data(sf, ft.find(INT_INTERFACE_ID))
integral_data = [(INT_INTERFACE_ID, interface_entities.flatten())]
integral_data +=[(INT_SHAPE_ID, compute_integration_domains(
IntegralType.interior_facet,
mesh.topology,
ft.find(shape_id),))]
integral_data +=[(INT_SQUARE_ID, compute_integration_domains(
IntegralType.interior_facet,
mesh.topology,
ft.find(square_id),))]
if shape_id > square_id:
shape_label = "-"
square_label = "+"
else:
shape_label = "+"
square_label = "-"
metadata = {"quadrature_degree": 15, "quadrature_scheme": "default"}
dx = Measure("dx", domain=mesh, subdomain_data=sf, metadata=metadata)
dx_shape = dx(shape_id)
dx_square = dx(square_id)
dS = Measure("dS", domain=mesh, subdomain_data=integral_data, metadata=metadata) # all interior facets
ds_outer = Measure("ds", domain=mesh, metadata=metadata) # outer boundary
dS_shape = dS(INT_INTERFACE_ID)
dS_I_shape = dS(INT_SHAPE_ID)
dS_I_square = dS(INT_SQUARE_ID)
if shape_id > square_id:
shape_label = "+"
square_label = "-"
else:
shape_label = "-"
square_label = "+"
print(f'label_shape = {shape_label} {shape_id=}')
print(f'label_square = {square_label} {square_id=}')
# =========================================================
# 6. FUNCTION SPACE
# =========================================================
Q = functionspace(mesh, ("DG", 2))
u = Function(Q)
nu_u = TestFunction(Q)
f = Function(Q)
u_exact = Function(Q)
# =========================================================
# 7. INTERPOLATION (subdomain-wise)
# =========================================================
def f_shape(x): return 2.0 * (-1 + x[0]**2 + 10*x[1]**2)
def f_square(x): return np.full_like(x[0], 6.0)
def u_sh(x): return 1 + x[0]**2 - 2*x[1]**2
def u_sq(x): return 1 + x[0]**2 + 2*x[1]**2
f.interpolate(f_shape, cells0=sf.find(shape_id))
f.interpolate(f_square, cells0=sf.find(square_id))
f.x.scatter_forward() # ensure f is up to date on all processes for assembly
u_exact.interpolate(u_sh, cells0=sf.find(shape_id))
u_exact.interpolate(u_sq, cells0=sf.find(square_id))
u_exact.x.scatter_forward() # ensure u_exact is up to date on all processes for assembly
# =========================================================
# 8. PARAMETERS AND GEOMETRIC QUANTITIES
# =========================================================
alpha = Constant(mesh,default_scalar_type(100))
n = FacetNormal(mesh)
h = CellDiameter(mesh) # NB: must be restricted ('.+') in dS integrals
# because CellEdgeVectors is discontinuous
i = indices(1)[0] # free UFL index (Einstein summation)
u.x.array[:] = 1.0 # initial guess (for computing h on dS)
h_loc = mesh.h(mesh.topology.dim, np.arange(mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32))
r_min = mesh.comm.allreduce(h_loc.min(), op=MPI.MIN)
# =========================================================
# 9. DG OPERATORS
# =========================================================
def average(expr):
"""Arithmetic average across an interior facet."""
return (expr('+') + expr('-')) / 2
def jump(f, n):
"""Vectorial jump of scalar f: f(+) n(+) + f(-) n(-)."""
return f('+')*n('+') + f('-')*n('-')
# =========================================================
# 10. VARIATIONAL FORM (exact transcription of original)
# =========================================================
# ---------- F_0 : bulk terms (from IBP) ----------
# shape region : nonlinear equation u * u_{,i} * (nu_u)_{,i} + f * nu_u
# square region : linear equation u_{,i} * (nu_u)_{,i} + f * nu_u
# Boundary/interface flux terms:
# outer boundary : - n_i u_{,i} nu_u
# interface, square side : - n_i u_{,i} nu_u (from linear BVP)
# interface, shape side : - n_i u u_{,i} nu_u (from nonlinear BVP)
def ufl_conditional_form(mesh, sf, form_a, form_b, tag_a, tag_b):
Q = functionspace(mesh, ("DG", 0))
cell_tag = Function(Q)
cell_tag.x.array[:] = sf.values # DG0 dofs are ordered by cell index
# build f as a pure UFL expression
result = conditional(
eq(cell_tag, tag_a),
form_a,
conditional(eq(cell_tag, tag_b),
form_b,
0.0
)
)
return result
# option 1 (not working)
i = indices(1)[0] # free UFL index (Einstein summation)
F_0 = (
u * u.dx(i) * nu_u.dx(i) + f * nu_u
) * dx_shape + (
u.dx(i) * nu_u.dx(i) + f * nu_u
) * dx_square
# option 2 (working)
F_01 = ufl_conditional_form(mesh,
sf,
u * u.dx(i) * nu_u.dx(i) + f * nu_u,
u.dx(i) * nu_u.dx(i) + f * nu_u,
shape_id,
square_id
) * dx
F_other = - n[i] * u.dx(i) * nu_u * ds_outer \
- n(square_label)[i] * (u(square_label)).dx(i) * nu_u(square_label) * dS_shape \
- n(shape_label)[i] * u(shape_label) * (u(shape_label)).dx(i) * nu_u(shape_label) * dS_shape
# ---------- F_I : interior DG (SIPG-style) ----------
# square interior: standard symmetric term -avg(grad u) · jump(nu_u, n)
# shape interior: nonlinear analog -avg(grad u) · jump(u nu_u, n)
# penalty on both: (alpha/h) jump(u,n) · jump(nu_u,n)
#
# NB: h must be restricted on dS integrals (CellEdgeVectors is discontinuous)
F_I = \
- average(u.dx(i)) * jump(nu_u, n)[i] * dS_I_square \
- average(u.dx(i)) * jump(u * nu_u, n)[i] * dS_I_shape \
+ alpha / r_min * (
jump(u, n)[i] * jump(nu_u, n)[i]
) * dS_I_square \
+ alpha / r_min * (
jump(u, n)[i] * jump(nu_u, n)[i]
) * dS_I_shape
# ---------- F_b : Nitsche boundary penalisation ----------
# outer boundary : (alpha/h) (u - u_exact) nu_u ds
# interface, square side : (alpha/h) (u - u_exact) nu_u dS
# interface, shape side : (alpha/h) (u - u_exact) nu_u dS
F_b = alpha / h * (u - u_exact) * nu_u * ds_outer \
+ alpha / r_min * (
(u(square_label) - u_exact(square_label)) * nu_u(square_label)
) * dS_shape \
+ alpha / r_min * (
(u(shape_label) - u_exact(shape_label)) * nu_u(shape_label)
) * dS_shape
F = F_0 + F_other + F_I + F_b
F2 = F_01 + F_other+ F_I + F_b
# =========================================================
# 11. SOLVE
# =========================================================
np.random.seed(0)
u.x.array[:] = np.random.rand(u.x.array.size) # random initial guess
from dolfinx.fem import assemble_vector, form, assemble_matrix
a = assemble_vector(form(F))
a2 = assemble_vector(form(F2))
np.testing.assert_allclose(a.array, a2.array)
J = derivative(F, u)
A = assemble_matrix(form(J))
A2 = assemble_matrix(form(derivative(F2, u)))
np.testing.assert_allclose(A.data, A2.data)
u.x.array[:] = 1.0
opts = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps", "snes_linesearch_type": "none", "snes_rtol": 1e-10,
"snes_error_if_not_converged": True, "snes_monitor": None}
solver = NonlinearProblem(F, u, petsc_options_prefix="org_solver", petsc_options=opts)
solver.solve()
d_array = u.x.array.copy()
u.x.array[:] = 1.0
solver_2 = NonlinearProblem(F2, u, petsc_options_prefix="option_solver", petsc_options=opts)
solver_2.solve()
print(np.max(np.abs(d_array - u.x.array)))
print(f'error_shape = {mesh.comm.allreduce(assemble_scalar(form((u - u_exact)**2 * dx_shape)), op=MPI.SUM)}')
print(f'error_square = {mesh.comm.allreduce(assemble_scalar(form((u - u_exact)**2 * dx_square)), op=MPI.SUM)}')
yielding:
label_shape = - shape_id=1
label_square = + square_id=2
0 SNES Function norm 1.073369114053e+03
1 SNES Function norm 4.162191337893e-01
2 SNES Function norm 2.028425869492e-03
3 SNES Function norm 8.088894392143e-06
4 SNES Function norm 1.585554376335e-10
0 SNES Function norm 1.073369114053e+03
1 SNES Function norm 4.162191337893e-01
2 SNES Function norm 2.028425869492e-03
3 SNES Function norm 8.088894392143e-06
4 SNES Function norm 1.585554376335e-10
0.0
error_shape = 8.189560766580326e-09
error_square = 8.465276585475065e-08