WIP Codim coupling - #1216
Draft
RemDelaporteMathurin wants to merge 14 commits into
Draft
Conversation
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #1216 +/- ##
==========================================
+ Coverage 95.19% 95.36% +0.16%
==========================================
Files 53 54 +1
Lines 3952 4355 +403
==========================================
+ Hits 3762 4153 +391
- Misses 190 202 +12 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Collaborator
Author
|
@jorgensd this is the PR of the codim implementation in FESTIM (currently limited to codim <=1 but I don't see why we couldn't do codim 2). Note that it doesn't use |
Contributor
|
It appears that the current internal subdomain formulation doesn't handle transient cases correctly; if we try to run the MWE below (a FESTIM adaptation of the pure dolfinx code found at https://github.com/ee-nn/FESTIM/blob/multi-dimensional-coupling/scratch_work/interface_transport_custom_solver.py) then we get a traceback which I believe points to the Full MWE"""
Modified from https://gist.github.com/RemDelaporteMathurin/d1a678b6b7439e339c8471e97cd31a39
Demonstrates a FESTIM codimension-1 internal interface with diffusion and coupling.
"""
from mpi4py import MPI
import dolfinx.mesh
import numpy as np
import festim as F
# Parameters
L = 4.0
x_int = L / 2
D_left = 2.0
D_right = 1.5
D_int = 1.0
k1 = 1.0
k2 = 1.0
k3 = 1.0
k4 = 1.0
c_int_max = 1.0
# Transient settings
dt = 0.1
T = 10.0
mesh = dolfinx.mesh.create_rectangle(
MPI.COMM_WORLD,
[np.array([0.0, 0.0]), np.array([L, 1.0])],
[20, 10],
cell_type=dolfinx.mesh.CellType.quadrilateral,
)
eps = 1e-12
left = F.VolumeSubdomain(
id=1,
material=F.Material(D_0=D_left, E_D=0.0),
locator=lambda x: x[0] <= x_int + eps,
name="left",
)
right = F.VolumeSubdomain(
id=2,
material=F.Material(D_0=D_right, E_D=0.0),
locator=lambda x: x[0] >= x_int - eps,
name="right",
)
gamma = F.VolumeSubdomain(
id=3,
material=F.Material(D_0=D_int, E_D=0.0),
dim=mesh.topology.dim - 1,
locator=lambda x: np.isclose(x[0], x_int),
name="interface",
)
left_boundary = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 0.0))
right_boundary = F.SurfaceSubdomain(id=5, locator=lambda x: np.isclose(x[0], L))
H_left = F.Species("c_left", subdomains=[left])
H_right = F.Species("c_right", subdomains=[right])
H_int = F.Species("c_int", subdomains=[gamma])
# Flux from the external left boundary into the left bulk.
left_boundary_flux = F.ParticleFluxBC(
subdomain=left_boundary,
species=H_left,
value=0.5,
)
# Interface coupling: bulk fluxes into the codim-1 interface, plus matching sources
# in the interface equation.
interface_sources = [
F.ParticleSource(
value=lambda c_int, c_left: -(
k1 * c_left * (1.0 - c_int / c_int_max) - k2 * c_int
),
species=H_int,
volume=gamma,
species_dependent_value={"c_int": H_int, "c_left": H_left},
),
F.ParticleSource(
value=lambda c_int, c_right: -(
k3 * c_right * (1.0 - c_int / c_int_max) - k4 * c_int
),
species=H_int,
volume=gamma,
species_dependent_value={"c_int": H_int, "c_right": H_right},
),
]
interface_fluxes = [
F.ParticleFluxBC(
subdomain=gamma,
species=H_left,
value=lambda c_int, c_left: k1 * c_left * (1.0 - c_int / c_int_max)
- k2 * c_int,
species_dependent_value={"c_int": H_int, "c_left": H_left},
),
F.ParticleFluxBC(
subdomain=gamma,
species=H_right,
value=lambda c_int, c_right: k3 * c_right * (1.0 - c_int / c_int_max)
- k4 * c_int,
species_dependent_value={"c_int": H_int, "c_right": H_right},
),
]
right_dirichlet = F.FixedConcentrationBC(
subdomain=right_boundary,
value=0.0,
species=H_right,
)
model = F.HydrogenTransportProblemDiscontinuous(
mesh=F.Mesh(mesh),
species=[H_left, H_right, H_int],
subdomains=[left, right, gamma, left_boundary, right_boundary],
sources=[*interface_sources],
boundary_conditions=[left_boundary_flux, *interface_fluxes, right_dirichlet],
temperature=500,
settings=F.Settings(
atol=1e-10,
rtol=1e-10,
transient=True,
final_time=T,
stepsize=dt,
),
exports=[
F.VTXSpeciesExport(filename="results/c_left.bp", field=H_left, subdomain=left),
F.VTXSpeciesExport(
filename="results/c_right.bp", field=H_right, subdomain=right
),
F.VTXSpeciesExport(filename="results/c_int.bp", field=H_int, subdomain=gamma),
],
)
model.initialise()
model.run()
c_left = H_left.subdomain_to_post_processing_solution[left].x.array
c_right = H_right.subdomain_to_post_processing_solution[right].x.array
c_int = H_int.subdomain_to_post_processing_solution[gamma].x.array
print("c_left range", c_left.min(), c_left.max())
print("c_right range", c_right.min(), c_right.max())
print("c_int range", c_int.min(), c_int.max())Full Traceback File "~/mwe_internal_interface.py", line 138, in <module>
model.initialise()
~~~~~~~~~~~~~~~~^^
File "~/src/festim/hydrogen_transport_problem.py", line 1398, in initialise
self.create_formulation()
~~~~~~~~~~~~~~~~~~~~~~~^^
File "~/src/festim/hydrogen_transport_problem.py", line 2284, in create_formulation
dolfinx.fem.form(g, entity_maps=entity_maps, jit_options=jit_options)
~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 476, in form
return _create_form(form)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 472, in _create_form
return list(map(lambda sub_form: _create_form(sub_form), form))
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 472, in <lambda>
return list(map(lambda sub_form: _create_form(sub_form), form))
~~~~~~~~~~~~^^^^^^^^^^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 468, in _create_form
return _form(form)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 384, in _form
ufcx_form, module, code = jit.ffcx_jit(
~~~~~~~~~~~~^
comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/jit.py", line 61, in mpi_jit
return local_jit(*args, **kwargs)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/jit.py", line 221, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/codegeneration/jit.py", line 244, in compile_forms
raise e
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/codegeneration/jit.py", line 224, in compile_forms
impl = _compile_objects(
decl,
...<9 lines>...
visualise=visualise,
)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/codegeneration/jit.py", line 349, in _compile_objects
code, _ = ffcx.compiler.compile_ufl_objects(
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
ufl_objects, namespace=module_name, options=options, visualise=visualise
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/compiler.py", line 116, in compile_ufl_objects
ir = compute_ir(analysis, _object_names, _namespace, options, visualise)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/ir/representation.py", line 261, in compute_ir
_compute_integral_ir(
~~~~~~~~~~~~~~~~~~~~^
fd,
^^^
...<4 lines>...
visualise,
^^^^^^^^^^
)
^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/ir/representation.py", line 470, in _compute_integral_ir
integral_ir = compute_integral_ir(
itg_data.domain.ufl_cell(),
...<5 lines>...
visualise,
)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/ir/integral.py", line 430, in compute_integral_ir
) = _compute_integral_ir(
~~~~~~~~~~~~~~~~~~~~^
expression,
^^^^^^^^^^^
...<7 lines>...
p,
^^
)
^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/ir/integral.py", line 178, in _compute_integral_ir
mt_table_reference = build_optimized_tables(
quadrature_rule,
...<8 lines>...
atol=p["table_atol"],
)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/ir/elementtables.py", line 569, in build_optimized_tables
tbl = clamp_table_small_numbers(t["array"], rtol=rtol, atol=atol)
^
UnboundLocalError: cannot access local variable 't' where it is not associated with a value
Although this MWE used an internal subdomain I believe that this is also an issue for surface subdomains since its use of |
RemDelaporteMathurin
marked this pull request as draft
August 3, 2026 17:45
Comment on lines
+152
to
+154
| # NOTE: the shape is the *geometric* dimension, not the topological one: on a | ||
| # codim-1 subdomain the mesh is a manifold (eg. a line in 2D) whose cells are | ||
| # 1D but whose points, and therefore velocities and gradients, are ambient |
Collaborator
Author
There was a problem hiding this comment.
Suggested change
| # NOTE: the shape is the *geometric* dimension, not the topological one: on a | |
| # codim-1 subdomain the mesh is a manifold (eg. a line in 2D) whose cells are | |
| # 1D but whose points, and therefore velocities and gradients, are ambient |
| if subdomain in solutions: | ||
| return solutions[subdomain] | ||
| if len(solutions) == 1: | ||
| return next(iter(solutions.values())) |
Collaborator
Author
There was a problem hiding this comment.
this is a way to get around the fact that solutions.values() is not a list
ee-nn
reviewed
Aug 3, 2026
|
|
||
| def convert_advection_term_to_fenics_objects(self): | ||
| """For each advection term convert the input value.""" | ||
|
|
Contributor
There was a problem hiding this comment.
Should this be added for consistency?
Suggested change
| t = self.subdomain_time(advec_term.subdomain) |
and then use t = t in advection_term.velocity_convert_input_value() on the line below
Contributor
There was a problem hiding this comment.
Actually I don't think this is needed because the time dependence isn't carried over the same way here
ee-nn
reviewed
Aug 3, 2026
| if isinstance(advec_term, AdvectionTerm): | ||
| for spe in advec_term.species: | ||
| V = spe.subdomain_to_function_space[advec_term.subdomain] | ||
| advec_term.velocity.convert_input_value(function_space=V, t=self.t) |
Contributor
There was a problem hiding this comment.
Suggested change
| advec_term.velocity.convert_input_value(function_space=V, t=t) |
ee-nn
reviewed
Aug 3, 2026
| if source.volume not in self.manifold_to_volumes: | ||
| return False | ||
| return not self.foreign_species(source, source.volume) | ||
|
|
Contributor
There was a problem hiding this comment.
Adding traps to a co-dimensional subdomain currently runs into a form compilation error because of a mesh mismatch (similar to the error that was fixed in a99b5a9). An updated create_implicit_species_value_fenics() function that overrides the one in the base HydrogenTransportProblem should fix this issue:
Suggested change
| def create_implicit_species_value_fenics(self): | |
| """For each implicit species, create the value_fenics. | |
| The density of an implicit species consumed by a reaction on a manifold | |
| subdomain appears in an integral over that manifold's submesh, so like every | |
| other coefficient of such an integral it has to be built there rather than on | |
| the parent mesh (see :meth:`create_submesh_time_constants`). | |
| """ | |
| species_to_mesh = {} | |
| for reaction in self.reactions: | |
| volume = reaction.volume | |
| if volume is not None and volume.codim(self.mesh.vdim) == 1: | |
| mesh, t = volume.submesh, self.subdomain_time(volume) | |
| else: | |
| mesh, t = self.mesh.mesh, self.t | |
| for reactant in reaction.reactant: | |
| if not isinstance(reactant, _species.ImplicitSpecies): | |
| continue | |
| # an implicit species shared by reactions on different meshes would be | |
| # built twice and keep only the last one, silently leaving a foreign | |
| # terminal in one of the two integrals | |
| previous = species_to_mesh.setdefault(id(reactant), mesh) | |
| if previous is not mesh: | |
| raise NotImplementedError( | |
| f"implicit species {reactant.name} is used by reactions on a " | |
| "codim-1 subdomain and on another subdomain, which are " | |
| "integrated over different meshes. Declare one implicit " | |
| "species per subdomain." | |
| ) | |
| reactant.create_value_fenics(mesh=mesh, t=t) |
Example Usage
"""
Modified from https://gist.github.com/RemDelaporteMathurin/d1a678b6b7439e339c8471e97cd31a39
Demonstrates a FESTIM codimension-1 internal interface with diffusion, coupling to the
bulk on both sides, and a trapping site living on the interface itself.
"""
from mpi4py import MPI
import dolfinx.mesh
import numpy as np
import festim as F
# Parameters
L = 4.0
x_int = L / 2
D_left = 2.0
D_right = 1.5
D_int = 1.0
k1 = 1.0
k2 = 1.0
k3 = 1.0
k4 = 1.0
c_int_max = 1.0
# Trapping on the interface. n_trap is a density of trap sites in the same units as
# c_int (see the note on units at the bottom); k_trap and p_detrap are the trapping and
# detrapping rate constants, taken temperature-independent here (E_k = E_p = 0).
n_trap = 0.5
k_trap = 5.0
p_detrap = 1.0
# Transient settings
dt = 1e-2
T_final = 10.0
mesh = dolfinx.mesh.create_rectangle(
MPI.COMM_WORLD,
[np.array([0.0, 0.0]), np.array([L, 1.0])],
[20, 10],
cell_type=dolfinx.mesh.CellType.quadrilateral,
)
eps = 1e-12
left = F.VolumeSubdomain(
id=1,
material=F.Material(D_0=D_left, E_D=0.0),
locator=lambda x: x[0] <= x_int + eps,
name="left",
)
right = F.VolumeSubdomain(
id=2,
material=F.Material(D_0=D_right, E_D=0.0),
locator=lambda x: x[0] >= x_int - eps,
name="right",
)
gamma = F.VolumeSubdomain(
id=3,
material=F.Material(D_0=D_int, E_D=0.0),
dim=mesh.topology.dim - 1,
locator=lambda x: np.isclose(x[0], x_int),
name="interface",
)
left_boundary = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 0.0))
right_boundary = F.SurfaceSubdomain(id=5, locator=lambda x: np.isclose(x[0], L))
H_left = F.Species("c_left", subdomains=[left])
H_right = F.Species("c_right", subdomains=[right])
H_int = F.Species("c_int", subdomains=[gamma])
# The trapped population lives on gamma like c_int does, and is immobile: it gets a time
# derivative and the reaction terms, but no diffusion. F.Trap cannot be used here -- it
# builds its trapped species without a `subdomains` list, which the discontinuous
# problem needs -- so the species, the empty sites and the reaction are written out.
H_trapped = F.Species("c_trapped", mobile=False, subdomains=[gamma])
empty_sites = F.ImplicitSpecies(n=n_trap, others=[H_trapped], name="empty_sites")
# R = k_trap * c_int * (n_trap - c_trapped) - p_detrap * c_trapped
# enters the c_int equation as a sink and the c_trapped equation as a source
trapping = F.Reaction(
reactant=[H_int, empty_sites],
product=H_trapped,
k_0=k_trap,
E_k=0.0,
p_0=p_detrap,
E_p=0.0,
volume=gamma,
)
# Flux from the external left boundary into the left bulk.
left_boundary_flux = F.ParticleFluxBC(
subdomain=left_boundary,
species=H_left,
value=0.5,
)
# Interface coupling. J is the rate at which the interface *gains* particles from that
# side, so the interface gets +J as a source and the bulk gets -J as a flux.
def J_left(c_int, c_left):
return k1 * c_left * (1.0 - c_int / c_int_max) - k2 * c_int
def J_right(c_int, c_right):
return k3 * c_right * (1.0 - c_int / c_int_max) - k4 * c_int
interface_sources = [
F.ParticleSource(
value=J_left,
species=H_int,
volume=gamma,
species_dependent_value={"c_int": H_int, "c_left": H_left},
),
F.ParticleSource(
value=J_right,
species=H_int,
volume=gamma,
species_dependent_value={"c_int": H_int, "c_right": H_right},
),
]
interface_fluxes = [
F.ParticleFluxBC(
subdomain=gamma,
species=H_left,
value=lambda c_int, c_left: -J_left(c_int, c_left),
species_dependent_value={"c_int": H_int, "c_left": H_left},
),
F.ParticleFluxBC(
subdomain=gamma,
species=H_right,
value=lambda c_int, c_right: -J_right(c_int, c_right),
species_dependent_value={"c_int": H_int, "c_right": H_right},
),
]
right_dirichlet = F.FixedConcentrationBC(
subdomain=right_boundary,
value=0.0,
species=H_right,
)
model = F.HydrogenTransportProblemDiscontinuous(
mesh=F.Mesh(mesh),
species=[H_left, H_right, H_int, H_trapped],
subdomains=[left, right, gamma, left_boundary, right_boundary],
reactions=[trapping],
sources=[*interface_sources],
boundary_conditions=[left_boundary_flux, *interface_fluxes, right_dirichlet],
temperature=500,
settings=F.Settings(
atol=1e-10,
rtol=1e-10,
transient=True,
final_time=T_final,
stepsize=dt,
),
exports=[
F.VTXSpeciesExport(filename="results/c_left.bp", field=H_left, subdomain=left),
F.VTXSpeciesExport(
filename="results/c_right.bp", field=H_right, subdomain=right
),
F.VTXSpeciesExport(filename="results/c_int.bp", field=H_int, subdomain=gamma),
F.VTXSpeciesExport(
filename="results/c_trapped.bp", field=H_trapped, subdomain=gamma
),
],
)
model.initialise()
model.run()
c_left = H_left.subdomain_to_post_processing_solution[left].x.array
c_right = H_right.subdomain_to_post_processing_solution[right].x.array
c_int = H_int.subdomain_to_post_processing_solution[gamma].x.array
c_trapped = H_trapped.subdomain_to_post_processing_solution[gamma].x.array
print("c_left range", c_left.min(), c_left.max())
print("c_right range", c_right.min(), c_right.max())
print("c_int range", c_int.min(), c_int.max())
print("c_trapped range", c_trapped.min(), c_trapped.max())
Collaborator
Author
There was a problem hiding this comment.
A trap or a reaction or an ImplicitSpecies?
Contributor
There was a problem hiding this comment.
Added an example. I didn't use an F.Trap instance, but a combination of F.Reaction and F.ImplicitSpecies akin to the Implicit Trapping sites tutorial
Collaborator
Author
There was a problem hiding this comment.
Why would the implicit species in the bulk interact with the manifold one?
Collaborator
Author
There was a problem hiding this comment.
I see now that this is not the case but I still struggle to see why you would need to define two implicit species...
Collaborator
Author
There was a problem hiding this comment.
I understand
Contributor
There was a problem hiding this comment.
@RemDelaporteMathurin it looks like you already merged this in but just wanted to follow up on our discussion this morning where you mentioned the above issue may be cropping up from Reaction rather than ImplicitSpecies. Per your suggestion, I modified the example posted above to use a decay reaction on the interface:
Modified Example with Decay
from mpi4py import MPI
import dolfinx.mesh
import numpy as np
import festim as F
# Parameters
L = 4.0
x_int = L / 2
D_left = 2.0
D_right = 1.5
D_int = 1.0
k1 = 1.0
k2 = 1.0
k3 = 1.0
k4 = 1.0
c_int_max = 1.0
# Transient settings
dt = 1e-2
T_final = 10.0
mesh = dolfinx.mesh.create_rectangle(
MPI.COMM_WORLD,
[np.array([0.0, 0.0]), np.array([L, 1.0])],
[20, 10],
cell_type=dolfinx.mesh.CellType.quadrilateral,
)
eps = 1e-12
left = F.VolumeSubdomain(
id=1,
material=F.Material(D_0=D_left, E_D=0.0),
locator=lambda x: x[0] <= x_int + eps,
name="left",
)
right = F.VolumeSubdomain(
id=2,
material=F.Material(D_0=D_right, E_D=0.0),
locator=lambda x: x[0] >= x_int - eps,
name="right",
)
gamma = F.VolumeSubdomain(
id=3,
material=F.Material(D_0=D_int, E_D=0.0),
dim=mesh.topology.dim - 1,
locator=lambda x: np.isclose(x[0], x_int),
name="interface",
)
left_boundary = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 0.0))
right_boundary = F.SurfaceSubdomain(id=5, locator=lambda x: np.isclose(x[0], L))
H_left = F.Species("c_left", subdomains=[left])
H_right = F.Species("c_right", subdomains=[right])
H_int = F.Species("c_int", subdomains=[gamma])
decay = F.Reaction(reactant=[H_int], k_0=1, E_k=0, volume=gamma)
# Flux from the external left boundary into the left bulk.
left_boundary_flux = F.ParticleFluxBC(
subdomain=left_boundary,
species=H_left,
value=0.5,
)
# Interface coupling. J is the rate at which the interface *gains* particles from that
# side, so the interface gets +J as a source and the bulk gets -J as a flux.
def J_left(c_int, c_left):
return k1 * c_left * (1.0 - c_int / c_int_max) - k2 * c_int
def J_right(c_int, c_right):
return k3 * c_right * (1.0 - c_int / c_int_max) - k4 * c_int
interface_sources = [
F.ParticleSource(
value=J_left,
species=H_int,
volume=gamma,
species_dependent_value={"c_int": H_int, "c_left": H_left},
),
F.ParticleSource(
value=J_right,
species=H_int,
volume=gamma,
species_dependent_value={"c_int": H_int, "c_right": H_right},
),
]
interface_fluxes = [
F.ParticleFluxBC(
subdomain=gamma,
species=H_left,
value=lambda c_int, c_left: -J_left(c_int, c_left),
species_dependent_value={"c_int": H_int, "c_left": H_left},
),
F.ParticleFluxBC(
subdomain=gamma,
species=H_right,
value=lambda c_int, c_right: -J_right(c_int, c_right),
species_dependent_value={"c_int": H_int, "c_right": H_right},
),
]
right_dirichlet = F.FixedConcentrationBC(
subdomain=right_boundary,
value=0.0,
species=H_right,
)
model = F.HydrogenTransportProblemDiscontinuous(
mesh=F.Mesh(mesh),
species=[H_left, H_right, H_int],
subdomains=[left, right, gamma, left_boundary, right_boundary],
reactions=[decay],
sources=[*interface_sources],
boundary_conditions=[left_boundary_flux, *interface_fluxes, right_dirichlet],
temperature=500,
settings=F.Settings(
atol=1e-10,
rtol=1e-10,
transient=True,
final_time=T_final,
stepsize=dt,
),
exports=[
F.VTXSpeciesExport(filename="results/c_left.bp", field=H_left, subdomain=left),
F.VTXSpeciesExport(
filename="results/c_right.bp", field=H_right, subdomain=right
),
F.VTXSpeciesExport(filename="results/c_int.bp", field=H_int, subdomain=gamma),
],
)
model.initialise()
model.run()
c_left = H_left.subdomain_to_post_processing_solution[left].x.array
c_right = H_right.subdomain_to_post_processing_solution[right].x.array
c_int = H_int.subdomain_to_post_processing_solution[gamma].x.array
print("c_left range", c_left.min(), c_left.max())
print("c_right range", c_right.min(), c_right.max())
print("c_int range", c_int.min(), c_int.max())
Then, I commented out the new create_implicit_species_value_fenics function:
FESTIM/src/festim/hydrogen_transport_problem.py
Lines 1943 to 1988 in c9347c4
Finally, I run the new simulation. I don't get any error and the resulting values appear to make sense:
c_left range 0.9099848067293339 1.3627224928594979
c_right range 0.0 0.17317388922589677
c_int range 0.26298103965231406 0.26298103965231406Running with the new create_implicit_species_value_fenics gives an identical result.
Meanwhile, running the previous ImplicitSpecies example without the new create_implicit_species_value_fenics gives an UnboundLocalError (this is what I forgot to include in my previous comment):
Full Traceback
Traceback (most recent call last):
File "~/mwe_internal_interfaces_traps.py", line 172, in <module>
model.initialise()
~~~~~~~~~~~~~~~~^^
File "~/src/festim/hydrogen_transport_problem.py", line 1408, in initialise
self.create_formulation()
~~~~~~~~~~~~~~~~~~~~~~~^^
File "~/src/festim/hydrogen_transport_problem.py", line 2533, in create_formulation
dolfinx.fem.form(g, entity_maps=entity_maps, jit_options=jit_options)
~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 476, in form
return _create_form(form)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 472, in _create_form
return list(map(lambda sub_form: _create_form(sub_form), form))
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 472, in <lambda>
return list(map(lambda sub_form: _create_form(sub_form), form))
~~~~~~~~~~~~^^^^^^^^^^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 468, in _create_form
return _form(form)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 384, in _form
ufcx_form, module, code = jit.ffcx_jit(
~~~~~~~~~~~~^
comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/jit.py", line 61, in mpi_jit
return local_jit(*args, **kwargs)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/dolfinx/jit.py", line 221, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/codegeneration/jit.py", line 244, in compile_forms
raise e
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/codegeneration/jit.py", line 224, in compile_forms
impl = _compile_objects(
decl,
...<9 lines>...
visualise=visualise,
)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/codegeneration/jit.py", line 349, in _compile_objects
code, _ = ffcx.compiler.compile_ufl_objects(
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
ufl_objects, namespace=module_name, options=options, visualise=visualise
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/compiler.py", line 116, in compile_ufl_objects
ir = compute_ir(analysis, _object_names, _namespace, options, visualise)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/ir/representation.py", line 261, in compute_ir
_compute_integral_ir(
~~~~~~~~~~~~~~~~~~~~^
fd,
^^^
...<4 lines>...
visualise,
^^^^^^^^^^
)
^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/ir/representation.py", line 470, in _compute_integral_ir
integral_ir = compute_integral_ir(
itg_data.domain.ufl_cell(),
...<5 lines>...
visualise,
)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/ir/integral.py", line 430, in compute_integral_ir
) = _compute_integral_ir(
~~~~~~~~~~~~~~~~~~~~^
expression,
^^^^^^^^^^^
...<7 lines>...
p,
^^
)
^
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/ir/integral.py", line 178, in _compute_integral_ir
mt_table_reference = build_optimized_tables(
quadrature_rule,
...<8 lines>...
atol=p["table_atol"],
)
File "~/anaconda3/envs/festim-dev/lib/python3.14/site-packages/ffcx/ir/elementtables.py", line 569, in build_optimized_tables
tbl = clamp_table_small_numbers(t["array"], rtol=rtol, atol=atol)
^
UnboundLocalError: cannot access local variable 't' where it is not associated with a valueSo it seems like it is only an issue with the ImplicitSpecies object
Collaborator
Author
Contributor
There was a problem hiding this comment.
No problem! I can give it another lookover myself later tonight (ended up doing unrelated work this afternoon...)
| subdomain (F.VolumeSubdomain): a subdomain of the geometry | ||
| """ | ||
| is_manifold = subdomain.codim(self.mesh.vdim) == 1 | ||
| if is_manifold and self.mesh.coordinate_system != CoordinateSystem.CARTESIAN: |
Comment on lines
+2040
to
+2041
| dx = self.subdomain_measure(subdomain) | ||
| dx_grad = dx |
Collaborator
Author
There was a problem hiding this comment.
maybe we don't need a separate dx_grad
Collaborator
I can't reproduce this issue using DOLFINx v0.10.0-r1, v0.11.0 or the main branch. |
jorgensd
reviewed
Aug 5, 2026
| from abc import ABC, abstractmethod | ||
|
|
||
|
|
||
| def compute_ordered_interior_facet_data( |
Collaborator
There was a problem hiding this comment.
Since you rely on Scifem, you could consider reducing this function quite a bit: https://scientificcomputing.github.io/scifem/_modules/scifem/mesh.html#compute_interface_data
I think you could reproduce the logic with subdomain plus and - if you know what value they have (i.e. how they align with compute interface data, and then do the switch.
Collaborator
Author
@jorgensd this specific issue has been fixed on this branch so that may be why you can't reproduce. Maybe @ee-nn can provide a pure dolfinx example? Edit: Here's a pure dolfinx example (0.11) Details
"""MWE: a fem.Constant from the parent mesh in a form integrated over a
codim-1 submesh fails to compile with an UnboundLocalError inside FFCx.
"""
from mpi4py import MPI
import dolfinx
import numpy as np
import ufl
print(f"dolfinx {dolfinx.__version__}")
# --- minimal reproducer ----------------------------------------------------
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8)
tdim = mesh.topology.dim
# codim-1 submesh: the line y = 0.5
facets = dolfinx.mesh.locate_entities(mesh, tdim - 1, lambda x: np.isclose(x[1], 0.5))
submesh, cell_map, v_map, n_map = dolfinx.mesh.create_submesh(mesh, tdim - 1, facets)
V = dolfinx.fem.functionspace(submesh, ("Lagrange", 1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
dx = ufl.Measure("dx", domain=submesh)
c_parent = dolfinx.fem.Constant(mesh, 1.0) # <-- lives on the PARENT mesh
c_sub = dolfinx.fem.Constant(submesh, 1.0) # <-- lives on the submesh
dolfinx.fem.form(ufl.inner(c_sub * u, v) * dx) # works
try:
dolfinx.fem.form(ufl.inner(c_parent * u, v) * dx) # UnboundLocalError
except UnboundLocalError as e:
print(f"parent-mesh constant on a codim-1 submesh: {type(e).__name__}: {e}")
# --- the full matrix -------------------------------------------------------
def build(codim: int, constant_on: str, entity_maps: bool):
dim = tdim - codim
if codim == 1:
ents = dolfinx.mesh.locate_entities(mesh, dim, lambda x: np.isclose(x[1], 0.5))
else:
ents = dolfinx.mesh.locate_entities(mesh, dim, lambda x: x[1] <= 0.5 + 1e-12)
sub, cmap, _, _ = dolfinx.mesh.create_submesh(mesh, dim, ents)
W = dolfinx.fem.functionspace(sub, ("Lagrange", 1))
uu, vv = ufl.TrialFunction(W), ufl.TestFunction(W)
c = dolfinx.fem.Constant(mesh if constant_on == "parent" else sub, 1.0)
a = ufl.inner(c * uu, vv) * ufl.Measure("dx", domain=sub)
kwargs = {"entity_maps": [cmap]} if entity_maps else {}
dolfinx.fem.form(a, **kwargs)
if __name__ == "__main__":
for codim in (0, 1):
for constant_on in ("sub", "parent"):
for entity_maps in (False, True):
label = (
f"codim={codim} constant on {constant_on:6s} "
f"entity_maps={entity_maps}"
)
try:
build(codim, constant_on, entity_maps)
print(f"{label}: OK")
except Exception as e:
print(f"{label}: {type(e).__name__}: {e}") |
RemDelaporteMathurin
force-pushed
the
rem/codim
branch
from
August 5, 2026 12:42
0039f22 to
2852b89
Compare
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Description
Summary
Related Issues
Motivation and Context
Type of Change
Testing
pytest)Code Quality Checklist
ruff format .)ruff check .)Documentation
Breaking Changes
Screenshots/Examples
Additional Notes