Coupled transient heat transfer and hydrogen transport models

src/festim/__init__.py

@@ -18,10 +18,11 @@
 )
 from .boundary_conditions.flux_bc import FluxBCBase, HeatFluxBC, ParticleFluxBC
 from .boundary_conditions.henrys_bc import HenrysBC
-from .boundary_conditions.flux_bc import FluxBCBase, ParticleFluxBC, HeatFluxBC
-from .boundary_conditions.surface_reaction import SurfaceReactionBC
-
 from .boundary_conditions.sieverts_bc import SievertsBC
+from .boundary_conditions.surface_reaction import SurfaceReactionBC
+from .coupled_heat_hydrogen_problem import (
+    CoupledTransientHeatTransferHydrogenTransport,
+)
 from .exports.average_surface import AverageSurface
 from .exports.average_volume import AverageVolume
 from .exports.maximum_surface import MaximumSurface
@@ -45,20 +46,20 @@
 )
 from .hydrogen_transport_problem import (
     HTransportProblemDiscontinuous,
-    HydrogenTransportProblem,
     HTransportProblemPenalty,
+    HydrogenTransportProblem,
 )
-from .problem_change_of_var import HydrogenTransportProblemDiscontinuousChangeVar
 from .initial_condition import InitialCondition, InitialTemperature
 from .material import Material
 from .mesh.mesh import Mesh
 from .mesh.mesh_1d import Mesh1D
 from .mesh.mesh_from_xdmf import MeshFromXDMF
 from .problem import ProblemBase
+from .problem_change_of_var import HydrogenTransportProblemDiscontinuousChangeVar
 from .reaction import Reaction
 from .settings import Settings
 from .source import HeatSource, ParticleSource, SourceBase
-from .species import ImplicitSpecies, Species, find_species_from_name, SpeciesChangeVar
+from .species import ImplicitSpecies, Species, SpeciesChangeVar, find_species_from_name
 from .stepsize import Stepsize
 from .subdomain.interface import Interface
 from .subdomain.surface_subdomain import SurfaceSubdomain, find_surface_from_id

src/festim/coupled_heat_hydrogen_problem.py

@@ -0,0 +1,178 @@
+import tqdm.autonotebook
+from dolfinx import fem
+
+from festim.heat_transfer_problem import HeatTransferProblem
+from festim.helpers import as_fenics_constant, nmm_interpolate
+from festim.hydrogen_transport_problem import (
+    HTransportProblemDiscontinuous,
+    HTransportProblemPenalty,
+    HydrogenTransportProblem,
+)
+from festim.problem_change_of_var import HydrogenTransportProblemDiscontinuousChangeVar
+
+
+class CoupledTransientHeatTransferHydrogenTransport:
+    """
+    Coupled heat transfer and hydrogen transport transient problem
+
+    Args:
+        heat_problem: the heat transfer problem
+        hydrogen_problem: the hydrogen transport problem
+
+    Attributes:
+        heat_problem: the heat transfer problem
+        hydrogen_problem: the hydrogen transport problem
+        non_matching_meshes: True if the meshes in the heat_problem and hydrogen_problem
+            are not matching
+
+    Examples:
+        .. highlight:: python
+        .. code-block:: python
+
+            import festim as F
+
+            my_heat_transfer_model = F.HeatTransferProblem(
+                mesh=F.Mesh(...),
+                subdomains=[F.Subdomain(...)],
+                ...
+            )
+
+            my_h_transport_model = F.HydrogenTransportProblem(
+                mesh=F.Mesh(...),
+                subdomains=[F.Subdomain(...)],
+                species=[F.Species(name="H"), F.Species(name="Trap")],
+                ...
+            )
+
+            coupled_problem = F.CoupledHeatTransferHydrogenTransport(
+                heat_problem=my_heat_transfer_model,
+                hydrogen_problem=my_h_transport_model,
+            )
+
+    """
+
+    heat_problem: HeatTransferProblem
+    hydrogen_problem: HydrogenTransportProblem
+
+    non_matching_meshes: bool
+
+    def __init__(
+        self,
+        heat_problem: HeatTransferProblem,
+        hydrogen_problem: HydrogenTransportProblem,
+    ) -> None:
+        self.heat_problem = heat_problem
+        self.hydrogen_problem = hydrogen_problem
+
+        if (
+            not self.heat_problem.settings.transient
+            or not self.hydrogen_problem.settings.transient
+        ):
+            raise TypeError(
+                "Both the heat and hydrogen problems must be set to transient"
+            )
+
+    @property
+    def heat_problem(self):
+        return self._heat_problem
+
+    @heat_problem.setter
+    def heat_problem(self, value):
+        if not isinstance(value, HeatTransferProblem):
+            raise TypeError("heat_problem must be a festim.HeatTransferProblem object")
+        value.show_progress_bar = False
+        self._heat_problem = value
+
+    @property
+    def hydrogen_problem(self):
+        return self._hydrogen_problem
+
+    @hydrogen_problem.setter
+    def hydrogen_problem(self, value):
+        if isinstance(
+            value,
+            HTransportProblemDiscontinuous
+            | HTransportProblemPenalty
+            | HydrogenTransportProblemDiscontinuousChangeVar,
+        ):
+            raise NotImplementedError(
+                "Coupled heat transfer - hydrogen transport simulations with "
+                "HydrogenTransportProblemDiscontinuousChangeVar, "
+                "HTransportProblemPenalty or"
+                "HydrogenTransportProblemDiscontinuousChangeVar, "
+                "not currently supported"
+            )
+        elif not isinstance(value, HydrogenTransportProblem):
+            raise TypeError(
+                "hydrogen_problem must be a festim.HydrogenTransportProblem object"
+            )
+        self._hydrogen_problem = value
+
+    @property
+    def non_matching_meshes(self):
+        return self.heat_problem.mesh.mesh != self.hydrogen_problem.mesh.mesh
+
+    def initialise(self):
+        # make sure both problems have the same initial time step and final time,
+        # use minimal initial value of the two and maximal final time of the two
+        min_initial_dt = min(
+            self.heat_problem.settings.stepsize.initial_value,
+            self.hydrogen_problem.settings.stepsize.initial_value,
+        )
+        self.heat_problem.settings.stepsize.initial_value = min_initial_dt
+        self.hydrogen_problem.settings.stepsize.initial_value = min_initial_dt
+
+        if (
+            self.heat_problem.settings.final_time
+            != self.hydrogen_problem.settings.final_time
+        ):
+            raise ValueError(
+                "Final time values in the heat transfer and hydrogen transport "
+                "model must be the same"
+            )
+
+        self.heat_problem.initialise()
+
+        if self.non_matching_meshes:
+            V = fem.functionspace(self.hydrogen_problem.mesh.mesh, ("P", 1))
+            T_func = fem.Function(V)
+
+            nmm_interpolate(T_func, self.heat_problem.u)
+
+            self.hydrogen_problem.temperature = T_func
+        else:
+            self.hydrogen_problem.temperature = self.heat_problem.u
+
+        self.hydrogen_problem.initialise()
+
+    def iterate(self):
+        self.heat_problem.iterate()
+
+        if self.non_matching_meshes:
+            nmm_interpolate(
+                self.hydrogen_problem.temperature_fenics, self.heat_problem.u
+            )
+
+        self.hydrogen_problem.iterate()
+
+        # use the same time step for both problems, use minimum of the two
+        next_dt_value = min(
+            float(self.hydrogen_problem.dt), float(self.heat_problem.dt)
+        )
+        self.heat_problem.dt.value = next_dt_value
+        self.hydrogen_problem.dt.value = next_dt_value
+
+    def run(self):
+        if self.hydrogen_problem.show_progress_bar:
+            self.hydrogen_problem.progress_bar = tqdm.autonotebook.tqdm(
+                desc=f"Solving {self.__class__.__name__}",
+                total=self.hydrogen_problem.settings.final_time,
+                unit_scale=True,
+            )
+
+        while self.hydrogen_problem.t.value < self.hydrogen_problem.settings.final_time:
+            self.iterate()
+
+        if self.hydrogen_problem.show_progress_bar:
+            self.hydrogen_problem.progress_bar.refresh()
+            self.hydrogen_problem.progress_bar.close()

src/festim/heat_transfer_problem.py

@@ -69,7 +69,7 @@ def initialise(self):
         if self.settings.transient:
             # TODO should raise error if no stepsize is provided
             # TODO Should this be an attribute of festim.Stepsize?
-            self.dt = helpers.as_fenics_constant(
+            self._dt = helpers.as_fenics_constant(
                 self.settings.stepsize.initial_value, self.mesh.mesh
             )
 

src/festim/helpers.py

@@ -287,3 +287,25 @@ def update(self, t: float):
     get_interpolation_points = lambda element: element.interpolation_points
 else:
     get_interpolation_points = lambda element: element.interpolation_points()
+
+
+def nmm_interpolate(f_out: fem.function, f_in: fem.function):
+    """Non Matching Mesh Interpolate: interpolate one function (f_in) from one mesh into
+    another function (f_out) with a mismatching mesh
+
+    args:
+        f_out: function to interpolate into
+        f_in: function to interpolate from
+
+    notes:
+    https://fenicsproject.discourse.group/t/gjk-error-in-interpolation-between-non-matching-second-ordered-3d-meshes/16086/6
+    """
+
+    dim = f_out.function_space.mesh.topology.dim
+    index_map = f_out.function_space.mesh.topology.index_map(dim)
+    ncells = index_map.size_local + index_map.num_ghosts
+    cells = np.arange(ncells, dtype=np.int32)
+    interpolation_data = fem.create_interpolation_data(
+        f_out.function_space, f_in.function_space, cells, padding=1e-11
+    )
+    f_out.interpolate_nonmatching(f_in, cells, interpolation_data=interpolation_data)

src/festim/hydrogen_transport_problem.py

@@ -208,7 +208,7 @@ def temperature_fenics(self, value):
     def temperature_time_dependent(self):
         if self.temperature is None:
             return False
-        if isinstance(self.temperature, fem.Constant):
+        if isinstance(self.temperature, fem.Constant | fem.Function):
             return False
         if callable(self.temperature):
             arguments = self.temperature.__code__.co_varnames
@@ -271,7 +271,7 @@ def initialise(self):
         if self.settings.transient:
             # TODO should raise error if no stepsize is provided
             # TODO Should this be an attribute of festim.Stepsize?
-            self.dt = as_fenics_constant(
+            self._dt = as_fenics_constant(
                 self.settings.stepsize.initial_value, self.mesh.mesh
             )
 
@@ -765,28 +765,28 @@ def update_time_dependent_values(self):
                 if isinstance(reactant, _species.ImplicitSpecies):
                     reactant.update_density(t=t)
 
-        if not self.temperature_time_dependent:
-            return
-
-        if isinstance(self.temperature_fenics, fem.Constant):
-            self.temperature_fenics.value = self.temperature(t=t)
-        elif isinstance(self.temperature_fenics, fem.Function):
-            self.temperature_fenics.interpolate(self.temperature_expr)
+        if (
+            isinstance(self.temperature, fem.Function)
+            or self.temperature_time_dependent
+        ):
+            for bc in self.boundary_conditions:
+                if isinstance(
+                    bc,
+                    boundary_conditions.FixedConcentrationBC
+                    | boundary_conditions.ParticleFluxBC,
+                ):
+                    if bc.temperature_dependent:
+                        bc.update(t=t)
 
-        for bc in self.boundary_conditions:
-            if isinstance(
-                bc,
-                (
-                    boundary_conditions.FixedConcentrationBC,
-                    boundary_conditions.ParticleFluxBC,
-                ),
-            ):
-                if bc.temperature_dependent:
-                    bc.update(t=t)
+            for source in self.sources:
+                if source.value.temperature_dependent:
+                    source.value.update(t=t)
 
-        for source in self.sources:
-            if source.value.temperature_dependent:
-                source.value.update(t=t)
+        if self.temperature_time_dependent:
+            if isinstance(self.temperature_fenics, fem.Constant):
+                self.temperature_fenics.value = self.temperature(t=t)
+            elif isinstance(self.temperature_fenics, fem.Function):
+                self.temperature_fenics.interpolate(self.temperature_expr)
 
     def post_processing(self):
         """Post processes the model"""
@@ -928,7 +928,7 @@ def initialise(self):
         if self.settings.transient:
             # TODO should raise error if no stepsize is provided
             # TODO Should this be an attribute of festim.Stepsize?
-            self.dt = as_fenics_constant(
+            self._dt = as_fenics_constant(
                 self.settings.stepsize.initial_value, self.mesh.mesh
             )
 

src/festim/problem.py

@@ -71,6 +71,10 @@ def volume_subdomains(self):
     def surface_subdomains(self):
         return [s for s in self.subdomains if isinstance(s, F.SurfaceSubdomain)]
 
+    @property
+    def dt(self):
+        return self._dt
+
     def define_meshtags_and_measures(self):
         """Defines the facet and volume meshtags of the model which are used
         to define the measures fo the model, dx and ds"""

src/festim/problem_change_of_var.py

@@ -1,13 +1,14 @@
-import festim.boundary_conditions
-from festim.hydrogen_transport_problem import HydrogenTransportProblem
-from festim import boundary_conditions, as_fenics_constant
-import festim
-from festim.helpers import get_interpolation_points
-import festim.species as _species
+from typing import List
+
 import ufl
 from dolfinx import fem
 
-from typing import List
+import festim
+import festim.boundary_conditions
+import festim.species as _species
+from festim import boundary_conditions
+from festim.helpers import as_fenics_constant, get_interpolation_points
+from festim.hydrogen_transport_problem import HydrogenTransportProblem
 
 
 class HydrogenTransportProblemDiscontinuousChangeVar(HydrogenTransportProblem):
@@ -23,7 +24,7 @@ def initialise(self):
         if self.settings.transient:
             # TODO should raise error if no stepsize is provided
             # TODO Should this be an attribute of festim.Stepsize?
-            self.dt = as_fenics_constant(
+            self._dt = as_fenics_constant(
                 self.settings.stepsize.initial_value, self.mesh.mesh
             )
 
@@ -79,6 +80,7 @@ def create_formulation(self):
 
             # hack enforce the concentration attribute of the species for all species
             # to be used in reaction.reaction_term
+
             for spe in self.species:
                 if spe.mobile:
                     K_S = reaction.volume.material.get_solubility_coefficient(