From a08bd81626a50f2a877dbab1b53761d82e243566 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 13:46:57 -0400 Subject: [PATCH 01/15] added custom quantity --- src/festim/__init__.py | 2 + src/festim/exports/__init__.py | 4 ++ src/festim/exports/custom_quantity.py | 60 ++++++++++++++++++++++++++ src/festim/exports/derived_quantity.py | 54 +++++++++++++++++++++++ src/festim/exports/surface_quantity.py | 49 ++++----------------- src/festim/exports/volume_quantity.py | 50 +++++---------------- 6 files changed, 139 insertions(+), 80 deletions(-) create mode 100644 src/festim/exports/custom_quantity.py create mode 100644 src/festim/exports/derived_quantity.py diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 0be6a96c6..b2095f861 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -26,6 +26,8 @@ ) from .exports.average_surface import AverageSurface from .exports.average_volume import AverageVolume +from .exports.custom_quantity import CustomQuantity +from .exports.derived_quantity import DerivedQuantity from .exports.maximum_surface import MaximumSurface from .exports.maximum_volume import MaximumVolume from .exports.minimum_surface import MinimumSurface diff --git a/src/festim/exports/__init__.py b/src/festim/exports/__init__.py index 28c309bca..ad3988734 100644 --- a/src/festim/exports/__init__.py +++ b/src/festim/exports/__init__.py @@ -1,5 +1,7 @@ from .average_surface import AverageSurface from .average_volume import AverageVolume +from .custom_quantity import CustomQuantity +from .derived_quantity import DerivedQuantity from .maximum_surface import MaximumSurface from .maximum_volume import MaximumVolume from .minimum_surface import MinimumSurface @@ -23,6 +25,8 @@ "AverageSurface", "AverageVolume", "CustomFieldExport", + "CustomQuantity", + "DerivedQuantity", "ExportBaseClass", "MaximumSurface", "MaximumVolume", diff --git a/src/festim/exports/custom_quantity.py b/src/festim/exports/custom_quantity.py new file mode 100644 index 000000000..8a337bbe6 --- /dev/null +++ b/src/festim/exports/custom_quantity.py @@ -0,0 +1,60 @@ +from typing import Callable +import ufl +from scifem import assemble_scalar +from dolfinx import fem + +from festim.exports.derived_quantity import DerivedQuantity +from festim.subdomain.surface_subdomain import SurfaceSubdomain +from festim.subdomain.volume_subdomain import VolumeSubdomain + + +class CustomQuantity(DerivedQuantity): + """Export CustomQuantity. + + Args: + expr: function that returns a UFL expression + subdomain: subdomain on which the quantity is evaluated + title: title of the exported quantity + filename: name of the file to which the quantity is exported + + Attributes: + expr: function that returns a UFL expression + subdomain: subdomain on which the quantity is evaluated + title: title of the exported quantity + filename: name of the file to which the quantity is exported + t: list of time values + data: list of values of the quantity + """ + + def __init__( + self, + expr: Callable, + subdomain: SurfaceSubdomain | VolumeSubdomain | int, + title: str = "Custom Quantity", + filename: str | None = None, + ) -> None: + super().__init__(filename=filename) + self.expr = expr + self.subdomain = subdomain + self._title = title + self.ufl_expr = None + + @property + def title(self): + return self._title + + def compute(self, measure: ufl.Measure, entity_maps: dict | None = None): + """Computes the value of the custom quantity and appends it to the data list. + + Args: + measure: volume or surface measure of the model + entity_maps: entity maps relating parent mesh and submesh + """ + if self.ufl_expr is None: + raise ValueError("The UFL expression has not been evaluated yet.") + + form = fem.form( + self.ufl_expr * measure(self.subdomain.id), entity_maps=entity_maps + ) + self.value = assemble_scalar(form) + self.data.append(self.value) diff --git a/src/festim/exports/derived_quantity.py b/src/festim/exports/derived_quantity.py new file mode 100644 index 000000000..475fe51e7 --- /dev/null +++ b/src/festim/exports/derived_quantity.py @@ -0,0 +1,54 @@ +import csv + + +class DerivedQuantity: + """Base class for all derived quantities. + + Attributes: + filename: name of the file to which the quantity is exported + t: list of time values + data: list of values of the quantity + """ + + filename: str | None + t: list[float] + data: list[float] + + def __init__(self, filename: str | None = None) -> None: + self.filename = filename + self.t = [] + self.data = [] + self._first_time_export = True + + @property + def title(self): + raise NotImplementedError("title must be implemented by subclasses") + + @property + def filename(self): + return self._filename + + @filename.setter + def filename(self, value): + if value is None: + self._filename = None + elif not isinstance(value, str): + raise TypeError("filename must be of type str") + elif not value.endswith(".csv") and not value.endswith(".txt"): + raise ValueError("filename must end with .csv or .txt") + self._filename = value + + def write(self, t): + """If the filename doesnt exist yet, create it and write the header, then append + the time and value to the file.""" + + if self.filename is not None: + if self._first_time_export: + header = ["t(s)", f"{self.title}"] + with open(self.filename, mode="w+", newline="") as file: + writer = csv.writer(file) + writer.writerow(header) + self._first_time_export = False + with open(self.filename, mode="a", newline="") as file: + writer = csv.writer(file) + writer.writerow([t, self.value]) diff --git a/src/festim/exports/surface_quantity.py b/src/festim/exports/surface_quantity.py index 515ab368e..e915aebd6 100644 --- a/src/festim/exports/surface_quantity.py +++ b/src/festim/exports/surface_quantity.py @@ -1,10 +1,9 @@ -import csv - +from festim.exports.derived_quantity import DerivedQuantity from festim.species import Species from festim.subdomain.surface_subdomain import SurfaceSubdomain -class SurfaceQuantity: +class SurfaceQuantity(DerivedQuantity): """Export SurfaceQuantity. Args: @@ -28,29 +27,14 @@ class SurfaceQuantity: data: list[float] def __init__( - self, field: Species, surface: SurfaceSubdomain, filename: str | None = None + self, + field: Species | str, + surface: SurfaceSubdomain | int, + filename: str | None = None, ) -> None: + super().__init__(filename=filename) self.field = field self.surface = surface - self.filename = filename - - self.t = [] - self.data = [] - self._first_time_export = True - - @property - def filename(self): - return self._filename - - @filename.setter - def filename(self, value): - if value is None: - self._filename = None - elif not isinstance(value, str): - raise TypeError("filename must be of type str") - elif not value.endswith(".csv") and not value.endswith(".txt"): - raise ValueError("filename must end with .csv or .txt") - self._filename = value @property def surface(self): @@ -60,6 +44,7 @@ def surface(self): def surface(self, value): if not isinstance(value, int | SurfaceSubdomain) or isinstance(value, bool): raise TypeError("surface should be an int or F.SurfaceSubdomain") + self._surface = value @property @@ -68,23 +53,7 @@ def field(self): @field.setter def field(self, value): - # check that field is festim.Species if not isinstance(value, Species | str): - raise TypeError("field must be of type festim.Species") + raise TypeError("field must be of type F.Species or str") self._field = value - - def write(self, t): - """If the filename doesnt exist yet, create it and write the header, then append - the time and value to the file.""" - - if self.filename is not None: - if self._first_time_export: - header = ["t(s)", f"{self.title}"] - with open(self.filename, mode="w+", newline="") as file: - writer = csv.writer(file) - writer.writerow(header) - self._first_time_export = False - with open(self.filename, mode="a", newline="") as file: - writer = csv.writer(file) - writer.writerow([t, self.value]) diff --git a/src/festim/exports/volume_quantity.py b/src/festim/exports/volume_quantity.py index 1daaaac6a..a44f64606 100644 --- a/src/festim/exports/volume_quantity.py +++ b/src/festim/exports/volume_quantity.py @@ -1,10 +1,9 @@ -import csv - from festim.species import Species from festim.subdomain.volume_subdomain import VolumeSubdomain +from festim.exports.derived_quantity import DerivedQuantity -class VolumeQuantity: +class VolumeQuantity(DerivedQuantity): """Export VolumeQuantity. Args: @@ -27,28 +26,15 @@ class VolumeQuantity: t: list[float] data: list[float] - def __init__(self, field, volume, filename: str | None = None) -> None: + def __init__( + self, + field: Species | str, + volume: VolumeSubdomain | int, + filename: str | None = None, + ) -> None: + super().__init__(filename=filename) self.field = field self.volume = volume - self.filename = filename - - self.t = [] - self.data = [] - self._first_time_export = True - - @property - def filename(self): - return self._filename - - @filename.setter - def filename(self, value): - if value is None: - self._filename = None - elif not isinstance(value, str): - raise TypeError("filename must be of type str") - elif not value.endswith(".csv") and not value.endswith(".txt"): - raise ValueError("filename must end with .csv or .txt") - self._filename = value @property def field(self): @@ -56,23 +42,7 @@ def field(self): @field.setter def field(self, value): - # check that field is festim.Species if not isinstance(value, Species | str): - raise TypeError("field must be of type festim.Species") + raise TypeError("field must be of type F.Species or str") self._field = value - - def write(self, t): - """If the filename doesnt exist yet, create it and write the header, then append - the time and value to the file.""" - - if self.filename is not None: - if self._first_time_export: - header = ["t(s)", f"{self.title}"] - with open(self.filename, mode="w+", newline="") as file: - writer = csv.writer(file) - writer.writerow(header) - self._first_time_export = False - with open(self.filename, mode="a", newline="") as file: - writer = csv.writer(file) - writer.writerow([t, self.value]) From 1adf3f2ae842a37d1c375e1c8a22c82e5ca446b3 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 13:54:55 -0400 Subject: [PATCH 02/15] integrated with hydrogen transport problem --- src/festim/hydrogen_transport_problem.py | 111 ++++++++++++++++++++--- 1 file changed, 98 insertions(+), 13 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 267754ae6..d341ea0c5 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -462,7 +462,12 @@ def initialise_exports(self): ) continue - elif isinstance(export, exports.SurfaceQuantity | exports.VolumeQuantity): + elif isinstance( + export, + exports.SurfaceQuantity + | exports.VolumeQuantity + | exports.CustomQuantity, + ): # raise not implemented error if the derived quantity don't match the # type of mesh eg. SurfaceFlux is used with cylindrical mesh if self.mesh.coordinate_system != CoordinateSystem.CARTESIAN: @@ -472,16 +477,17 @@ def initialise_exports(self): ) # if name of species is given then replace with species object - if isinstance(export.field, list): - for idx, field in enumerate(export.field): - if isinstance(field, str): - export.field[idx] = _species.find_species_from_name( - field, self.species - ) - elif isinstance(export.field, str): - export.field = _species.find_species_from_name( - export.field, self.species - ) + if hasattr(export, "field"): + if isinstance(export.field, list): + for idx, field in enumerate(export.field): + if isinstance(field, str): + export.field[idx] = _species.find_species_from_name( + field, self.species + ) + elif isinstance(export.field, str): + export.field = _species.find_species_from_name( + export.field, self.species + ) # Initialize XDMFFile for writer if isinstance(export, exports.XDMFExport): @@ -521,10 +527,38 @@ def initialise_exports(self): export.volume_meshtags = self.volume_meshtags export.facet_meshtags = self.facet_meshtags # reset the data and time for SurfaceQuantity and VolumeQuantity - if isinstance(export, exports.SurfaceQuantity | exports.VolumeQuantity): + if isinstance( + export, + exports.SurfaceQuantity + | exports.VolumeQuantity + | exports.CustomQuantity, + ): export.t = [] export.data = [] + if isinstance(export, exports.CustomQuantity): + import ufl + import festim as F + + # <-- this won't work in the Discontinuous case we need to use species.subdomain_to_post_processing_solution instead + # we know the subdomain from the surface_to_subdomain attribute of the export, we can use that to get the post processing solution to use in the expression + kwargs = { + species.name: species.post_processing_solution + for species in self.species + } + kwargs["n"] = ufl.FacetNormal(self.mesh.mesh) + kwargs["T"] = self.temperature_fenics + + D_kwargs = { + f"D_{sp.name}": self.define_D_global(sp)[0] for sp in self.species + } + kwargs.update(D_kwargs) + kwargs["D"] = {sp.name: D_kwargs[f"D_{sp.name}"] for sp in self.species} + if len(self.species) == 1: + kwargs["D"] = kwargs[f"D_{self.species[0].name}"] + kwargs["x"] = ufl.SpatialCoordinate(self.mesh.mesh) + export.ufl_expr = export.expr(**kwargs) + def _get_temperature_field_as_function(self) -> dolfinx.fem.Function: """Based on the type of the temperature_fenics attribute, converts it as a Function to be used in VTX export. @@ -1031,6 +1065,27 @@ def post_processing(self): # if filename given write export data to file if export.filename is not None: export.write(t=float(self.t)) + elif isinstance(export, exports.CustomQuantity): + import festim as F + + if isinstance(export.subdomain, F.SurfaceSubdomain): + is_surface = True + elif type(export.subdomain) == int and any( + s.id == export.subdomain for s in self.surface_subdomains + ): + is_surface = True + else: + is_surface = False + measure = self.ds if is_surface else self.dx + export.compute(measure) + + # update export data + export.t.append(float(self.t)) + + # if filename given write export data to file + if export.filename is not None: + export.write(t=float(self.t)) + if isinstance(export, exports.XDMFExport): export.write(float(self.t)) @@ -1820,7 +1875,37 @@ def post_processing(self): else: export.compute() - if isinstance(export, exports.SurfaceQuantity | exports.VolumeQuantity): + elif isinstance(export, exports.CustomQuantity): + import festim as F + + if isinstance(export.subdomain, F.SurfaceSubdomain): + is_surface = True + elif type(export.subdomain) == int and any( + s.id == export.subdomain for s in self.surface_subdomains + ): + is_surface = True + else: + is_surface = False + measure = self.ds if is_surface else self.dx + + # getting entity_maps + try: + from dolfinx.mesh import EntityMap # noqa: F401 + + entity_maps = [sd.cell_map for sd in self.volume_subdomains] + except ImportError: + entity_maps = { + sd.submesh: sd.parent_to_submesh + for sd in self.volume_subdomains + } + export.compute(measure, entity_maps=entity_maps) + + if isinstance( + export, + exports.SurfaceQuantity + | exports.VolumeQuantity + | exports.CustomQuantity, + ): # update export data export.t.append(float(self.t)) From 75c94d77a178ad2a5ece6eb9506b7566f82255d7 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 14:21:19 -0400 Subject: [PATCH 03/15] integration with discontinuous --- src/festim/exports/custom_quantity.py | 2 +- src/festim/hydrogen_transport_problem.py | 63 +++++++++++++++++------- 2 files changed, 47 insertions(+), 18 deletions(-) diff --git a/src/festim/exports/custom_quantity.py b/src/festim/exports/custom_quantity.py index 8a337bbe6..94664ddc1 100644 --- a/src/festim/exports/custom_quantity.py +++ b/src/festim/exports/custom_quantity.py @@ -29,7 +29,7 @@ class CustomQuantity(DerivedQuantity): def __init__( self, expr: Callable, - subdomain: SurfaceSubdomain | VolumeSubdomain | int, + subdomain: SurfaceSubdomain | VolumeSubdomain, title: str = "Custom Quantity", filename: str | None = None, ) -> None: diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index d341ea0c5..18109438e 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1068,14 +1068,7 @@ def post_processing(self): elif isinstance(export, exports.CustomQuantity): import festim as F - if isinstance(export.subdomain, F.SurfaceSubdomain): - is_surface = True - elif type(export.subdomain) == int and any( - s.id == export.subdomain for s in self.surface_subdomains - ): - is_surface = True - else: - is_surface = False + is_surface = isinstance(export.subdomain, F.SurfaceSubdomain) measure = self.ds if is_surface else self.dx export.compute(measure) @@ -1770,10 +1763,53 @@ def initialise_exports(self): export.D = D # reset the data and time for SurfaceQuantity and VolumeQuantity - if isinstance(export, exports.SurfaceQuantity | exports.VolumeQuantity): + if isinstance( + export, + exports.SurfaceQuantity + | exports.VolumeQuantity + | exports.CustomQuantity, + ): export.t = [] export.data = [] + if isinstance(export, exports.CustomQuantity): + import ufl + import festim as F + + volume = ( + export.subdomain + if not isinstance(export.subdomain, F.SurfaceSubdomain) + else self.surface_to_volume[ + export.subdomain + if isinstance(export.subdomain, F.SurfaceSubdomain) + else next( + s + for s in self.surface_subdomains + if s.id == export.subdomain + ) + ] + ) + + kwargs = { + species.name: species.subdomain_to_post_processing_solution[volume] + for species in self.species + } + kwargs["n"] = ufl.FacetNormal(self.mesh.mesh) + kwargs["T"] = self.temperature_fenics + + D_kwargs = { + f"D_{sp.name}": volume.material.get_diffusion_coefficient( + self.mesh.mesh, self.temperature_fenics, sp + ) + for sp in self.species + } + kwargs.update(D_kwargs) + kwargs["D"] = {sp.name: D_kwargs[f"D_{sp.name}"] for sp in self.species} + if len(self.species) == 1: + kwargs["D"] = kwargs[f"D_{self.species[0].name}"] + kwargs["x"] = ufl.SpatialCoordinate(self.mesh.mesh) + export.ufl_expr = export.expr(**kwargs) + def post_processing(self): # update post-processing solutions (for each species in each subdomain) # with new solution @@ -1878,14 +1914,7 @@ def post_processing(self): elif isinstance(export, exports.CustomQuantity): import festim as F - if isinstance(export.subdomain, F.SurfaceSubdomain): - is_surface = True - elif type(export.subdomain) == int and any( - s.id == export.subdomain for s in self.surface_subdomains - ): - is_surface = True - else: - is_surface = False + is_surface = isinstance(export.subdomain, F.SurfaceSubdomain) measure = self.ds if is_surface else self.dx # getting entity_maps From 6fadd59d6870c60ac573547594095139d5e31f12 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 14:25:25 -0400 Subject: [PATCH 04/15] minor refactoring --- src/festim/exports/maximum_surface.py | 2 ++ src/festim/exports/maximum_volume.py | 2 ++ src/festim/exports/minimum_surface.py | 2 ++ src/festim/exports/minimum_volume.py | 2 ++ src/festim/hydrogen_transport_problem.py | 32 ++++++------------------ 5 files changed, 15 insertions(+), 25 deletions(-) diff --git a/src/festim/exports/maximum_surface.py b/src/festim/exports/maximum_surface.py index 1a29f5d7e..a65ca92be 100644 --- a/src/festim/exports/maximum_surface.py +++ b/src/festim/exports/maximum_surface.py @@ -19,6 +19,8 @@ class MaximumSurface(SurfaceQuantity): see `festim.SurfaceQuantity` """ + facet_meshtags: dolfinx.mesh.MeshTags + @property def title(self): return f"Maximum {self.field.name} surface {self.surface.id}" diff --git a/src/festim/exports/maximum_volume.py b/src/festim/exports/maximum_volume.py index d3bb257e3..2e3aa35a5 100644 --- a/src/festim/exports/maximum_volume.py +++ b/src/festim/exports/maximum_volume.py @@ -19,6 +19,8 @@ class MaximumVolume(VolumeQuantity): see `festim.VolumeQuantity` """ + volume_meshtags: dolfinx.mesh.MeshTags + @property def title(self): return f"Maximum {self.field.name} volume {self.volume.id}" diff --git a/src/festim/exports/minimum_surface.py b/src/festim/exports/minimum_surface.py index 1d2e934d5..468532c5f 100644 --- a/src/festim/exports/minimum_surface.py +++ b/src/festim/exports/minimum_surface.py @@ -19,6 +19,8 @@ class MinimumSurface(SurfaceQuantity): see `festim.SurfaceQuantity` """ + facet_meshtags: dolfinx.mesh.MeshTags + @property def title(self): return f"Minimum {self.field.name} surface {self.surface.id}" diff --git a/src/festim/exports/minimum_volume.py b/src/festim/exports/minimum_volume.py index e17ce5aff..a2cfc560c 100644 --- a/src/festim/exports/minimum_volume.py +++ b/src/festim/exports/minimum_volume.py @@ -19,6 +19,8 @@ class MinimumVolume(VolumeQuantity): see `festim.VolumeQuantity` """ + volume_meshtags: dolfinx.mesh.MeshTags + @property def title(self): return f"Minimum {self.field.name} volume {self.volume.id}" diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 18109438e..491f000d8 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -517,31 +517,17 @@ def initialise_exports(self): # add the global D to the export export.D = D export.D_expr = D_expr - if isinstance( - export, - exports.MaximumVolume - | exports.MaximumSurface - | exports.MinimumVolume - | exports.MinimumSurface, - ): + if isinstance(export, exports.MaximumVolume | exports.MinimumVolume): export.volume_meshtags = self.volume_meshtags + if isinstance(export, exports.MaximumSurface | exports.MinimumSurface): export.facet_meshtags = self.facet_meshtags + # reset the data and time for SurfaceQuantity and VolumeQuantity - if isinstance( - export, - exports.SurfaceQuantity - | exports.VolumeQuantity - | exports.CustomQuantity, - ): + if isinstance(export, exports.DerivedQuantity): export.t = [] export.data = [] if isinstance(export, exports.CustomQuantity): - import ufl - import festim as F - - # <-- this won't work in the Discontinuous case we need to use species.subdomain_to_post_processing_solution instead - # we know the subdomain from the surface_to_subdomain attribute of the export, we can use that to get the post processing solution to use in the expression kwargs = { species.name: species.post_processing_solution for species in self.species @@ -1763,17 +1749,13 @@ def initialise_exports(self): export.D = D # reset the data and time for SurfaceQuantity and VolumeQuantity - if isinstance( - export, - exports.SurfaceQuantity - | exports.VolumeQuantity - | exports.CustomQuantity, - ): + if isinstance(export, exports.DerivedQuantity): export.t = [] export.data = [] if isinstance(export, exports.CustomQuantity): import ufl + import festim as F volume = ( @@ -1892,7 +1874,7 @@ def post_processing(self): elif isinstance(export, exports.VolumeQuantity): if isinstance(export, exports.TotalVolume | exports.AverageVolume): try: - from dolfinx.mesh import EntityMap # noqa: F401 + from dolfinx.mesh import EntityMap entity_maps = [sd.cell_map for sd in self.volume_subdomains] except ImportError: From 6396a084928689c6d6802e673e89836d70eccc6c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 14:31:22 -0400 Subject: [PATCH 05/15] test --- test/test_derived_quantities.py | 362 ++++++++++++++++++++++++++++++++ 1 file changed, 362 insertions(+) diff --git a/test/test_derived_quantities.py b/test/test_derived_quantities.py index edfd6ae37..101115cbf 100644 --- a/test/test_derived_quantities.py +++ b/test/test_derived_quantities.py @@ -1,6 +1,10 @@ import os +from mpi4py import MPI + +import numpy as np import pytest +import ufl import festim as F @@ -46,3 +50,361 @@ def test_title(quantity, expected_title, tmp_path): quantity.value = 1 assert quantity.title == expected_title + + +class TestCustomQuantity: + """Test suite for CustomQuantity export""" + + def test_custom_quantity_title_volume(self): + """Test that CustomQuantity has correct title for volume subdomain""" + + def expr(**kwargs): + return kwargs.get("A", 0) + + quantity = F.CustomQuantity( + expr=expr, subdomain=vol_1, title="My Custom Quantity" + ) + assert quantity.title == "My Custom Quantity" + + def test_custom_quantity_title_surface(self): + """Test that CustomQuantity has correct title for surface subdomain""" + + def expr(**kwargs): + return kwargs.get("A", 0) + + quantity = F.CustomQuantity(expr=expr, subdomain=surf_1, title="Surface Custom") + assert quantity.title == "Surface Custom" + + def test_custom_quantity_subdomain_volume(self): + """Test that CustomQuantity stores volume subdomain correctly""" + + def expr(**kwargs): + return 1 + + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1) + assert quantity.subdomain == vol_1 + assert quantity.subdomain.id == 1 + + def test_custom_quantity_subdomain_surface(self): + """Test that CustomQuantity stores surface subdomain correctly""" + + def expr(**kwargs): + return 1 + + quantity = F.CustomQuantity(expr=expr, subdomain=surf_2) + assert quantity.subdomain == surf_2 + assert quantity.subdomain.id == 2 + + def test_custom_quantity_accepts_kwargs(self): + """Test that expression callable receives kwargs""" + received_kwargs = {} + + def expr(**kwargs): + received_kwargs.update(kwargs) + return 1 + + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1) + quantity.expr(**{"A": 1, "B": 2, "n": None}) + assert received_kwargs == {"A": 1, "B": 2, "n": None} + + def test_custom_quantity_data_append(self): + """Test that values are appended to data list""" + + def expr(**kwargs): + return 1 + + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1) + assert quantity.data == [] + quantity.value = 5 + quantity.data.append(quantity.value) + assert quantity.data == [5] + quantity.value = 10 + quantity.data.append(quantity.value) + assert quantity.data == [5, 10] + + def test_custom_quantity_t_list(self): + """Test that time values are stored in t list""" + + def expr(**kwargs): + return 1 + + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1) + assert quantity.t == [] + quantity.t.append(0.0) + quantity.t.append(0.1) + assert quantity.t == [0.0, 0.1] + + def test_custom_quantity_filename(self, tmp_path): + """Test that filename can be set and retrieved""" + + def expr(**kwargs): + return 1 + + filename = os.path.join(tmp_path, "custom.csv") + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1, filename=filename) + assert quantity.filename == filename + + def test_custom_quantity_expr_stored(self): + """Test that expression callable is stored correctly""" + + def my_expr(**kwargs): + return kwargs.get("value", 0) + + quantity = F.CustomQuantity(expr=my_expr, subdomain=vol_1) + assert quantity.expr == my_expr + assert quantity.expr(value=42) == 42 + + +class TestCustomQuantityWithHydrogenTransportProblem: + """Integration tests for CustomQuantity with HydrogenTransportProblem""" + + def test_custom_quantity_volume_integration(self): + """Test CustomQuantity with volume subdomain in HydrogenTransportProblem""" + # BUILD + material = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + vol = F.VolumeSubdomain(id=1, material=material) + 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)) + + from dolfinx.mesh import create_unit_square + + mesh = create_unit_square(MPI.COMM_WORLD, 5, 5) + + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh(mesh) + my_model.subdomains = [vol, top, bottom] + my_model.temperature = 300 + my_model.settings = F.Settings( + final_time=0.02, atol=1e-6, rtol=1e-6, stepsize=0.01 + ) + + species = F.Species("A") + my_model.species = [species] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC(species=species, subdomain=top, value=1), + F.FixedConcentrationBC(species=species, subdomain=bottom, value=0), + ] + + # Custom quantity that computes the total concentration + def total_conc(**kwargs): + return kwargs["A"] + + custom_qty = F.CustomQuantity( + expr=total_conc, subdomain=vol, title="Total concentration" + ) + + my_model.exports = [custom_qty] + + # RUN + my_model.initialise() + my_model.run() + + # TEST + assert len(custom_qty.data) > 0 + assert len(custom_qty.t) > 0 + assert len(custom_qty.t) == len(custom_qty.data) + # Check that values are non-negative (concentration) + assert all(v >= 0 for v in custom_qty.data) + + def test_custom_quantity_surface_integration(self): + """Test CustomQuantity with surface subdomain in HydrogenTransportProblem""" + # BUILD + material = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + vol = F.VolumeSubdomain(id=1, material=material) + 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)) + + from dolfinx.mesh import create_unit_square + + mesh = create_unit_square(MPI.COMM_WORLD, 5, 5) + + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh(mesh) + my_model.subdomains = [vol, top, bottom] + my_model.temperature = 300 + my_model.settings = F.Settings( + final_time=0.02, atol=1e-6, rtol=1e-6, stepsize=0.01 + ) + + species = F.Species("A") + my_model.species = [species] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC(species=species, subdomain=top, value=1), + F.FixedConcentrationBC(species=species, subdomain=bottom, value=0), + ] + + # Custom quantity on surface + def surface_flux(**kwargs): + D = kwargs["D_A"] + c = kwargs["A"] + n = kwargs["n"] + return -D * ufl.dot(ufl.grad(c), n) + + custom_qty = F.CustomQuantity( + expr=surface_flux, subdomain=top, title="Surface flux" + ) + + my_model.exports = [custom_qty] + + # RUN + my_model.initialise() + my_model.run() + + # TEST + assert len(custom_qty.data) > 0 + assert len(custom_qty.t) > 0 + assert len(custom_qty.t) == len(custom_qty.data) + + def test_custom_quantity_with_multiple_species(self): + """Test CustomQuantity with multiple species""" + # BUILD + material = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + vol = F.VolumeSubdomain(id=1, material=material) + 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)) + + from dolfinx.mesh import create_unit_square + + mesh = create_unit_square(MPI.COMM_WORLD, 5, 5) + + my_model = F.HydrogenTransportProblem() + my_model.mesh = F.Mesh(mesh) + my_model.subdomains = [vol, top, bottom] + my_model.temperature = 300 + my_model.settings = F.Settings( + final_time=0.02, atol=1e-6, rtol=1e-6, stepsize=0.01 + ) + + speciesA = F.Species("A") + speciesB = F.Species("B") + my_model.species = [speciesA, speciesB] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC(species=speciesA, subdomain=top, value=1), + F.FixedConcentrationBC(species=speciesA, subdomain=bottom, value=0), + F.FixedConcentrationBC(species=speciesB, subdomain=top, value=0), + F.FixedConcentrationBC(species=speciesB, subdomain=bottom, value=1), + ] + + # Custom quantity combining both species + def combined_conc(**kwargs): + return kwargs["A"] + kwargs["B"] + + custom_qty = F.CustomQuantity( + expr=combined_conc, subdomain=vol, title="Combined concentration" + ) + + my_model.exports = [custom_qty] + + # RUN + my_model.initialise() + my_model.run() + + # TEST + assert len(custom_qty.data) > 0 + assert len(custom_qty.t) > 0 + # Combined concentration should be <= 2 (max of each species) + assert all(0 <= v <= 2 for v in custom_qty.data) + + +class TestCustomQuantityWithDiscontinuousProblem: + """Integration tests for CustomQuantity with HydrogenTransportProblemDiscontinuous""" + + def test_custom_quantity_discontinuous_volume(self): + """Test CustomQuantity with volume subdomain in discontinuous problem""" + # BUILD + material = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + vol = F.VolumeSubdomain(id=1, material=material) + 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)) + + from dolfinx.mesh import create_unit_square + + mesh = create_unit_square(MPI.COMM_WORLD, 5, 5) + + my_model = F.HydrogenTransportProblemDiscontinuous() + my_model.mesh = F.Mesh(mesh) + my_model.subdomains = [vol, top, bottom] + my_model.temperature = 300 + my_model.settings = F.Settings( + final_time=0.02, atol=1e-6, rtol=1e-6, stepsize=0.01 + ) + + species = F.Species("A", subdomains=[vol]) + my_model.species = [species] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC(species=species, subdomain=top, value=1), + F.FixedConcentrationBC(species=species, subdomain=bottom, value=0), + ] + + # Custom quantity + def total_conc(**kwargs): + return kwargs["A"] + + custom_qty = F.CustomQuantity( + expr=total_conc, subdomain=vol, title="Total concentration" + ) + + my_model.exports = [custom_qty] + + # RUN + my_model.initialise() + my_model.run() + + # TEST + assert len(custom_qty.data) > 0 + assert len(custom_qty.t) > 0 + assert len(custom_qty.t) == len(custom_qty.data) + + def test_custom_quantity_discontinuous_surface(self): + """Test CustomQuantity with surface subdomain in discontinuous problem""" + # BUILD + material = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + vol = F.VolumeSubdomain(id=1, material=material) + 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)) + + from dolfinx.mesh import create_unit_square + + mesh = create_unit_square(MPI.COMM_WORLD, 5, 5) + + my_model = F.HydrogenTransportProblemDiscontinuous() + my_model.mesh = F.Mesh(mesh) + my_model.subdomains = [vol, top, bottom] + my_model.temperature = 300 + my_model.settings = F.Settings( + final_time=0.02, atol=1e-6, rtol=1e-6, stepsize=0.01 + ) + + species = F.Species("A", subdomains=[vol]) + my_model.species = [species] + + my_model.boundary_conditions = [ + F.FixedConcentrationBC(species=species, subdomain=top, value=1), + F.FixedConcentrationBC(species=species, subdomain=bottom, value=0), + ] + + # Custom quantity on surface + def surface_flux(**kwargs): + D = kwargs["D_A"] + c = kwargs["A"] + n = kwargs["n"] + return -D * ufl.dot(ufl.grad(c), n) + + custom_qty = F.CustomQuantity( + expr=surface_flux, subdomain=top, title="Surface flux" + ) + + my_model.exports = [custom_qty] + + # RUN + my_model.initialise() + my_model.run() + + # TEST + assert len(custom_qty.data) > 0 + assert len(custom_qty.t) > 0 + assert len(custom_qty.t) == len(custom_qty.data) From 7fb33cd7ee95f128c05031ec823f453a52bc0952 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 14:35:57 -0400 Subject: [PATCH 06/15] more refactoring --- src/festim/hydrogen_transport_problem.py | 14 ++------------ 1 file changed, 2 insertions(+), 12 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 491f000d8..5bee150a9 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -462,12 +462,7 @@ def initialise_exports(self): ) continue - elif isinstance( - export, - exports.SurfaceQuantity - | exports.VolumeQuantity - | exports.CustomQuantity, - ): + elif isinstance(export, exports.DerivedQuantity): # raise not implemented error if the derived quantity don't match the # type of mesh eg. SurfaceFlux is used with cylindrical mesh if self.mesh.coordinate_system != CoordinateSystem.CARTESIAN: @@ -1911,12 +1906,7 @@ def post_processing(self): } export.compute(measure, entity_maps=entity_maps) - if isinstance( - export, - exports.SurfaceQuantity - | exports.VolumeQuantity - | exports.CustomQuantity, - ): + if isinstance(export, exports.DerivedQuantity): # update export data export.t.append(float(self.t)) From 793655a7e45cd834d7f4a5eb459ac3997f98963a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 14:36:54 -0400 Subject: [PATCH 07/15] ruff --- src/festim/exports/custom_quantity.py | 5 +++-- src/festim/exports/volume_quantity.py | 2 +- test/test_derived_quantities.py | 4 +++- test/test_vtx.py | 4 ++-- 4 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/festim/exports/custom_quantity.py b/src/festim/exports/custom_quantity.py index 94664ddc1..c49395fc1 100644 --- a/src/festim/exports/custom_quantity.py +++ b/src/festim/exports/custom_quantity.py @@ -1,7 +1,8 @@ -from typing import Callable +from collections.abc import Callable + import ufl -from scifem import assemble_scalar from dolfinx import fem +from scifem import assemble_scalar from festim.exports.derived_quantity import DerivedQuantity from festim.subdomain.surface_subdomain import SurfaceSubdomain diff --git a/src/festim/exports/volume_quantity.py b/src/festim/exports/volume_quantity.py index a44f64606..153a6c3f8 100644 --- a/src/festim/exports/volume_quantity.py +++ b/src/festim/exports/volume_quantity.py @@ -1,6 +1,6 @@ +from festim.exports.derived_quantity import DerivedQuantity from festim.species import Species from festim.subdomain.volume_subdomain import VolumeSubdomain -from festim.exports.derived_quantity import DerivedQuantity class VolumeQuantity(DerivedQuantity): diff --git a/test/test_derived_quantities.py b/test/test_derived_quantities.py index 101115cbf..1c6505bda 100644 --- a/test/test_derived_quantities.py +++ b/test/test_derived_quantities.py @@ -310,7 +310,9 @@ def combined_conc(**kwargs): class TestCustomQuantityWithDiscontinuousProblem: - """Integration tests for CustomQuantity with HydrogenTransportProblemDiscontinuous""" + """ + Integration tests for CustomQuantity with HydrogenTransportProblemDiscontinuous + """ def test_custom_quantity_discontinuous_volume(self): """Test CustomQuantity with volume subdomain in discontinuous problem""" diff --git a/test/test_vtx.py b/test/test_vtx.py index 9c8bf27c5..b372d96e6 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -411,8 +411,8 @@ def test_custom_field_not_implemented_error(expression): @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 ReactionRateExport export functionality for different directions, product formats, - and reaction configurations. + Test ReactionRateExport export functionality for different directions, product + formats, and reaction configurations. """ if p_0 == 0.0 and direction == "backward": pytest.skip( From a007328fc77a9cd1c58f7a234270eb693e6c786f Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 14:45:51 -0400 Subject: [PATCH 08/15] documentation and usage --- src/festim/exports/custom_quantity.py | 55 ++++++++++++++++++++++++++- 1 file changed, 54 insertions(+), 1 deletion(-) diff --git a/src/festim/exports/custom_quantity.py b/src/festim/exports/custom_quantity.py index c49395fc1..6e0b41d48 100644 --- a/src/festim/exports/custom_quantity.py +++ b/src/festim/exports/custom_quantity.py @@ -10,7 +10,8 @@ class CustomQuantity(DerivedQuantity): - """Export CustomQuantity. + """ + Export CustomQuantity. Args: expr: function that returns a UFL expression @@ -25,6 +26,58 @@ class CustomQuantity(DerivedQuantity): filename: name of the file to which the quantity is exported t: list of time values data: list of values of the quantity + + Usage: + + .. testcode:: + + import numpy as np + import festim as F + + material = F.Material(D_0=1, E_D=0) + volume = F.VolumeSubdomain(id=1, material=material) + surface = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1)) + + def total_concentration(**kwargs): + return kwargs["A"] + kwargs["B"] + + quantity = F.CustomQuantity( + expr=total_concentration, + subdomain=volume, + title="Total quantity", + ) + + surface_quantity = F.CustomQuantity( + expr=lambda **kwargs: -kwargs["D_A"] * ufl.dot( + ufl.grad(kwargs["A"]), kwargs["n"] + ), + subdomain=surface, + title="Surface flux", + ) + + For a surface quantity, the returned UFL expression can represent a flux such as + + .. math:: + + q = -D\,\nabla c \cdot n + + and FESTIM will assemble + + .. math:: + + Q = \int_{\Gamma} q\,\mathrm{d}\Gamma + + over the selected surface subdomain. + + The expression returned by ``expr`` is treated as an integrand and assembled over + the chosen subdomain. + + .. math:: + + Q = \int_{\Omega} q\,\mathrm{d}\Omega + + where ``q`` is the UFL expression returned by ``expr`` and ``\Omega`` is either a + surface or a volume subdomain. """ def __init__( From bd2cd3823abc1091206447bd0c9927902fabcc7f Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 14:51:04 -0400 Subject: [PATCH 09/15] docs for the expr --- src/festim/exports/custom_quantity.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/src/festim/exports/custom_quantity.py b/src/festim/exports/custom_quantity.py index 6e0b41d48..4d1df7679 100644 --- a/src/festim/exports/custom_quantity.py +++ b/src/festim/exports/custom_quantity.py @@ -10,7 +10,7 @@ class CustomQuantity(DerivedQuantity): - """ + r""" Export CustomQuantity. Args: @@ -55,6 +55,23 @@ def total_concentration(**kwargs): title="Surface flux", ) + The callable passed to ``expr`` receives keyword arguments assembled by the + problem class. Common entries are: + + ``A``, ``B``, ... + Concentrations of the species present in the problem (here A and B). + ``n`` + The facet normal on the selected surface subdomain. + ``T`` + The temperature field. + ``D_A``, ``D_B``, ... + Species-specific diffusion coefficients. + ``D`` + The diffusion coefficient data, either a single field for one species or + a dictionary keyed by species name when several species are present. + ``x`` + The spatial coordinate (x[0], x[1], x[2]). + For a surface quantity, the returned UFL expression can represent a flux such as .. math:: From ec16a28d63a7b12285a58b5c0ddd2f86f010f137 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 14:55:15 -0400 Subject: [PATCH 10/15] abstract method --- src/festim/exports/derived_quantity.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/festim/exports/derived_quantity.py b/src/festim/exports/derived_quantity.py index 475fe51e7..8b2625209 100644 --- a/src/festim/exports/derived_quantity.py +++ b/src/festim/exports/derived_quantity.py @@ -1,7 +1,8 @@ import csv +from abc import ABC, abstractmethod -class DerivedQuantity: +class DerivedQuantity(ABC): """Base class for all derived quantities. Attributes: @@ -21,8 +22,13 @@ def __init__(self, filename: str | None = None) -> None: self._first_time_export = True @property + @abstractmethod def title(self): - raise NotImplementedError("title must be implemented by subclasses") + pass + + @abstractmethod + def compute(self, *args, **kwargs): + pass @property def filename(self): From 4cb577d4cbfa7edc2fe27933ed0f8cd90301e12d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 15:07:17 -0400 Subject: [PATCH 11/15] Volume and Surface Quantity are also base classes + adapted tests --- src/festim/exports/surface_quantity.py | 11 +++ src/festim/exports/volume_quantity.py | 11 +++ test/test_derived_quantities.py | 105 ++++++++++++++++++++++++- test/test_surface_quantity.py | 35 --------- test/test_volume_quantity.py | 23 ------ 5 files changed, 125 insertions(+), 60 deletions(-) diff --git a/src/festim/exports/surface_quantity.py b/src/festim/exports/surface_quantity.py index e915aebd6..1c2d0b06e 100644 --- a/src/festim/exports/surface_quantity.py +++ b/src/festim/exports/surface_quantity.py @@ -1,3 +1,5 @@ +from abc import abstractmethod + from festim.exports.derived_quantity import DerivedQuantity from festim.species import Species from festim.subdomain.surface_subdomain import SurfaceSubdomain @@ -36,6 +38,15 @@ def __init__( self.field = field self.surface = surface + @property + @abstractmethod + def title(self): + pass + + @abstractmethod + def compute(self, *args, **kwargs): + pass + @property def surface(self): return self._surface diff --git a/src/festim/exports/volume_quantity.py b/src/festim/exports/volume_quantity.py index 153a6c3f8..7ebcace4c 100644 --- a/src/festim/exports/volume_quantity.py +++ b/src/festim/exports/volume_quantity.py @@ -1,3 +1,5 @@ +from abc import abstractmethod + from festim.exports.derived_quantity import DerivedQuantity from festim.species import Species from festim.subdomain.volume_subdomain import VolumeSubdomain @@ -36,6 +38,15 @@ def __init__( self.field = field self.volume = volume + @property + @abstractmethod + def title(self): + pass + + @abstractmethod + def compute(self, *args, **kwargs): + pass + @property def field(self): return self._field diff --git a/test/test_derived_quantities.py b/test/test_derived_quantities.py index 1c6505bda..fd6945b29 100644 --- a/test/test_derived_quantities.py +++ b/test/test_derived_quantities.py @@ -26,8 +26,6 @@ max_surface = F.MaximumSurface(mobile_H, surface=surf_2, filename=results) avg_surface = F.AverageSurface(mobile_D, surface=surf_1, filename=results) avg_vol = F.AverageVolume(mobile_H, volume=vol_2, filename=results) -surf_quant = F.SurfaceQuantity(mobile_H, surface=surf_1, filename=results) -vol_quant = F.VolumeQuantity(mobile_H, volume=vol_1, filename=results) @pytest.mark.parametrize( @@ -309,6 +307,109 @@ def combined_conc(**kwargs): assert all(0 <= v <= 2 for v in custom_qty.data) +class TestDerivedQuantityFilename: + """Test suite for DerivedQuantity filename property validation""" + + def test_filename_none(self): + """Test that filename can be set to None""" + + def expr(**kwargs): + return 1 + + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1, filename=None) + assert quantity.filename is None + + def test_filename_csv(self, tmp_path): + """Test that filename accepts .csv extension""" + + def expr(**kwargs): + return 1 + + filename = os.path.join(tmp_path, "output.csv") + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1, filename=filename) + assert quantity.filename == filename + + def test_filename_txt(self, tmp_path): + """Test that filename accepts .txt extension""" + + def expr(**kwargs): + return 1 + + filename = os.path.join(tmp_path, "output.txt") + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1, filename=filename) + assert quantity.filename == filename + + def test_filename_invalid_extension(self): + """Test that filename rejects invalid extensions""" + + def expr(**kwargs): + return 1 + + with pytest.raises(ValueError, match="filename must end with .csv or .txt"): # noqa: RUF043 + F.CustomQuantity(expr=expr, subdomain=vol_1, filename="output.dat") + + def test_filename_invalid_type(self): + """Test that filename rejects non-string types""" + + def expr(**kwargs): + return 1 + + with pytest.raises(TypeError, match="filename must be of type str"): + F.CustomQuantity(expr=expr, subdomain=vol_1, filename=123) + + def test_filename_property_setter_none(self): + """Test setting filename property to None after initialization""" + + def expr(**kwargs): + return 1 + + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1) + quantity.filename = None + assert quantity.filename is None + + def test_filename_property_setter_csv(self, tmp_path): + """Test setting filename property to a .csv file""" + + def expr(**kwargs): + return 1 + + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1) + filename = os.path.join(tmp_path, "new_output.csv") + quantity.filename = filename + assert quantity.filename == filename + + def test_filename_property_setter_txt(self, tmp_path): + """Test setting filename property to a .txt file""" + + def expr(**kwargs): + return 1 + + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1) + filename = os.path.join(tmp_path, "new_output.txt") + quantity.filename = filename + assert quantity.filename == filename + + def test_filename_property_setter_invalid_extension(self): + """Test that setting invalid extension raises ValueError""" + + def expr(**kwargs): + return 1 + + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1) + with pytest.raises(ValueError, match="filename must end with .csv or .txt"): # noqa: RUF043 + quantity.filename = "output.json" + + def test_filename_property_setter_invalid_type(self): + """Test that setting non-string filename raises TypeError""" + + def expr(**kwargs): + return 1 + + quantity = F.CustomQuantity(expr=expr, subdomain=vol_1) + with pytest.raises(TypeError, match="filename must be of type str"): + quantity.filename = 42 + + class TestCustomQuantityWithDiscontinuousProblem: """ Integration tests for CustomQuantity with HydrogenTransportProblemDiscontinuous diff --git a/test/test_surface_quantity.py b/test/test_surface_quantity.py index c2943bcaf..a94378d44 100644 --- a/test/test_surface_quantity.py +++ b/test/test_surface_quantity.py @@ -101,29 +101,6 @@ def test_write_overwrite(tmp_path): assert file_length == expected_length -def test_filename_setter_raises_TypeError(): - """Test that a TypeError is raised when the filename is not a string.""" - - with pytest.raises(TypeError, match="filename must be of type str"): - F.SurfaceQuantity( - filename=1, - field=F.Species("test"), - surface=F.SurfaceSubdomain1D(id=1, x=0), - ) - - -def test_filename_setter_raises_ValueError(tmp_path): - """Test that a ValueError is raised when the filename does not end with .csv or - .txt.""" - - with pytest.raises(ValueError): - F.SurfaceQuantity( - filename=os.path.join(tmp_path, "my_export.xdmf"), - field=F.Species("test"), - surface=F.SurfaceSubdomain1D(id=1, x=0), - ) - - def test_field_setter_raises_TypeError(): """Test that a TypeError is raised when the field is not a F.Species.""" @@ -151,15 +128,3 @@ def test_writer(tmp_path, value): expected_length = i + 2 assert file_length == expected_length - - -def test_surface_setter_raises_TypeError(): - """Test that a TypeError is raised when the surface is not a F.SurfaceSubdomain.""" - - with pytest.raises( - TypeError, match=r"surface should be an int or F.SurfaceSubdomain" - ): - F.SurfaceQuantity( - field=F.Species("H"), - surface="1", - ) diff --git a/test/test_volume_quantity.py b/test/test_volume_quantity.py index 7f530ef58..eb3a44478 100644 --- a/test/test_volume_quantity.py +++ b/test/test_volume_quantity.py @@ -6,29 +6,6 @@ dummy_volume = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=dummy_mat) -def test_filename_setter_raises_TypeError(): - """Test that a TypeError is raised when the filename is not a string.""" - - with pytest.raises(TypeError, match="filename must be of type str"): - F.VolumeQuantity( - filename=1, - field=F.Species("test"), - volume=dummy_volume, - ) - - -def test_filename_setter_raises_ValueError(): - """Test that a ValueError is raised when the filename does not end with .csv or - .txt.""" - - with pytest.raises(ValueError): - F.VolumeQuantity( - filename="my_export.xdmf", - field=F.Species("test"), - volume=dummy_volume, - ) - - def test_field_setter_raises_TypeError(): """Test that a TypeError is raised when the field is not a F.Species.""" From 7f1ca11771fe6acb23ec197769f0d35db90b3884 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 15:47:22 -0400 Subject: [PATCH 12/15] fix for time-dependent temperature --- src/festim/hydrogen_transport_problem.py | 57 +++++++++------------ test/system_tests/test_surface_reactions.py | 2 +- 2 files changed, 24 insertions(+), 35 deletions(-) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 5bee150a9..993656e15 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -135,6 +135,8 @@ class HydrogenTransportProblem(problem.ProblemBase): """ _temperature_as_function: fem.Function + _species_to_D_global: dict[_species.Species, fem.Function] + _species_to_D_global_expr: dict[_species.Species, fem.Expression] def __init__( self, @@ -190,6 +192,8 @@ def __init__( self._element_immobile = element_immobile self._temperature_as_function = None + self._species_to_D_global = None + self._species_to_D_global_expr = None @property def temperature(self): @@ -494,24 +498,21 @@ def initialise_exports(self): export.t = [] # compute diffusivity function for surface fluxes - spe_to_D_global = {} # links species to global D function - spe_to_D_global_expr = {} # links species to D expression + spe_to_D_global_func_expr = { + spe: self.define_D_global(spe) for spe in self.species if spe.mobile + } + self._species_to_D_global_expr = { + k: v[1] for k, v in spe_to_D_global_func_expr.items() + } # links species to D expression + self._species_to_D_global = { + k: v[0] for k, v in spe_to_D_global_func_expr.items() + } # links species to global D function for export in self.exports: if isinstance(export, exports.SurfaceQuantity): - if export.field in spe_to_D_global: - # if already computed then use the same D - D = spe_to_D_global[export.field] - D_expr = spe_to_D_global_expr[export.field] - else: - # compute D and add it to the dict - D, D_expr = self.define_D_global(export.field) - spe_to_D_global[export.field] = D - spe_to_D_global_expr[export.field] = D_expr - # add the global D to the export - export.D = D - export.D_expr = D_expr + export.D = self._species_to_D_global.get(export.field) + export.D_expr = self._species_to_D_global_expr.get(export.field) if isinstance(export, exports.MaximumVolume | exports.MinimumVolume): export.volume_meshtags = self.volume_meshtags if isinstance(export, exports.MaximumSurface | exports.MinimumSurface): @@ -530,8 +531,9 @@ def initialise_exports(self): kwargs["n"] = ufl.FacetNormal(self.mesh.mesh) kwargs["T"] = self.temperature_fenics + # NOTE we need to change our D_global approach D_kwargs = { - f"D_{sp.name}": self.define_D_global(sp)[0] for sp in self.species + f"D_{sp.name}": self._species_to_D_global[sp] for sp in self.species } kwargs.update(D_kwargs) kwargs["D"] = {sp.name: D_kwargs[f"D_{sp.name}"] for sp in self.species} @@ -973,13 +975,8 @@ def post_processing(self): if self.temperature_time_dependent: # update global D if temperature time dependent or internal # variables time dependent - species_not_updated = self.species.copy() # make a copy of the species - for export in self.exports: - if isinstance(export, exports.SurfaceFlux): - # if the D of the species has not been updated yet - if export.field in species_not_updated: - export.D.interpolate(export.D_expr) - species_not_updated.remove(export.field) + for spe, D_global in self._species_to_D_global.items(): + D_global.interpolate(self._species_to_D_global_expr[spe]) for export in self.exports: # skip if it isn't time to export @@ -1047,9 +1044,7 @@ def post_processing(self): if export.filename is not None: export.write(t=float(self.t)) elif isinstance(export, exports.CustomQuantity): - import festim as F - - is_surface = isinstance(export.subdomain, F.SurfaceSubdomain) + is_surface = isinstance(export.subdomain, _subdomain.SurfaceSubdomain) measure = self.ds if is_surface else self.dx export.compute(measure) @@ -1749,16 +1744,12 @@ def initialise_exports(self): export.data = [] if isinstance(export, exports.CustomQuantity): - import ufl - - import festim as F - volume = ( export.subdomain - if not isinstance(export.subdomain, F.SurfaceSubdomain) + if not isinstance(export.subdomain, _subdomain.SurfaceSubdomain) else self.surface_to_volume[ export.subdomain - if isinstance(export.subdomain, F.SurfaceSubdomain) + if isinstance(export.subdomain, _subdomain.SurfaceSubdomain) else next( s for s in self.surface_subdomains @@ -1889,9 +1880,7 @@ def post_processing(self): export.compute() elif isinstance(export, exports.CustomQuantity): - import festim as F - - is_surface = isinstance(export.subdomain, F.SurfaceSubdomain) + is_surface = isinstance(export.subdomain, _subdomain.SurfaceSubdomain) measure = self.ds if is_surface else self.dx # getting entity_maps diff --git a/test/system_tests/test_surface_reactions.py b/test/system_tests/test_surface_reactions.py index 7efaa2bad..dedc4013c 100644 --- a/test/system_tests/test_surface_reactions.py +++ b/test/system_tests/test_surface_reactions.py @@ -9,7 +9,7 @@ class FluxFromSurfaceReaction(F.SurfaceFlux): def __init__(self, reaction: F.SurfaceReactionBC): super().__init__( - F.Species(), # just a dummy species here + F.Species(mobile=False), # just a dummy species here reaction.subdomain, ) self.reaction = reaction.flux_bcs[0] From 15baccff5d6a8b68dbd9a675543fea094f284781 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 15:50:15 -0400 Subject: [PATCH 13/15] comment --- src/festim/hydrogen_transport_problem.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 993656e15..bca04f763 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -498,6 +498,8 @@ def initialise_exports(self): export.t = [] # compute diffusivity function for surface fluxes + # TODO: probably a better way to handle things would be to follow what's done in + # https://jsdokken.com/dolfinx-tutorial/chapter3/subdomains.html spe_to_D_global_func_expr = { spe: self.define_D_global(spe) for spe in self.species if spe.mobile } From 78526817f188fcf3761aca926d9653d194919be7 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Apr 2026 15:52:29 -0400 Subject: [PATCH 14/15] added new TODO --- src/festim/hydrogen_transport_problem.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index bca04f763..4d1fc829a 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -977,6 +977,9 @@ def post_processing(self): if self.temperature_time_dependent: # update global D if temperature time dependent or internal # variables time dependent + # TODO: honestly, we probably don't need to do this at all + # SurfaceFlux quantities should use ufl.Expr for D instead of a fem.Function + for spe, D_global in self._species_to_D_global.items(): D_global.interpolate(self._species_to_D_global_expr[spe]) From 67828016a801f06e5efaf14368263d2d9d76126e Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 8 May 2026 13:08:00 -0400 Subject: [PATCH 15/15] added time in kwargs --- src/festim/hydrogen_transport_problem.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 4d1fc829a..2a0a54162 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -531,6 +531,7 @@ def initialise_exports(self): for species in self.species } kwargs["n"] = ufl.FacetNormal(self.mesh.mesh) + kwargs["t"] = self.t kwargs["T"] = self.temperature_fenics # NOTE we need to change our D_global approach @@ -1768,6 +1769,7 @@ def initialise_exports(self): for species in self.species } kwargs["n"] = ufl.FacetNormal(self.mesh.mesh) + kwargs["t"] = self.t kwargs["T"] = self.temperature_fenics D_kwargs = {