From 82dbd1ffb77bec8b5389254e188e49d7892f48af Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 16:15:07 -0400 Subject: [PATCH 01/27] make sure post processing functions have a name --- src/festim/hydrogen_transport_problem.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 14d9c35b4..ee37837d2 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -646,6 +646,7 @@ def assign_functions_to_species(self): spe.sub_function_space = self.function_space.sub(idx) spe.sub_function = self.u.sub(idx) # TODO add this to discontinuous class spe.post_processing_solution = self.u.sub(idx).collapse() + spe.post_processing_solution.name = spe.name spe.collapsed_function_space, spe.map_sub_to_main_solution = ( self.function_space.sub(idx).collapse() ) From 14dcf9580d127503c90bb2ed2348c4829b403334 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 17:22:02 -0400 Subject: [PATCH 02/27] added test --- test/test_vtx.py | 81 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/test/test_vtx.py b/test/test_vtx.py index edbb2e531..335a9bda9 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -179,3 +179,84 @@ def test_filename_temp_raises_error_when_wrong_type(): """Test that the filename attribute for VTXTemperature export raises an error if the extension is not .bp""" with pytest.raises(TypeError): F.VTXTemperatureExport(1) + + +@pytest.mark.parametrize( + "expression", + [ + lambda x: x[0] + x[1] * 2, + lambda T: T + 1, + lambda c_A, c_B: c_A + c_B, + lambda c_A, c_B, x: c_A * c_B + x[0], + lambda c_A, T, x: c_A * T + x[0], + ], +) +def test_custom_field(tmp_path, expression): + """ + Test custom field export functionality. + This test checks that a custom field can be created with various types of + expressions. + """ + + my_model = F.HydrogenTransportProblem() + + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + + vol = F.VolumeSubdomain(id=1, material=mat) + + top = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1)) + bottom = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[1], 0)) + left = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) + right = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 1)) + + my_model.subdomains = [vol, top, bottom, left, right] + + dolfinx_mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) + my_model.mesh = F.Mesh(dolfinx_mesh) + + A = F.Species("A") + B = F.Species("B") + C = F.Species("C") + D = F.Species("D") + + my_model.species = [A, B, C, D] + + my_model.boundary_conditions = ( + [ + F.FixedConcentrationBC(species=A, subdomain=top, value=1), + F.FixedConcentrationBC(species=B, subdomain=left, value=1), + ] + + [ + F.FixedConcentrationBC(species=C, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + + [ + F.FixedConcentrationBC(species=D, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + ) + + my_model.reactions = [ + F.Reaction( + reactant=[A, B], product=[C], k_0=1, E_k=0, p_0=0, E_p=0, volume=vol + ), + F.Reaction(reactant=[C], product=[D], k_0=0.1, E_k=0, p_0=0, E_p=0, volume=vol), + ] + + my_model.temperature = 300 + + my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) + + custom_field = F.CustomField( + filename=tmp_path / "custom_field.bp", + expression=expression, + species_dependent_value={"c_A": A, "c_B": B}, + ) + + my_model.exports = [ + custom_field, + ] + + my_model.initialise() + + my_model.run() From 9f672f0925eb539a9ba783cc53f3ee4c9004e457 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 17:22:14 -0400 Subject: [PATCH 03/27] initial implementation --- src/festim/__init__.py | 7 ++- src/festim/exports/__init__.py | 3 +- src/festim/exports/vtx.py | 55 ++++++++++++++++++++++++ src/festim/helpers.py | 2 + src/festim/hydrogen_transport_problem.py | 21 +++++++++ 5 files changed, 86 insertions(+), 2 deletions(-) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 823d519d9..11db36aec 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -35,7 +35,12 @@ from .exports.total_surface import TotalSurface from .exports.total_volume import TotalVolume from .exports.volume_quantity import VolumeQuantity -from .exports.vtx import ExportBaseClass, VTXSpeciesExport, VTXTemperatureExport +from .exports.vtx import ( + ExportBaseClass, + VTXSpeciesExport, + VTXTemperatureExport, + CustomField, +) from .exports.xdmf import XDMFExport from .heat_transfer_problem import HeatTransferProblem from .helpers import ( diff --git a/src/festim/exports/__init__.py b/src/festim/exports/__init__.py index 31a287244..57c0e822f 100644 --- a/src/festim/exports/__init__.py +++ b/src/festim/exports/__init__.py @@ -10,12 +10,13 @@ from .total_surface import TotalSurface from .total_volume import TotalVolume from .volume_quantity import VolumeQuantity -from .vtx import ExportBaseClass, VTXSpeciesExport, VTXTemperatureExport +from .vtx import ExportBaseClass, VTXSpeciesExport, VTXTemperatureExport, CustomField from .xdmf import XDMFExport __all__ = [ "AverageSurface", "AverageVolume", + "CustomField", "ExportBaseClass", "MaximumSurface", "MaximumVolume", diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 7159360fd..cf2005d66 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -2,9 +2,11 @@ from pathlib import Path from dolfinx import fem, io +import ufl from festim.species import Species from festim.subdomain.volume_subdomain import VolumeSubdomain +import inspect class ExportBaseClass: @@ -172,3 +174,56 @@ def get_functions(self) -> list[fem.Function]: field.subdomain_to_post_processing_solution[self._subdomain] ) return outfiles + + +class CustomField(ExportBaseClass): + function: fem.Function + writer: io.VTXWriter + dolfinx_expression: fem.Expression + + def __init__( + self, + filename, + expression, + species_dependent_value=None, + times=None, + subdomain: VolumeSubdomain = None, + checkpoint=False, + ): + super().__init__( + filename=filename, + times=times, + ext=".bp", + ) + self.expression = expression + self.species_dependent_value = species_dependent_value or {} + self.checkpoint = checkpoint + self.subdomain = subdomain + + def set_dolfinx_expression( + self, + temperature: fem.Constant | fem.Function, + species: list, + time: fem.Constant, + ): + # check + arguments = inspect.signature(self.expression).parameters + kwargs = {} + if "t" in arguments: + kwargs["t"] = time + if "x" in arguments: + x = ufl.SpatialCoordinate(time._ufl_domain) + kwargs["x"] = x + if "T" in arguments: + kwargs["T"] = temperature + # check if there are other arguments and if they are in species_dependent_value + for arg in arguments: + if arg in self.species_dependent_value: + spe = self.species_dependent_value[arg] + kwargs[arg] = spe.post_processing_solution + # TODO handle discontinuous case and use self.subdomain + + self.dolfinx_expression = fem.Expression( + self.expression(**kwargs), + self.function.function_space.element.interpolation_points, + ) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 0a1ba5290..236d31159 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -34,6 +34,7 @@ def as_fenics_constant( ) +# TODO change this to accept species dependent values def as_mapped_function( value: Callable, function_space: fem.FunctionSpace | None = None, @@ -68,6 +69,7 @@ def as_mapped_function( return value(**kwargs) +# TODO change this to accept species dependent values def as_fenics_interp_expr_and_function( value: Callable, function_space: dolfinx.fem.function.FunctionSpace, diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index ee37837d2..f29b1f38f 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -449,6 +449,22 @@ def initialise_exports(self): else: adios4dolfinx.write_mesh(export.filename, mesh=self.mesh.mesh) + elif isinstance(export, exports.CustomField): + export.function = fem.Function(self.V_CG_1) + export.set_dolfinx_expression( + temperature=self.temperature_fenics, + species=self.species, + time=self.t, + ) + + export.writer = dolfinx.io.VTXWriter( + comm=export.function.function_space.mesh.comm, + filename=export.filename, + output=export.function, + engine="BP5", + ) + continue + elif isinstance(export, exports.SurfaceQuantity | exports.VolumeQuantity): # raise not implemented error if the derived quantity don't match the # type of mesh eg. SurfaceFlux is used with cylindrical mesh @@ -629,6 +645,7 @@ def define_function_spaces(self, element_degree: int = 1): ) self.V_DG_0 = fem.functionspace(self.mesh.mesh, element_DG0) self.V_DG_1 = fem.functionspace(self.mesh.mesh, element_DG1) + self.V_CG_1 = fem.functionspace(self.mesh.mesh, ("CG", 1)) self.u = fem.Function(self.function_space) self.u_n = fem.Function(self.function_space) @@ -982,6 +999,10 @@ def post_processing(self): self._get_temperature_field_as_function() ) export.writer.write(float(self.t)) + elif isinstance(export, exports.CustomField): + # update internal function + export.function.interpolate(export.dolfinx_expression) + export.writer.write(float(self.t)) # TODO if export type derived quantity if isinstance(export, exports.SurfaceQuantity): From 78e491a383882fb4f794b2c16de56da647b50c4c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 17:28:42 -0400 Subject: [PATCH 04/27] backwards compatibility --- src/festim/exports/vtx.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index cf2005d66..03d9d87c6 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -5,6 +5,7 @@ import ufl from festim.species import Species +from festim.helpers import get_interpolation_points from festim.subdomain.volume_subdomain import VolumeSubdomain import inspect @@ -225,5 +226,5 @@ def set_dolfinx_expression( self.dolfinx_expression = fem.Expression( self.expression(**kwargs), - self.function.function_space.element.interpolation_points, + get_interpolation_points(self.function.function_space.element), ) From d30f07cbef22bfd60b47ee52163fd39a269a4e0a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 21:00:24 -0400 Subject: [PATCH 05/27] support for mixed domain --- src/festim/exports/vtx.py | 22 +++++-- src/festim/hydrogen_transport_problem.py | 35 +++++++++-- test/test_vtx.py | 80 ++++++++++++++++++++++++ 3 files changed, 127 insertions(+), 10 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 03d9d87c6..1cbce66e4 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -216,14 +216,28 @@ def set_dolfinx_expression( x = ufl.SpatialCoordinate(time._ufl_domain) kwargs["x"] = x if "T" in arguments: - kwargs["T"] = temperature + if ( + isinstance(temperature, fem.Function) + and self.subdomain.sub_T is not None + ): + kwargs["T"] = ( + self.subdomain.sub_T + ) # NOTE I'm not sure that sub_T is updated at every time step + else: + kwargs["T"] = temperature # check if there are other arguments and if they are in species_dependent_value for arg in arguments: if arg in self.species_dependent_value: spe = self.species_dependent_value[arg] - kwargs[arg] = spe.post_processing_solution - # TODO handle discontinuous case and use self.subdomain - + if spe.subdomain_to_post_processing_solution: + kwargs[arg] = spe.subdomain_to_post_processing_solution[ + self.subdomain + ] + else: + kwargs[arg] = spe.post_processing_solution + assert kwargs[arg] is not None, ( + f"Argument {arg} not found in species_dependent_value" + ) self.dolfinx_expression = fem.Expression( self.expression(**kwargs), get_interpolation_points(self.function.function_space.element), diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index f29b1f38f..4a529315a 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1117,6 +1117,7 @@ def __init__( self.interfaces = interfaces or [] self.surface_to_volume = surface_to_volume or {} self.subdomain_to_species = {} # maps subdomain to species defined in it + self.subdomain_to_V_CG1 = {} @property def method_interface(self): @@ -1357,6 +1358,10 @@ def define_function_spaces( u = dolfinx.fem.Function(V) u_n = dolfinx.fem.Function(V) + self.subdomain_to_V_CG1[subdomain] = dolfinx.fem.functionspace( + subdomain.submesh, ("CG", 1) + ) + # store attributes in the subdomain object subdomain.u = u subdomain.u_n = u_n @@ -1684,7 +1689,22 @@ def initialise_exports(self): export.filename, mesh=functions[0].function_space.mesh, ) + elif isinstance(export, exports.CustomField): + # need to find an appropriate function space on the right submesh + V = self.subdomain_to_V_CG1[export.subdomain] + export.function = fem.Function(V) + export.set_dolfinx_expression( + temperature=self.temperature_fenics, # need to pass the right temperature + species=self.species, + time=self.t, + ) + export.writer = dolfinx.io.VTXWriter( + comm=export.function.function_space.mesh.comm, + filename=export.filename, + output=export.function, + engine="BP5", + ) # compute diffusivity function for surface fluxes # for the discontinuous case, we don't use D_global as in # HydrogenTransportProblem @@ -1726,11 +1746,11 @@ def post_processing(self): continue # handle VTX exports if isinstance(export, exports.ExportBaseClass): - if not isinstance(export, exports.VTXSpeciesExport): - raise NotImplementedError( - f"Export type {type(export)} not implemented" - ) - if isinstance(export, exports.VTXSpeciesExport): + if isinstance(export, exports.CustomField): + # update internal function + export.function.interpolate(export.dolfinx_expression) + export.writer.write(float(self.t)) + elif isinstance(export, exports.VTXSpeciesExport): if export._checkpoint: for species in export.field: post_processing_solution = ( @@ -1746,7 +1766,10 @@ def post_processing(self): ) else: export.writer.write(float(self.t)) - + else: + raise NotImplementedError( + f"Export type {type(export)} not implemented" + ) # handle derived quantities if isinstance(export, exports.SurfaceQuantity): if isinstance( diff --git a/test/test_vtx.py b/test/test_vtx.py index 335a9bda9..9daf49118 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -260,3 +260,83 @@ def test_custom_field(tmp_path, expression): my_model.initialise() my_model.run() + + +@pytest.mark.parametrize( + "expression", + [ + lambda x: x[0] + x[1] * 2, + lambda T: T + 1, + lambda c_A, c_B: c_A + c_B, + lambda c_A, T: c_A * T, + # FIXME: the following don't work (easily) because of https://github.com/FEniCS/dolfinx/issues/3207 + # lambda c_A, c_B, x: c_A * c_B + x[0], + # lambda c_A, T, x: c_A * T + x[0], + # lambda T, x: T + x[0], + ], +) +def test_custom_field_discontinuous(tmp_path, expression): + """ + Test custom field export functionality. + This test checks that a custom field can be created with various types of + expressions. + """ + + my_model = F.HydrogenTransportProblemDiscontinuous() + + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + + vol = F.VolumeSubdomain(id=1, material=mat) + + top = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1)) + bottom = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[1], 0)) + left = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) + right = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 1)) + + my_model.subdomains = [vol, top, bottom, left, right] + + my_model.surface_to_volume = {top: vol, bottom: vol, left: vol, right: vol} + + dolfinx_mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) + my_model.mesh = F.Mesh(dolfinx_mesh) + + A = F.Species("A", subdomains=my_model.volume_subdomains) + B = F.Species("B", subdomains=my_model.volume_subdomains) + C = F.Species("C", subdomains=my_model.volume_subdomains) + D = F.Species("D", subdomains=my_model.volume_subdomains) + + my_model.species = [A, B, C, D] + + my_model.boundary_conditions = ( + [ + F.FixedConcentrationBC(species=A, subdomain=top, value=1), + F.FixedConcentrationBC(species=B, subdomain=left, value=1), + ] + + [ + F.FixedConcentrationBC(species=C, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + + [ + F.FixedConcentrationBC(species=D, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + ) + + my_model.temperature = lambda x: 300 + 100 * x[0] + + my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) + + custom_field = F.CustomField( + filename=tmp_path / "custom_field.bp", + expression=expression, + species_dependent_value={"c_A": A, "c_B": B}, + subdomain=vol, + ) + + my_model.exports = [ + custom_field, + ] + + my_model.initialise() + + my_model.run() From 00193d3ece222b6a4ac852c0a4ec001d06c9ad21 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 10:52:15 -0400 Subject: [PATCH 06/27] docstrings --- src/festim/exports/vtx.py | 45 +++++++++++++++++++++++++++++++++------ 1 file changed, 38 insertions(+), 7 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 1cbce66e4..ef927b6de 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -1,13 +1,15 @@ +import inspect import warnings +from collections.abc import Callable from pathlib import Path +from typing import Union from dolfinx import fem, io import ufl -from festim.species import Species from festim.helpers import get_interpolation_points +from festim.species import Species from festim.subdomain.volume_subdomain import VolumeSubdomain -import inspect class ExportBaseClass: @@ -178,18 +180,47 @@ def get_functions(self) -> list[fem.Function]: class CustomField(ExportBaseClass): + """Export a custom field to a VTX file + + Args: + filename: The name of the output file + expression: A function evaluating the custom field. Positional + arguments of the function can be "t" (time), "x" (spatial coordinate), + "T" (temperature), or any key from the `species_dependent_value` dictionary. + species_dependent_value: A dictionary mapping argument names + in `expression` to Species objects. Defaults to None. + times: if provided, the field will be exported at these timesteps. Otherwise + exports at all timesteps. Defaults to None. + subdomain: The volume subdomain on which the custom + field is evaluated. Defaults to None. + checkpoint: If True, the export will be a checkpoint file using + adios4dolfinx and won't be readable by ParaView. Default is False. + + Attributes: + filename: The name of the output file + expression: A function evaluating the custom field. + species_dependent_value: A dictionary mapping argument names to Species objects. + subdomain: The volume subdomain on which the custom field is evaluated. + checkpoint: If True, the export will be a checkpoint file. + times: if provided, the field will be exported at these timesteps. Otherwise + exports at all timesteps. + function: the function containing the custom field values + writer: The VTXWriter object used to write the file + dolfinx_expression: the dolfinx expression used to evaluate the function + """ + function: fem.Function writer: io.VTXWriter dolfinx_expression: fem.Expression def __init__( self, - filename, - expression, - species_dependent_value=None, - times=None, + filename: Union[str, Path], + expression: Callable, + species_dependent_value: Union[dict[str, Species], None] = None, + times: Union[list[float], list[int], None] = None, subdomain: VolumeSubdomain = None, - checkpoint=False, + checkpoint: bool = False, ): super().__init__( filename=filename, From 598103380e86ef38c019789993e9e3d585228d4d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 10:53:09 -0400 Subject: [PATCH 07/27] ruff --- src/festim/exports/vtx.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index ef927b6de..9b74dc95e 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -4,8 +4,8 @@ from pathlib import Path from typing import Union -from dolfinx import fem, io import ufl +from dolfinx import fem, io from festim.helpers import get_interpolation_points from festim.species import Species From dc5c21f8e984e3059b133519ba028775bc2d1132 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 10:55:45 -0400 Subject: [PATCH 08/27] removed unused argument --- src/festim/exports/vtx.py | 14 +++++++++++++- src/festim/hydrogen_transport_problem.py | 2 -- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 9b74dc95e..a3baa29b7 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -212,6 +212,10 @@ class CustomField(ExportBaseClass): function: fem.Function writer: io.VTXWriter dolfinx_expression: fem.Expression + expression: Callable + species_dependent_value: dict[str, Species] + subdomain: VolumeSubdomain + checkpoint: bool def __init__( self, @@ -235,9 +239,17 @@ def __init__( def set_dolfinx_expression( self, temperature: fem.Constant | fem.Function, - species: list, time: fem.Constant, ): + """ + Set the dolfinx expression used to evaluate the custom field. This is done by + evaluating the user-provided expression with the appropriate arguments and using + the result to create a dolfinx expression. + + Args: + temperature: The temperature field to use in the expression + time: The time to use in the expression + """ # check arguments = inspect.signature(self.expression).parameters kwargs = {} diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 4a529315a..b065d8ded 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -453,7 +453,6 @@ def initialise_exports(self): export.function = fem.Function(self.V_CG_1) export.set_dolfinx_expression( temperature=self.temperature_fenics, - species=self.species, time=self.t, ) @@ -1695,7 +1694,6 @@ def initialise_exports(self): export.function = fem.Function(V) export.set_dolfinx_expression( temperature=self.temperature_fenics, # need to pass the right temperature - species=self.species, time=self.t, ) From d660bcac6cc8ba53341bb6a34b3ff0cedeef13a0 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 11:06:14 -0400 Subject: [PATCH 09/27] added comments --- src/festim/exports/vtx.py | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index a3baa29b7..20ac72767 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -250,8 +250,16 @@ def set_dolfinx_expression( temperature: The temperature field to use in the expression time: The time to use in the expression """ - # check + # check if we are in a mixed domain/discontinuous case + mixed_domain = any( + spe.subdomain_to_post_processing_solution + for spe in self.species_dependent_value.values() + ) or (self.subdomain.sub_T if self.subdomain else None) + + # get the arguments of the user-provided expression arguments = inspect.signature(self.expression).parameters + + # create a dictionary mapping the arguments to the appropriate values kwargs = {} if "t" in arguments: kwargs["t"] = time @@ -259,20 +267,18 @@ def set_dolfinx_expression( x = ufl.SpatialCoordinate(time._ufl_domain) kwargs["x"] = x if "T" in arguments: - if ( - isinstance(temperature, fem.Function) - and self.subdomain.sub_T is not None - ): - kwargs["T"] = ( - self.subdomain.sub_T - ) # NOTE I'm not sure that sub_T is updated at every time step + if isinstance(temperature, fem.Function) and mixed_domain: + # fem.Function in mixed domain/discontinuous case, use sub_T + # NOTE I'm not sure that sub_T is updated at every time step + kwargs["T"] = self.subdomain.sub_T else: + # else use the provided temperature kwargs["T"] = temperature # check if there are other arguments and if they are in species_dependent_value for arg in arguments: if arg in self.species_dependent_value: spe = self.species_dependent_value[arg] - if spe.subdomain_to_post_processing_solution: + if mixed_domain: kwargs[arg] = spe.subdomain_to_post_processing_solution[ self.subdomain ] @@ -281,6 +287,9 @@ def set_dolfinx_expression( assert kwargs[arg] is not None, ( f"Argument {arg} not found in species_dependent_value" ) + + # evaluate the user-provided expression with the appropriate arguments and create a + # dolfinx.fem.Expression self.dolfinx_expression = fem.Expression( self.expression(**kwargs), get_interpolation_points(self.function.function_space.element), From 6f647488e671755b4c9f3cd501416ec9be456483 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 11:25:50 -0400 Subject: [PATCH 10/27] x is the same domain as self.function --- src/festim/exports/vtx.py | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 20ac72767..4b8fa1f22 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -264,7 +264,7 @@ def set_dolfinx_expression( if "t" in arguments: kwargs["t"] = time if "x" in arguments: - x = ufl.SpatialCoordinate(time._ufl_domain) + x = ufl.SpatialCoordinate(self.function.function_space.mesh) kwargs["x"] = x if "T" in arguments: if isinstance(temperature, fem.Function) and mixed_domain: @@ -288,9 +288,34 @@ def set_dolfinx_expression( f"Argument {arg} not found in species_dependent_value" ) + self.check_valid_inputs(kwargs, mixed_domain) + # evaluate the user-provided expression with the appropriate arguments and create a # dolfinx.fem.Expression self.dolfinx_expression = fem.Expression( self.expression(**kwargs), get_interpolation_points(self.function.function_space.element), ) + + def check_valid_inputs(self, kwargs: dict, mixed_domain: bool): + """ + Check if we are in the mixed domain/discontinuous case and if the user-provided + expression is valid in this case. + dolfinx.fem.Expression does not support co-dim 0 submeshes and time is defined + on the parent mesh, so we cannot have time-dependent custom fields in the mixed + domain/discontinuous case. + + When https://github.com/FEniCS/dolfinx/issues/3207 is resolved we should be + able to support this + """ + + # check the domain of all kwargs and check that they are the same + + if mixed_domain and "t" in kwargs: + raise NotImplementedError( + "Time-dependent custom fields are not implemented in the case of a " + "mixed domain/discontinuous case." + "dolfinx.fem.Expression does not support co-dim 0 submeshes and time is" + "defined on the parent mesh." + "See https://github.com/FEniCS/dolfinx/issues/3207 for more details." + ) From 3c8c002ba4abc7712fe3df4d72f46b030541beab Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 11:26:04 -0400 Subject: [PATCH 11/27] more tests now pass + check for error --- test/test_vtx.py | 71 +++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 67 insertions(+), 4 deletions(-) diff --git a/test/test_vtx.py b/test/test_vtx.py index 9daf49118..d6452ad1a 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -269,10 +269,9 @@ def test_custom_field(tmp_path, expression): lambda T: T + 1, lambda c_A, c_B: c_A + c_B, lambda c_A, T: c_A * T, - # FIXME: the following don't work (easily) because of https://github.com/FEniCS/dolfinx/issues/3207 - # lambda c_A, c_B, x: c_A * c_B + x[0], - # lambda c_A, T, x: c_A * T + x[0], - # lambda T, x: T + x[0], + lambda c_A, c_B, x: c_A * c_B + x[0], + lambda c_A, T, x: c_A * T + x[0], + lambda T, x: T + x[0], ], ) def test_custom_field_discontinuous(tmp_path, expression): @@ -340,3 +339,67 @@ def test_custom_field_discontinuous(tmp_path, expression): my_model.initialise() my_model.run() + + +@pytest.mark.parametrize( + "expression", + [lambda c_A, t: c_A * t, lambda t: t], +) +def test_custom_field_not_implemented_error(expression): + my_model = F.HydrogenTransportProblemDiscontinuous() + + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + + vol = F.VolumeSubdomain(id=1, material=mat) + + top = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1)) + bottom = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[1], 0)) + left = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) + right = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 1)) + + my_model.subdomains = [vol, top, bottom, left, right] + + my_model.surface_to_volume = {top: vol, bottom: vol, left: vol, right: vol} + + dolfinx_mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) + my_model.mesh = F.Mesh(dolfinx_mesh) + + A = F.Species("A", subdomains=my_model.volume_subdomains) + B = F.Species("B", subdomains=my_model.volume_subdomains) + C = F.Species("C", subdomains=my_model.volume_subdomains) + D = F.Species("D", subdomains=my_model.volume_subdomains) + + my_model.species = [A, B, C, D] + + my_model.boundary_conditions = ( + [ + F.FixedConcentrationBC(species=A, subdomain=top, value=1), + F.FixedConcentrationBC(species=B, subdomain=left, value=1), + ] + + [ + F.FixedConcentrationBC(species=C, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + + [ + F.FixedConcentrationBC(species=D, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + ) + + my_model.temperature = lambda x: 300 + 100 * x[0] + + my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) + + custom_field = F.CustomField( + filename="custom_field.bp", + expression=expression, + species_dependent_value={"c_A": A, "c_B": B}, + subdomain=vol, + ) + + my_model.exports = [ + custom_field, + ] + + with pytest.raises(NotImplementedError): + my_model.initialise() From c6d16d748e58f4512a016082fa4868f17963333e Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 11:55:51 -0400 Subject: [PATCH 12/27] pseudo code --- src/festim/__init__.py | 1 + src/festim/exports/__init__.py | 9 ++++++++- src/festim/exports/vtx.py | 28 ++++++++++++++++++++++++++++ 3 files changed, 37 insertions(+), 1 deletion(-) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 11db36aec..a93434547 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -40,6 +40,7 @@ VTXSpeciesExport, VTXTemperatureExport, CustomField, + ReactionRate, ) from .exports.xdmf import XDMFExport from .heat_transfer_problem import HeatTransferProblem diff --git a/src/festim/exports/__init__.py b/src/festim/exports/__init__.py index 57c0e822f..628428ef7 100644 --- a/src/festim/exports/__init__.py +++ b/src/festim/exports/__init__.py @@ -10,7 +10,13 @@ from .total_surface import TotalSurface from .total_volume import TotalVolume from .volume_quantity import VolumeQuantity -from .vtx import ExportBaseClass, VTXSpeciesExport, VTXTemperatureExport, CustomField +from .vtx import ( + ExportBaseClass, + VTXSpeciesExport, + VTXTemperatureExport, + CustomField, + ReactionRate, +) from .xdmf import XDMFExport __all__ = [ @@ -23,6 +29,7 @@ "MinimumSurface", "MinimumVolume", "Profile1DExport", + "ReactionRate", "SurfaceFlux", "SurfaceQuantity", "TotalSurface", diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 4b8fa1f22..9f2871109 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -10,6 +10,7 @@ from festim.helpers import get_interpolation_points from festim.species import Species from festim.subdomain.volume_subdomain import VolumeSubdomain +from festim.reaction import Reaction class ExportBaseClass: @@ -319,3 +320,30 @@ def check_valid_inputs(self, kwargs: dict, mixed_domain: bool): "defined on the parent mesh." "See https://github.com/FEniCS/dolfinx/issues/3207 for more details." ) + + +class ReactionRate(CustomField): + def __init__( + self, + reaction: Reaction, + filename: str | Path, + times: list[float] | None = None, + subdomain: VolumeSubdomain | None = None, + checkpoint: bool = False, + ): + # for example + def expression(T, c_a, c_b, c_c): + return reaction.reaction_term( + T, reactant_concentrations=[c_a, c_b], product_concentrations=[c_c] + ) + + super().__init__( + filename=filename, + expression=expression, + species_dependent_value={ + spe.name: spe for spe in reaction.reactant + reaction.product + }, + times=times, + subdomain=subdomain, + checkpoint=checkpoint, + ) From 1ae4c444aa0ad55d0b119074f2a4b39c057dbc68 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 13:00:57 -0400 Subject: [PATCH 13/27] removed the need for `override_solution` --- src/festim/hydrogen_transport_problem.py | 24 ------------------ src/festim/reaction.py | 32 +++++++++++++++++++----- src/festim/species.py | 17 +++++++++++++ 3 files changed, 43 insertions(+), 30 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index b065d8ded..ede66eb55 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1467,9 +1467,6 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain): if reaction.volume != subdomain: continue - # TODO remove - # temporarily overide the solution to the one of the subdomain - self.override_solution_attributes(reaction) # reactant for reactant in reaction.reactant: if isinstance(reactant, _species.Species): @@ -1592,27 +1589,6 @@ def create_formulation(self): }, ) - def override_solution_attributes(self, reaction: _reaction.Reaction): - """ - Reaction.reaction_term() relies on the .solution attribute of the species - however, in the discontinuous class, this attribute doesn't really make sense - since there is one solution per subdomain. - Therefore we temporarily override the .solution attribute based on the reactants, - products, and `others` if there are implicit species - """ - list_of_species_to_override = reaction.reactant + reaction.product - - # check if we have implicit species: - for reactant in reaction.reactant: - if isinstance(reactant, _species.ImplicitSpecies): - for other_spe in reactant.others: - if other_spe not in list_of_species_to_override: - list_of_species_to_override.append(other_spe) - - for species in list_of_species_to_override: - if isinstance(species, _species.Species): - species.solution = species.subdomain_to_solution[reaction.volume] - def create_solver(self): if Version(dolfinx.__version__) == Version("0.9.0"): self.solver = BlockedNewtonSolver( diff --git a/src/festim/reaction.py b/src/festim/reaction.py index 336b18f80..fb24a56d4 100644 --- a/src/festim/reaction.py +++ b/src/festim/reaction.py @@ -134,6 +134,24 @@ def reaction_term( The reaction term to be used in a formulation. """ + # make sure products is a list + products = self.product if isinstance(self.product, list) else [self.product] + + # detect if mixed_domain + mixed_domain = any( + isinstance(reactant, _Species) and reactant.subdomain_to_solution != {} + for reactant in self.reactant + ) or any( + isinstance(product, _Species) and product.subdomain_to_solution != {} + for product in products + ) + + def get_concentration(species): + if mixed_domain: + return species.concentration_submesh(self.volume) + else: + return species.concentration + if self.product == []: if self.p_0 is not None: raise ValueError( @@ -155,8 +173,6 @@ def reaction_term( "E_p cannot be None when reaction products are present." ) - products = self.product if isinstance(self.product, list) else [self.product] - # reaction rates k = self.k_0 * exp(-self.E_k / (_k_B * temperature)) @@ -173,18 +189,22 @@ def reaction_term( assert len(reactant_concentrations) == len(reactants) for i, reactant in enumerate(reactants): if reactant_concentrations[i] is None: - reactant_concentrations[i] = reactant.concentration + reactant_concentrations[i] = get_concentration(reactant) else: - reactant_concentrations = [reactant.concentration for reactant in reactants] + reactant_concentrations = [ + get_concentration(reactant) for reactant in reactants + ] # if product_concentrations is provided, use these concentrations if product_concentrations is not None: assert len(product_concentrations) == len(products) for i, product in enumerate(products): if product_concentrations[i] is None: - product_concentrations[i] = product.concentration + product_concentrations[i] = get_concentration(product) else: - product_concentrations = [product.concentration for product in products] + product_concentrations = [ + get_concentration(product) for product in products + ] # multiply all concentrations to be used in the term product_of_reactants = reactant_concentrations[0] diff --git a/src/festim/species.py b/src/festim/species.py index 9b9aa412b..5cc11b728 100644 --- a/src/festim/species.py +++ b/src/festim/species.py @@ -111,6 +111,12 @@ def __str__(self) -> str: def concentration(self): return self.solution + def concentration_submesh(self, subdomain: _VolumeSubdomain): + assert subdomain in self.subdomains, ( + f"Species {self.name} has no solution on subdomain {subdomain}." + ) + return self.subdomain_to_solution[subdomain] + @property def legacy(self) -> bool: """ @@ -174,6 +180,17 @@ def concentration(self): ) return self.value_fenics - sum([other.solution for other in self.others]) + def concentration_submesh(self, subdomain: _VolumeSubdomain): + if len(self.others) > 0: + for other in self.others: + assert other.subdomain_to_solution[subdomain], ( + f"Cannot compute concentration of {self.name} because {other.name}" + + f" has no solution on subdomain {subdomain}." + ) + return self.value_fenics - sum( + [other.subdomain_to_solution[subdomain] for other in self.others] + ) + def create_value_fenics(self, mesh, t: fem.Constant): """Creates the value of the density as a fenics object and sets it to self.value_fenics. From 88c1477baa41c1deb3084f1e44129edd5372db72 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 13:34:19 -0400 Subject: [PATCH 14/27] added expression --- src/festim/exports/vtx.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 9f2871109..a8a041b8d 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -8,6 +8,7 @@ from dolfinx import fem, io from festim.helpers import get_interpolation_points +from festim import k_B as _k_B from festim.species import Species from festim.subdomain.volume_subdomain import VolumeSubdomain from festim.reaction import Reaction @@ -331,11 +332,14 @@ def __init__( subdomain: VolumeSubdomain | None = None, checkpoint: bool = False, ): - # for example - def expression(T, c_a, c_b, c_c): - return reaction.reaction_term( - T, reactant_concentrations=[c_a, c_b], product_concentrations=[c_c] - ) + def expression(T, reactants, products): + k = reaction.k_0 * ufl.exp(-reaction.E_k / (_k_B * T)) + if reaction.p_0 and reaction.E_p: + p = reaction.p_0 * ufl.exp(-reaction.E_p / (_k_B * T)) + elif reaction.p_0: + p = reaction.p_0 + + return k * ufl.product(reactants) - p * ufl.product(products) super().__init__( filename=filename, From a63c369a6066f9557f81b4042d36cba7a81eb607 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 13:35:38 -0400 Subject: [PATCH 15/27] added not implemented error for implicit species --- src/festim/exports/vtx.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 4b8fa1f22..d2572f3cb 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -8,7 +8,7 @@ from dolfinx import fem, io from festim.helpers import get_interpolation_points -from festim.species import Species +from festim.species import Species, ImplicitSpecies from festim.subdomain.volume_subdomain import VolumeSubdomain @@ -278,6 +278,11 @@ def set_dolfinx_expression( for arg in arguments: if arg in self.species_dependent_value: spe = self.species_dependent_value[arg] + if isinstance(spe, ImplicitSpecies): + raise NotImplementedError( + "Custom fields depending on implicit species are not" + "implemented yet." + ) if mixed_domain: kwargs[arg] = spe.subdomain_to_post_processing_solution[ self.subdomain From 80f58bee2df566f01f543b639b5d6a7d4984e18d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 13:59:32 -0400 Subject: [PATCH 16/27] working with a hack --- src/festim/exports/vtx.py | 41 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 38 insertions(+), 3 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 2b60e4101..5fb6da3fc 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -337,14 +337,49 @@ def __init__( subdomain: VolumeSubdomain | None = None, checkpoint: bool = False, ): - def expression(T, reactants, products): + + reactant_names = [reactant.name for reactant in reaction.reactant] + if isinstance(reaction.product, list): + product_names = [product.name for product in reaction.product] + else: + product_names = [reaction.product.name] + + def expression(T, **kwargs): + _reactant_names = [kwargs[name] for name in reactant_names] + _product_names = [kwargs[name] for name in product_names] k = reaction.k_0 * ufl.exp(-reaction.E_k / (_k_B * T)) if reaction.p_0 and reaction.E_p: p = reaction.p_0 * ufl.exp(-reaction.E_p / (_k_B * T)) elif reaction.p_0: p = reaction.p_0 - - return k * ufl.product(reactants) - p * ufl.product(products) + else: + p = 0.0 + + return k * ufl.product(_reactant_names) - p * ufl.product(_product_names) + + # generate the __signature__ of the expression function so that it has the + # correct arguments for the user-provided expression. + # This is needed as set_dolfinx_expression() checks + # the arguments of expression and it would otherwise look for kwargs + sig_params = [inspect.Parameter("T", inspect.Parameter.POSITIONAL_OR_KEYWORD)] + # Use dict.fromkeys to preserve order and remove duplicates + for name in dict.fromkeys(reactant_names + product_names): + sig_params.append( + inspect.Parameter(name, inspect.Parameter.POSITIONAL_OR_KEYWORD) + ) + expression.__signature__ = inspect.Signature(sig_params) + + assert inspect.signature(expression).parameters.keys() == { + "T", + *reactant_names, + *product_names, + }, ( + "The expression for the reaction rate is automatically generated based on the " + "reaction provided. The arguments of the expression must be T (temperature) and " + "the names of the reactants and products. The current expression has arguments " + f"{inspect.signature(expression).parameters.keys()} but should have arguments " + f"T and {reactant_names + product_names}." + ) super().__init__( filename=filename, From 5fb83bcc9e45bb9659115a69e3f92a67dc624e51 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 14:05:46 -0400 Subject: [PATCH 17/27] refactoring --- src/festim/exports/vtx.py | 40 ++++++++++++++++++++++++--------------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 5fb6da3fc..d3eb42732 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -357,10 +357,31 @@ def expression(T, **kwargs): return k * ufl.product(_reactant_names) - p * ufl.product(_product_names) - # generate the __signature__ of the expression function so that it has the - # correct arguments for the user-provided expression. - # This is needed as set_dolfinx_expression() checks - # the arguments of expression and it would otherwise look for kwargs + self.override_signature(expression, reactant_names, product_names) + + super().__init__( + filename=filename, + expression=expression, + species_dependent_value={ + spe.name: spe for spe in reaction.reactant + reaction.product + }, + times=times, + subdomain=subdomain, + checkpoint=checkpoint, + ) + + def override_signature( + self, expression: Callable, reactant_names: list[str], product_names: list[str] + ): + """ + Override the signature of the expression function. This is needed to ensure that + the expression has the correct arguments for set_dolfinx_expression(). + + Args: + expression: The user-provided expression for the reaction rate. The arguments + of the expression must be T (temperature) and the names of the reactants + and products. + """ sig_params = [inspect.Parameter("T", inspect.Parameter.POSITIONAL_OR_KEYWORD)] # Use dict.fromkeys to preserve order and remove duplicates for name in dict.fromkeys(reactant_names + product_names): @@ -380,14 +401,3 @@ def expression(T, **kwargs): f"{inspect.signature(expression).parameters.keys()} but should have arguments " f"T and {reactant_names + product_names}." ) - - super().__init__( - filename=filename, - expression=expression, - species_dependent_value={ - spe.name: spe for spe in reaction.reactant + reaction.product - }, - times=times, - subdomain=subdomain, - checkpoint=checkpoint, - ) From e803ab1edc18539988abdfcdcb08a7594ae167b6 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 14:56:50 -0400 Subject: [PATCH 18/27] added direction --- src/festim/exports/vtx.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index d3eb42732..37f4c4f72 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -336,6 +336,7 @@ def __init__( times: list[float] | None = None, subdomain: VolumeSubdomain | None = None, checkpoint: bool = False, + direction: str = "both", ): reactant_names = [reactant.name for reactant in reaction.reactant] @@ -355,7 +356,15 @@ def expression(T, **kwargs): else: p = 0.0 - return k * ufl.product(_reactant_names) - p * ufl.product(_product_names) + forward = k * ufl.product(_reactant_names) + backward = p * ufl.product(_product_names) + + if direction == "forward": + return forward + elif direction == "backward": + return backward + else: + return forward - backward self.override_signature(expression, reactant_names, product_names) From 78f648c10121b46afa1c0b35329696078c39a60c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 12:51:14 -0400 Subject: [PATCH 19/27] renamed export --- src/festim/__init__.py | 2 +- src/festim/exports/__init__.py | 9 +++++++-- src/festim/exports/vtx.py | 2 +- src/festim/hydrogen_transport_problem.py | 8 ++++---- test/test_vtx.py | 6 +++--- 5 files changed, 16 insertions(+), 11 deletions(-) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 84aec9c37..392c35621 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -39,7 +39,7 @@ ExportBaseClass, VTXSpeciesExport, VTXTemperatureExport, - CustomField, + CustomFieldExport, ) from .exports.xdmf import XDMFExport from .heat_transfer_problem import HeatTransferProblem diff --git a/src/festim/exports/__init__.py b/src/festim/exports/__init__.py index 57c0e822f..e18892474 100644 --- a/src/festim/exports/__init__.py +++ b/src/festim/exports/__init__.py @@ -10,13 +10,18 @@ from .total_surface import TotalSurface from .total_volume import TotalVolume from .volume_quantity import VolumeQuantity -from .vtx import ExportBaseClass, VTXSpeciesExport, VTXTemperatureExport, CustomField +from .vtx import ( + ExportBaseClass, + VTXSpeciesExport, + VTXTemperatureExport, + CustomFieldExport, +) from .xdmf import XDMFExport __all__ = [ "AverageSurface", "AverageVolume", - "CustomField", + "CustomFieldExport", "ExportBaseClass", "MaximumSurface", "MaximumVolume", diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index d2572f3cb..30293869e 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -179,7 +179,7 @@ def get_functions(self) -> list[fem.Function]: return outfiles -class CustomField(ExportBaseClass): +class CustomFieldExport(ExportBaseClass): """Export a custom field to a VTX file Args: diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 2115d0762..2491da42c 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -448,7 +448,7 @@ def initialise_exports(self): else: adios4dolfinx.write_mesh(export.filename, mesh=self.mesh.mesh) - elif isinstance(export, exports.CustomField): + elif isinstance(export, exports.CustomFieldExport): export.function = fem.Function(self.V_CG_1) export.set_dolfinx_expression( temperature=self.temperature_fenics, @@ -997,7 +997,7 @@ def post_processing(self): self._get_temperature_field_as_function() ) export.writer.write(float(self.t)) - elif isinstance(export, exports.CustomField): + elif isinstance(export, exports.CustomFieldExport): # update internal function export.function.interpolate(export.dolfinx_expression) export.writer.write(float(self.t)) @@ -1707,7 +1707,7 @@ def initialise_exports(self): export.filename, mesh=functions[0].function_space.mesh, ) - elif isinstance(export, exports.CustomField): + elif isinstance(export, exports.CustomFieldExport): # need to find an appropriate function space on the right submesh V = self.subdomain_to_V_CG1[export.subdomain] export.function = fem.Function(V) @@ -1763,7 +1763,7 @@ def post_processing(self): continue # handle VTX exports if isinstance(export, exports.ExportBaseClass): - if isinstance(export, exports.CustomField): + if isinstance(export, exports.CustomFieldExport): # update internal function export.function.interpolate(export.dolfinx_expression) export.writer.write(float(self.t)) diff --git a/test/test_vtx.py b/test/test_vtx.py index 737bbdeac..da7923f3f 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -246,7 +246,7 @@ def test_custom_field(tmp_path, expression): my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) - custom_field = F.CustomField( + custom_field = F.CustomFieldExport( filename=tmp_path / "custom_field.bp", expression=expression, species_dependent_value={"c_A": A, "c_B": B}, @@ -324,7 +324,7 @@ def test_custom_field_discontinuous(tmp_path, expression): my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) - custom_field = F.CustomField( + custom_field = F.CustomFieldExport( filename=tmp_path / "custom_field.bp", expression=expression, species_dependent_value={"c_A": A, "c_B": B}, @@ -389,7 +389,7 @@ def test_custom_field_not_implemented_error(expression): my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) - custom_field = F.CustomField( + custom_field = F.CustomFieldExport( filename="custom_field.bp", expression=expression, species_dependent_value={"c_A": A, "c_B": B}, From 1eeaaca3f79b247bf9b11cb9a8a0e68bcbceb700 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 12:54:54 -0400 Subject: [PATCH 20/27] removed unused tests --- test/test_reaction.py | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/test/test_reaction.py b/test/test_reaction.py index 397d82faa..6481b874d 100644 --- a/test/test_reaction.py +++ b/test/test_reaction.py @@ -426,22 +426,3 @@ def test_product_setter_raise_error_E_p_no_product(): for subdomain in my_model.volume_subdomains: my_model.define_function_spaces(subdomain) - - -@pytest.mark.parametrize("reaction", [reac1, reac2]) -def test_override_solution_attributes(reaction): - """ - Tests the HydrogenTransportProblemDiscontinuous.override_solution_attributes method - Checks that the .solution attribute is the expected one based on the volume of the - reaction - """ - - # RUN - my_model.override_solution_attributes(reaction) - - # TEST - relevant_species = reaction.reactant + reaction.product + empty_traps.others - for species in relevant_species: - if isinstance(species, F.Species): - expected_solution = species.subdomain_to_solution[reaction.volume] - assert species.solution == expected_solution From 505209eb0f4f10303a1798b794e4e605cbf343f6 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 13:14:03 -0400 Subject: [PATCH 21/27] added mixed_domain property --- src/festim/exports/vtx.py | 32 +++++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 80b4f9d2f..b2d398922 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -238,6 +238,25 @@ def __init__( self.checkpoint = checkpoint self.subdomain = subdomain + @property + def mixed_domain(self) -> bool: + """ + Check if we are in a mixed domain/discontinuous case. This is the case if at least + one of the species in species_dependent_value is defined on a subdomain or if the + custom field is defined on a subdomain. + + Returns: + True if we are in a mixed domain/discontinuous case, False otherwise. + """ + all_explicit_species = [ + spe + for spe in self.species_dependent_value.values() + if isinstance(spe, Species) + ] + return any( + spe.subdomain_to_post_processing_solution for spe in all_explicit_species + ) or (self.subdomain.sub_T if self.subdomain else None) + def set_dolfinx_expression( self, temperature: fem.Constant | fem.Function, @@ -252,12 +271,6 @@ def set_dolfinx_expression( temperature: The temperature field to use in the expression time: The time to use in the expression """ - # check if we are in a mixed domain/discontinuous case - mixed_domain = any( - spe.subdomain_to_post_processing_solution - for spe in self.species_dependent_value.values() - ) or (self.subdomain.sub_T if self.subdomain else None) - # get the arguments of the user-provided expression arguments = inspect.signature(self.expression).parameters @@ -269,7 +282,7 @@ def set_dolfinx_expression( x = ufl.SpatialCoordinate(self.function.function_space.mesh) kwargs["x"] = x if "T" in arguments: - if isinstance(temperature, fem.Function) and mixed_domain: + if isinstance(temperature, fem.Function) and self.mixed_domain: # fem.Function in mixed domain/discontinuous case, use sub_T # NOTE I'm not sure that sub_T is updated at every time step kwargs["T"] = self.subdomain.sub_T @@ -285,7 +298,7 @@ def set_dolfinx_expression( "Custom fields depending on implicit species are not" "implemented yet." ) - if mixed_domain: + if self.mixed_domain: kwargs[arg] = spe.subdomain_to_post_processing_solution[ self.subdomain ] @@ -295,7 +308,7 @@ def set_dolfinx_expression( f"Argument {arg} not found in species_dependent_value" ) - self.check_valid_inputs(kwargs, mixed_domain) + self.check_valid_inputs(kwargs, self.mixed_domain) # evaluate the user-provided expression with the appropriate arguments and create a # dolfinx.fem.Expression @@ -391,6 +404,7 @@ def override_signature( of the expression must be T (temperature) and the names of the reactants and products. """ + breakpoint() sig_params = [inspect.Parameter("T", inspect.Parameter.POSITIONAL_OR_KEYWORD)] # Use dict.fromkeys to preserve order and remove duplicates for name in dict.fromkeys(reactant_names + product_names): From a2e3130063b1050c199c2b18dfb16fae0f769ee4 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 13:14:18 -0400 Subject: [PATCH 22/27] no breakpoint --- src/festim/exports/vtx.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index b2d398922..b53ddbb07 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -404,7 +404,6 @@ def override_signature( of the expression must be T (temperature) and the names of the reactants and products. """ - breakpoint() sig_params = [inspect.Parameter("T", inspect.Parameter.POSITIONAL_OR_KEYWORD)] # Use dict.fromkeys to preserve order and remove duplicates for name in dict.fromkeys(reactant_names + product_names): From f9d736220f09253672c3276c1146f81e249a8bc6 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 13:21:41 -0400 Subject: [PATCH 23/27] work for implicit species --- src/festim/exports/vtx.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index b53ddbb07..3e76b01bc 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -294,16 +294,17 @@ def set_dolfinx_expression( if arg in self.species_dependent_value: spe = self.species_dependent_value[arg] if isinstance(spe, ImplicitSpecies): - raise NotImplementedError( - "Custom fields depending on implicit species are not" - "implemented yet." - ) - if self.mixed_domain: - kwargs[arg] = spe.subdomain_to_post_processing_solution[ - self.subdomain - ] + if self.mixed_domain: + kwargs[arg] = spe.concentration_submesh(self.subdomain) + else: + kwargs[arg] = spe.concentration else: - kwargs[arg] = spe.post_processing_solution + if self.mixed_domain: + kwargs[arg] = spe.subdomain_to_post_processing_solution[ + self.subdomain + ] + else: + kwargs[arg] = spe.post_processing_solution assert kwargs[arg] is not None, ( f"Argument {arg} not found in species_dependent_value" ) From ed5a7511f7e3074206200bdfcb6884fc719e0a22 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 13:25:27 -0400 Subject: [PATCH 24/27] direction first --- src/festim/exports/vtx.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 3e76b01bc..ea4c6186c 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -347,10 +347,10 @@ def __init__( self, reaction: Reaction, filename: str | Path, + direction: str = "both", times: list[float] | None = None, subdomain: VolumeSubdomain | None = None, checkpoint: bool = False, - direction: str = "both", ): reactant_names = [reactant.name for reactant in reaction.reactant] From 5390a106d8383e9fcd64c8c12c13d199fbc9a6d0 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 13:32:10 -0400 Subject: [PATCH 25/27] refactoring --- src/festim/exports/vtx.py | 39 +++++++++++++++++++++------------------ 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index ea4c6186c..10c543b0f 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -289,36 +289,39 @@ def set_dolfinx_expression( else: # else use the provided temperature kwargs["T"] = temperature + # check if there are other arguments and if they are in species_dependent_value for arg in arguments: if arg in self.species_dependent_value: - spe = self.species_dependent_value[arg] - if isinstance(spe, ImplicitSpecies): - if self.mixed_domain: - kwargs[arg] = spe.concentration_submesh(self.subdomain) - else: - kwargs[arg] = spe.concentration - else: - if self.mixed_domain: - kwargs[arg] = spe.subdomain_to_post_processing_solution[ - self.subdomain - ] - else: - kwargs[arg] = spe.post_processing_solution + kwargs[arg] = self._get_species_function( + self.species_dependent_value[arg] + ) assert kwargs[arg] is not None, ( f"Argument {arg} not found in species_dependent_value" ) - self.check_valid_inputs(kwargs, self.mixed_domain) + self.check_valid_inputs(kwargs) - # evaluate the user-provided expression with the appropriate arguments and create a - # dolfinx.fem.Expression + # evaluate the user-provided expression with the appropriate arguments and + # create a dolfinx.fem.Expression self.dolfinx_expression = fem.Expression( self.expression(**kwargs), get_interpolation_points(self.function.function_space.element), ) - def check_valid_inputs(self, kwargs: dict, mixed_domain: bool): + def _get_species_function(self, spe: Species): + if isinstance(spe, ImplicitSpecies): + if self.mixed_domain: + return spe.concentration_submesh(self.subdomain) + else: + return spe.concentration + else: + if self.mixed_domain: + return spe.subdomain_to_post_processing_solution[self.subdomain] + else: + return spe.post_processing_solution + + def check_valid_inputs(self, kwargs: dict): """ Check if we are in the mixed domain/discontinuous case and if the user-provided expression is valid in this case. @@ -332,7 +335,7 @@ def check_valid_inputs(self, kwargs: dict, mixed_domain: bool): # check the domain of all kwargs and check that they are the same - if mixed_domain and "t" in kwargs: + if self.mixed_domain and "t" in kwargs: raise NotImplementedError( "Time-dependent custom fields are not implemented in the case of a " "mixed domain/discontinuous case." From 7c42b2d132c424b65577203103e59183e238f1fe Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 13:39:18 -0400 Subject: [PATCH 26/27] tests + make sure reaction products is a list --- src/festim/exports/vtx.py | 4 +- test/test_vtx.py | 107 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 110 insertions(+), 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 10c543b0f..5b4ddbbac 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -385,11 +385,13 @@ def expression(T, **kwargs): self.override_signature(expression, reactant_names, product_names) + reaction_products = reaction.product if isinstance(reaction.product, list) else [reaction.product] + super().__init__( filename=filename, expression=expression, species_dependent_value={ - spe.name: spe for spe in reaction.reactant + reaction.product + spe.name: spe for spe in reaction.reactant + reaction_products }, times=times, subdomain=subdomain, diff --git a/test/test_vtx.py b/test/test_vtx.py index da7923f3f..cf5cd3373 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -402,3 +402,110 @@ def test_custom_field_not_implemented_error(expression): with pytest.raises(NotImplementedError): my_model.initialise() + + +@pytest.mark.parametrize("direction", ["both", "forward", "backward"]) +@pytest.mark.parametrize("product_type", ["list", "single"]) +@pytest.mark.parametrize("p_0, E_p", [(0.01, 0.05), (0.01, 0.0), (0.0, 0.0)]) +def test_reaction_rate_export(tmp_path, direction, product_type, p_0, E_p): + """ + Test ReactionRate export functionality for different directions, product formats, + and reaction configurations. + """ + if p_0 == 0.0 and direction == "backward": + pytest.skip( + "Backward direction export not supported when backward reaction is disabled" + ) + my_model = F.HydrogenTransportProblem() + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + vol = F.VolumeSubdomain(id=1, material=mat) + top = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1)) + bottom = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[1], 0)) + left = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) + right = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 1)) + + my_model.subdomains = [vol, top, bottom, left, right] + + dolfinx_mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) + my_model.mesh = F.Mesh(dolfinx_mesh) + + A = F.Species("A") + B = F.Species("B") + C = F.Species("C") + + my_model.species = [A, B, C] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC(species=A, subdomain=top, value=1), + F.FixedConcentrationBC(species=B, subdomain=left, value=1), + F.FixedConcentrationBC(species=C, subdomain=bottom, value=0), + ] + + reaction = F.Reaction( + reactant=[A, B], + product=[C] if product_type == "list" else C, + k_0=1, + E_k=0.1, + p_0=p_0, + E_p=E_p, + volume=vol, + ) + + my_model.reactions = [reaction] + + my_model.temperature = 300 + + my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) + + reaction_rate_export = F.ReactionRate( + filename=tmp_path / f"reaction_rate_{direction}.bp", + reaction=reaction, + direction=direction, + ) + + my_model.exports = [reaction_rate_export] + + my_model.initialise() + my_model.run() + + +def test_reaction_rate_override_signature(): + """ + Test that ReactionRate signature override correctly updates signatures. + """ + mat = F.Material(D_0=1, E_D=0) + vol = F.VolumeSubdomain(id=1, material=mat) + A = F.Species("A") + B = F.Species("B") + reaction = F.Reaction( + reactant=[A], product=[B], k_0=1, E_k=0, p_0=0, E_p=0, volume=vol + ) + + rr = F.ReactionRate(reaction=reaction, filename="dummy.bp") + + def my_expression(x, y): + return x + y + + rr.override_signature(my_expression, ["A"], ["B"]) + import inspect + + sig = inspect.signature(my_expression) + assert set(sig.parameters.keys()) == {"T", "A", "B"} + + +def test_export_base_class_times_and_extension(tmp_path): + """ + Test that ExportBaseClass sorts times and warns when wrong extension is given. + """ + with pytest.warns(UserWarning, match="does not have .bp extension"): + export = F.ExportBaseClass( + filename=tmp_path / "wrong_extension.txt", ext=".bp", times=[3.0, 1.0, 2.0] + ) + + assert export.filename.suffix == ".bp" + assert export.times == [1.0, 2.0, 3.0] + + +def test_export_base_class_no_times(tmp_path): + export = F.ExportBaseClass(filename=tmp_path / "correct.bp", ext=".bp", times=None) + assert export.times is None From 7c45fe833a7daf655c14841f3988ecfbc16a2b54 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 13:40:26 -0400 Subject: [PATCH 27/27] test is closer to actual use case --- src/festim/exports/vtx.py | 6 +++++- test/test_vtx.py | 4 ++-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 5b4ddbbac..ecd262f5f 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -385,7 +385,11 @@ def expression(T, **kwargs): self.override_signature(expression, reactant_names, product_names) - reaction_products = reaction.product if isinstance(reaction.product, list) else [reaction.product] + reaction_products = ( + reaction.product + if isinstance(reaction.product, list) + else [reaction.product] + ) super().__init__( filename=filename, diff --git a/test/test_vtx.py b/test/test_vtx.py index cf5cd3373..554788d20 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -483,8 +483,8 @@ def test_reaction_rate_override_signature(): rr = F.ReactionRate(reaction=reaction, filename="dummy.bp") - def my_expression(x, y): - return x + y + def my_expression(**kwargs): + return kwargs.get("x", 0) + kwargs.get("y", 0) rr.override_signature(my_expression, ["A"], ["B"]) import inspect