Value class to allow species-dependent inputs

src/festim/boundary_conditions/flux_bc.py

@@ -151,10 +151,12 @@ class ParticleFluxBC(FluxBCBase):
             is applied
         value (float, fem.Constant, callable): the value of the particle flux
         species (festim.Species): the species to which the flux is applied
-        species_dependent_value (dict): a dictionary containing the species
-        that the value. Example: {"name": species}
-            where "name" is the variable name in the callable value and species
-            is a festim.Species object.
+        species_dependent_value (dict): a dictionary mapping the argument names in a
+            callable value to festim.Species objects, allowing the flux to depend on the
+            concentration of other species. Example: {"c1": species1} where "c1" is the
+            argument name in the callable value and species1 is a festim.Species
+            object. Ignored if value is not callable. Defaults to an empty dict. Not
+            supported by festim.HydrogenTransportProblemDiscontinuousChangeVar.
 
     Attributes:
         subdomain (festim.SurfaceSubdomain): the surface subdomain where the
@@ -166,10 +168,8 @@ class ParticleFluxBC(FluxBCBase):
             fenics format
         bc_expr (fem.Expression): the expression of the particle flux that is used to
             update the value_fenics
-        species_dependent_value (dict): a dictionary containing the species
-        that the value. Example: {"name": species}
-            where "name" is the variable name in the callable value and species
-            is a festim.Species object.
+        species_dependent_value (dict): a dictionary mapping the argument names in a
+            callable value to festim.Species objects
 
 
     Examples:
@@ -193,10 +193,10 @@ class ParticleFluxBC(FluxBCBase):
             species="H", species_dependent_value={"c1": species1})
     """
 
-    def __init__(self, subdomain, value, species, species_dependent_value={}):
+    def __init__(self, subdomain, value, species, species_dependent_value=None):
         super().__init__(subdomain=subdomain, value=value)
         self.species = species
-        self.species_dependent_value = species_dependent_value
+        self.species_dependent_value = species_dependent_value or {}
         self._volume_subdomain = None
 
     def create_value_fenics(self, mesh, temperature, t: fem.Constant):
@@ -218,7 +218,12 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant):
         elif callable(self.value):
             arguments = self.value.__code__.co_varnames
 
-            if "t" in arguments and "x" not in arguments and "T" not in arguments:
+            if (
+                "t" in arguments
+                and "x" not in arguments
+                and "T" not in arguments
+                and not self.species_dependent_value
+            ):
                 # only t is an argument
                 if not isinstance(self.value(t=float(t)), (float, int)):
                     raise ValueError(
@@ -237,9 +242,9 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant):
                     kwargs["T"] = temperature
 
                 for name, species in self.species_dependent_value.items():
-                    if species.concentration:
+                    if species.concentration is not None:
                         kwargs[name] = species.concentration
-                    else:  # probably in discontinuous case
+                    else:  # discontinuous case: one solution per subdomain
                         kwargs[name] = species.subdomain_to_solution[
                             self._volume_subdomain
                         ]

src/festim/helpers.py

@@ -1,4 +1,5 @@
 from collections.abc import Callable
+from typing import TYPE_CHECKING
 
 from mpi4py import MPI
 
@@ -7,6 +8,10 @@
 import ufl
 from dolfinx import fem
 
+if TYPE_CHECKING:
+    from festim.species import Species
+    from festim.subdomain.volume_subdomain import VolumeSubdomain
+
 
 def as_fenics_constant(
     value: float | int | fem.Constant, mesh: dolfinx.mesh.Mesh
@@ -33,21 +38,26 @@ def as_fenics_constant(
         )
 
 
-# TODO change this to accept species dependent values
 def as_mapped_function(
     value: Callable,
     function_space: fem.FunctionSpace | None = None,
     t: fem.Constant | None = None,
     temperature: fem.Function | fem.Constant | ufl.core.expr.Expr | None = None,
+    species_dependent_value: dict[str, "Species"] | None = None,
+    subdomain: "VolumeSubdomain | None" = None,
 ) -> ufl.core.expr.Expr:
-    """Maps a user given callable function to the mesh, time or temperature within
-    festim as needed.
+    """Maps a user given callable function to the mesh, time, temperature or the
+    concentration of other species within festim as needed.
 
     Args:
         value: the callable to convert
         function_space: the function space of the domain, optional
         t: the time, optional
         temperature: the temperature, optional
+        species_dependent_value: a dictionary mapping the argument names in the callable
+            ``value`` to festim.Species objects, optional
+        subdomain: the volume subdomain on which the value is evaluated. Only needed in
+            the discontinuous case to select the correct species solution, optional
 
     Returns:
         The mapped function
@@ -65,32 +75,48 @@ def as_mapped_function(
     if "T" in arguments:
         kwargs["T"] = temperature
 
+    for name, species in (species_dependent_value or {}).items():
+        if species.concentration is not None:
+            kwargs[name] = species.concentration
+        else:  # discontinuous case: the species has one solution per subdomain
+            kwargs[name] = species.subdomain_to_solution[subdomain]
+
     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,
     t: fem.Constant | None = None,
     temperature: fem.Function | fem.Constant | ufl.core.expr.Expr | None = None,
+    species_dependent_value: dict[str, "Species"] | None = None,
+    subdomain: "VolumeSubdomain | None" = None,
 ) -> tuple[fem.Expression, fem.Function]:
-    """Takes a user given callable function, maps the function to the mesh, time or
-    temperature within festim as needed. Then creates the fenics interpolation
-    expression and function objects.
+    """Takes a user given callable function, maps the function to the mesh, time,
+    temperature or the concentration of other species within festim as needed. Then
+    creates the fenics interpolation expression and function objects.
 
     Args:
         value: the callable to convert
         function_space: The function space to interpolate function over
         t: the time, optional
         temperature: the temperature, optional
+        species_dependent_value: a dictionary mapping the argument names in the callable
+            ``value`` to festim.Species objects, optional
+        subdomain: the volume subdomain on which the value is evaluated. Only needed in
+            the discontinuous case to select the correct species solution, optional
 
     Returns:
         fenics interpolation expression, fenics function
     """
 
     mapped_function = as_mapped_function(
-        value=value, function_space=function_space, t=t, temperature=temperature
+        value=value,
+        function_space=function_space,
+        t=t,
+        temperature=temperature,
+        species_dependent_value=species_dependent_value,
+        subdomain=subdomain,
     )
 
     fenics_interpolation_expression = fem.Expression(
@@ -110,15 +136,26 @@ class Value:
 
     Args:
         input_value: The value of the user input
+        species_dependent_value: A dictionary mapping the argument names in a callable
+            ``input_value`` to festim.Species objects. This allows the value to depend
+            on the concentration of other species. Example: ``{"c1": species1}`` where
+            ``"c1"`` is the argument name in the callable ``input_value`` and
+            ``species1`` is a festim.Species object. Ignored if ``input_value`` is not
+            callable. Defaults to None.
 
     Attributes:
         input_value : The value of the user input
+        species_dependent_value : A dictionary mapping the argument names in a callable
+            ``input_value`` to festim.Species objects. An empty dict if the value does
+            not depend on other species
         fenics_interpolation_expression : The expression of the user input that is used
             to update the `fenics_object`
         fenics_object : The value of the user input in fenics format
         explicit_time_dependent : True if the user input value is explicitly time
             dependent
         temperature_dependent : True if the user input value is temperature dependent
+        species_dependent : True if the user input value depends on the concentration of
+            other species
     """
 
     input_value: (
@@ -130,15 +167,20 @@ class Value:
         | ufl.core.expr.Expr
         | fem.Function
     )
+    species_dependent_value: dict[str, "Species"] | None
 
     ufl_expression: ufl.core.expr.Expr
     fenics_interpolation_expression: fem.Expression
     fenics_object: fem.Function | fem.Constant | ufl.core.expr.Expr
     explicit_time_dependent: bool
     temperature_dependent: bool
+    species_dependent: bool
 
-    def __init__(self, input_value):
+    def __init__(
+        self, input_value, species_dependent_value: dict[str, "Species"] | None = None
+    ):
         self.input_value = input_value
+        self.species_dependent_value = species_dependent_value or {}
 
         self.ufl_expression = None
         self.fenics_interpolation_expression = None
@@ -187,6 +229,14 @@ def explicit_time_dependent(self) -> bool:
         else:
             return False
 
+    @property
+    def species_dependent(self) -> bool:
+        """Returns true if the value given depends on the concentration of other
+        species."""
+        if not callable(self.input_value):
+            return False
+        return bool(self.species_dependent_value)
+
     @property
     def temperature_dependent(self) -> bool:
         """Returns true if the value given is temperature dependent."""
@@ -206,6 +256,7 @@ def convert_input_value(
         t: fem.Constant | None = None,
         temperature: fem.Function | fem.Constant | ufl.core.expr.Expr | None = None,
         up_to_ufl_expr: bool | None = False,
+        subdomain: "VolumeSubdomain | None" = None,
     ):
         """Converts a user given value to a relevent fenics object depending on the type
         of the value provided.
@@ -216,6 +267,9 @@ def convert_input_value(
             temperature: the temperature, optional
             up_to_ufl_expr: if True, the value is only mapped to a function if the input
                 is callable, not interpolated or converted to a function, optional
+            subdomain: the volume subdomain on which the value is evaluated. Only needed
+                in the discontinuous case to select the correct species solution,
+                optional
         """
         if isinstance(
             self.input_value, fem.Constant | fem.Function | ufl.core.expr.Expr
@@ -233,7 +287,12 @@ def convert_input_value(
         elif callable(self.input_value):
             args = self.input_value.__code__.co_varnames
             # if only t is an argument, create constant object
-            if "t" in args and "x" not in args and "T" not in args:
+            if (
+                "t" in args
+                and "x" not in args
+                and "T" not in args
+                and not self.species_dependent_value
+            ):
                 if not isinstance(self.input_value(t=float(t)), float | int):
                     raise ValueError(
                         "self.value should return a float or an int, not "
@@ -250,6 +309,8 @@ def convert_input_value(
                     function_space=function_space,
                     t=t,
                     temperature=temperature,
+                    species_dependent_value=self.species_dependent_value,
+                    subdomain=subdomain,
                 )
 
             else:
@@ -259,6 +320,8 @@ def convert_input_value(
                         function_space=function_space,
                         t=t,
                         temperature=temperature,
+                        species_dependent_value=self.species_dependent_value,
+                        subdomain=subdomain,
                     )
                 )
 
@@ -389,3 +452,6 @@ def KSPMonitor(ksp, iter, rnorm):
     dolfinx.log.log(dolfinx.log.LogLevel.DEBUG, f"KSP {iter=}, {_residual0=:.5e}")
     if MPI.COMM_WORLD.rank == 0:
         dolfinx.log.log(dolfinx.log.LogLevel.DEBUG, f"KSP {iter=} {rnorm=:.5e}")
+
+
+lambda c, t: 1 + c + t

src/festim/hydrogen_transport_problem.py

@@ -1453,15 +1453,15 @@ def convert_source_input_values_to_fenics_objects(self):
         for source in self.sources:
             # create value_fenics for all F.ParticleSource objects
             if isinstance(source, _source.ParticleSource):
-                for subdomain in source.species.subdomains:
-                    V = source.species.subdomain_to_function_space[subdomain]
+                V = source.species.subdomain_to_function_space[source.volume]
 
-                    source.value.convert_input_value(
-                        function_space=V,
-                        t=self.t,
-                        temperature=self.temperature_fenics,
-                        up_to_ufl_expr=True,
-                    )
+                source.value.convert_input_value(
+                    function_space=V,
+                    t=self.t,
+                    temperature=self.temperature_fenics,
+                    up_to_ufl_expr=True,
+                    subdomain=source.volume,
+                )
 
     def convert_advection_term_to_fenics_objects(self):
         """For each advection term convert the input value."""
@@ -1996,13 +1996,24 @@ def initialise(self):
                     f"{type(bc)} not implemented for "
                     f"HydrogenTransportProblemDiscontinuousChangeVar"
                 )
+            # with the change of variable, the solution of a mobile species is the
+            # chemical potential and not the concentration, so a species-dependent
+            # value would silently be given the wrong quantity
             if isinstance(bc, boundary_conditions.ParticleFluxBC):
                 if bc.species_dependent_value:
                     raise ValueError(
                         f"{type(bc)} concentration-dependent not implemented for "
                         f"HydrogenTransportProblemDiscontinuousChangeVar"
                     )
 
+        for source in self.sources:
+            if isinstance(source, _source.ParticleSource):
+                if source.value.species_dependent:
+                    raise ValueError(
+                        f"{type(source)} concentration-dependent not implemented for "
+                        f"HydrogenTransportProblemDiscontinuousChangeVar"
+                    )
+
         super().initialise()
 
     def create_formulation(self):