Merge pull request #904 from RemDelaporteMathurin/change-of-var-method

Change of variable method for interface discontinuities

Author: RemDelaporteMathurin

pyproject.toml

@@ -13,7 +13,7 @@ description = "Finite element simulations of hydrogen transport"
 readme = "README/md"
 requires-python = ">=3.9"
 license = { file = "LICENSE" }
-dependencies = ["tqdm", "scifem>=0.2.8"]
+dependencies = ["tqdm", "scifem>=0.2.13"]
 classifiers = [
     "Natural Language :: English",
     "Topic :: Scientific/Engineering",

src/festim/__init__.py

@@ -42,6 +42,7 @@
     HydrogenTransportProblem,
     HTransportProblemPenalty,
 )
+from .problem_change_of_var import HydrogenTransportProblemDiscontinuousChangeVar
 from .initial_condition import InitialCondition, InitialTemperature
 from .material import Material
 from .mesh.mesh import Mesh
@@ -51,7 +52,7 @@
 from .reaction import Reaction
 from .settings import Settings
 from .source import HeatSource, ParticleSource, SourceBase
-from .species import ImplicitSpecies, Species, find_species_from_name
+from .species import ImplicitSpecies, Species, find_species_from_name, SpeciesChangeVar
 from .stepsize import Stepsize
 from .subdomain.interface import Interface
 from .subdomain.surface_subdomain import SurfaceSubdomain, find_surface_from_id

src/festim/boundary_conditions/dirichlet_bc.py

@@ -5,9 +5,12 @@
 import ufl
 from dolfinx import fem
 from dolfinx import mesh as _mesh
+import ufl.core
+import ufl.core.expr
 
 from festim import helpers
 from festim import subdomain as _subdomain
+from festim.species import Species
 
 
 class DirichletBCBase:
@@ -112,7 +115,7 @@ def define_surface_subdomain_dofs(
             )
         if mesh.topology.dim - 1 != facet_meshtags.dim:
             raise ValueError(
-                f"Meshtags of dimension {facet_meshtags.dim}, expected {mesh.topology.dim-1}"
+                f"Meshtags of dimension {facet_meshtags.dim}, expected {mesh.topology.dim - 1}"
             )
         bc_dofs = fem.locate_dofs_topological(
             function_space, facet_meshtags.dim, facet_meshtags.find(self.subdomain.id)
@@ -124,14 +127,16 @@ def update(self, t: float):
         """Updates the boundary condition value
 
         Args:
-            t (float): the time
+            t: the time
         """
         if callable(self.value):
             arguments = self.value.__code__.co_varnames
             if isinstance(self.value_fenics, fem.Constant) and "t" in arguments:
                 self.value_fenics.value = self.value(t=t)
             else:
                 self.value_fenics.interpolate(self.bc_expr)
+        elif self.bc_expr is not None:
+            self.value_fenics.interpolate(self.bc_expr)
 
 
 class FixedConcentrationBC(DirichletBCBase):
@@ -163,13 +168,13 @@ class FixedConcentrationBC(DirichletBCBase):
 
     """
 
-    species: str
+    species: Species
 
     def __init__(
         self,
         subdomain: _subdomain.SurfaceSubdomain,
         value: np.ndarray | fem.Constant | int | float | Callable,
-        species: str,
+        species: Species,
     ):
         self.species = species
         super().__init__(subdomain, value)
@@ -191,6 +196,7 @@ def create_value(
         function_space: fem.FunctionSpace,
         temperature: float | fem.Constant,
         t: float | fem.Constant,
+        K_S: fem.Function = None,
     ):
         """Creates the value of the boundary condition as a fenics object and sets it to
         self.value_fenics.
@@ -200,9 +206,11 @@ def create_value(
         expression of the function is stored in `bc_expr`.
 
         Args:
-            function_space (dolfinx.fem.FunctionSpace): the function space
+            function_space: the function space
             temperature: The temperature
-            t (dolfinx.fem.Constant): the time
+            t: the time
+            K_S: The solubility of the species. If provided, the value of the boundary condition
+                is divided by K_S (change of variable method).
         """
         mesh = function_space.mesh
         x = ufl.SpatialCoordinate(mesh)
@@ -241,6 +249,54 @@ def create_value(
                 )
                 self.value_fenics.interpolate(self.bc_expr)
 
+        # if K_S is provided, divide the value by K_S (change of variable method)
+        if K_S is not None:
+            if isinstance(self.value, (int, float)):
+                val_as_cst = helpers.as_fenics_constant(mesh=mesh, value=self.value)
+                self.bc_expr = fem.Expression(
+                    val_as_cst / K_S,
+                    function_space.element.interpolation_points(),
+                )
+                self.value_fenics = fem.Function(function_space)
+                self.value_fenics.interpolate(self.bc_expr)
+
+            elif callable(self.value):
+                arguments = self.value.__code__.co_varnames
+
+                if "t" in arguments and "x" not in arguments and "T" not in arguments:
+                    # only t is an argument
+
+                    # check that value returns a ufl expression
+                    if not isinstance(self.value(t=t), (ufl.core.expr.Expr)):
+                        raise ValueError(
+                            "self.value should return a ufl expression"
+                            + f"{type(self.value(t=t))} "
+                        )
+
+                    self.bc_expr = fem.Expression(
+                        self.value(t=t) / K_S,
+                        function_space.element.interpolation_points(),
+                    )
+                    self.value_fenics = fem.Function(function_space)
+                    self.value_fenics.interpolate(self.bc_expr)
+                else:
+                    self.value_fenics = fem.Function(function_space)
+                    kwargs = {}
+                    if "t" in arguments:
+                        kwargs["t"] = t
+                    if "x" in arguments:
+                        kwargs["x"] = x
+                    if "T" in arguments:
+                        kwargs["T"] = temperature
+
+                    # store the expression of the boundary condition
+                    # to update the value_fenics later
+                    self.bc_expr = fem.Expression(
+                        self.value(**kwargs) / K_S,
+                        function_space.element.interpolation_points(),
+                    )
+                    self.value_fenics.interpolate(self.bc_expr)
+
 
 # alias for FixedConcentrationBC
 DirichletBC = FixedConcentrationBC

src/festim/hydrogen_transport_problem.py

@@ -8,7 +8,7 @@
 import tqdm.autonotebook
 import ufl
 from dolfinx import fem
-from scifem import NewtonSolver
+from scifem import BlockedNewtonSolver
 
 import festim.boundary_conditions
 import festim.problem
@@ -167,6 +167,9 @@ def __init__(
         self.temperature_fenics = None
         self._vtxfiles: list[dolfinx.io.VTXWriter] = []
 
+        self._element_for_traps = "DG"
+        self.petcs_options = petsc_options
+
     @property
     def temperature(self):
         return self._temperature
@@ -463,14 +466,25 @@ def define_function_spaces(self):
             degree,
             basix.LagrangeVariant.equispaced,
         )
+        element_DG = basix.ufl.element(
+            "DG",
+            self.mesh.mesh.basix_cell(),
+            degree,
+            basix.LagrangeVariant.equispaced,
+        )
 
         if not self.multispecies:
             element = element_CG
         else:
             elements = []
             for spe in self.species:
                 if isinstance(spe, _species.Species):
-                    elements.append(element_CG)
+                    if spe.mobile:
+                        elements.append(element_CG)
+                    elif self._element_for_traps == "DG":
+                        elements.append(element_DG)
+                    else:
+                        elements.append(element_CG)
             element = basix.ufl.mixed_element(elements)
 
         self.function_space = fem.functionspace(self.mesh.mesh, element)
@@ -730,7 +744,8 @@ def create_formulation(self):
                         # check reactions
                         for reaction in self.reactions:
                             if vol == reaction.volume:
-                                not_defined_in_volume.remove(vol)
+                                if vol in not_defined_in_volume:
+                                    not_defined_in_volume.remove(vol)
 
                     # add c = 0 to formulation where needed
                     for vol in not_defined_in_volume:
@@ -1213,18 +1228,32 @@ def mixed_term(u, v, n):
         self.forms = dolfinx.fem.form(
             [subdomain.F for subdomain in self.volume_subdomains],
             entity_maps=entity_maps,
+            jit_options={
+                "cffi_extra_compile_args": ["-O3", "-march=native"],
+                "cffi_libraries": ["m"],
+            },
+        )
+        self.J = dolfinx.fem.form(
+            J,
+            entity_maps=entity_maps,
+            jit_options={
+                "cffi_extra_compile_args": ["-O3", "-march=native"],
+                "cffi_libraries": ["m"],
+            },
         )
-        self.J = dolfinx.fem.form(J, entity_maps=entity_maps)
 
     def create_solver(self):
-        self.solver = NewtonSolver(
+        self.solver = BlockedNewtonSolver(
             self.forms,
-            self.J,
             [subdomain.u for subdomain in self.volume_subdomains],
+            J=self.J,
             bcs=self.bc_forms,
-            max_iterations=self.settings.max_iterations,
             petsc_options=self.petsc_options,
         )
+        self.solver.max_iterations = self.settings.max_iterations
+        self.solver.convergence_criterion = self.settings.convergence_criterion
+        self.solver.atol = self.settings.atol
+        self.solver.rtol = self.settings.rtol
 
     def create_flux_values_fenics(self):
         """For each particle flux create the ``value_fenics`` attribute"""
@@ -1283,8 +1312,8 @@ def iterate(self):
 
         self.update_time_dependent_values()
 
-        # solve main problem
-        self.solver.solve(self.settings.atol, self.settings.rtol)
+        # Solve main problem
+        self.solver.solve()
 
         # post processing
         self.post_processing()
@@ -1312,7 +1341,7 @@ def run(self):
                 self.progress_bar.refresh()  # refresh progress bar to show 100%
         else:
             # Solve steady-state
-            self.solver.solve(self.settings.rtol)
+            self.solver.solve()
             self.post_processing()
 
     def __del__(self):

src/festim/material.py

@@ -121,6 +121,58 @@ def get_E_D(self, species=None):
         else:
             raise TypeError("E_D must be either a float, int or a dict")
 
+    def get_K_S_0(self, species=None) -> float:
+        """Returns the pre-exponential factor of the solubility coefficient
+        Args:
+            species: the species we want the pre-exponential
+                factor of the solubility coefficient of. Only needed if K_S_0 is a dict.
+        Returns:
+            the pre-exponential factor of the solubility coefficient
+        """
+
+        if isinstance(self.K_S_0, (float, int)):
+            return self.K_S_0
+
+        elif isinstance(self.K_S_0, dict):
+            if species is None:
+                raise ValueError("species must be provided if K_S_0 is a dict")
+
+            if species in self.K_S_0:
+                return self.K_S_0[species]
+            elif species.name in self.K_S_0:
+                return self.K_S_0[species.name]
+            else:
+                raise ValueError(f"{species} is not in K_S_0 keys")
+
+        else:
+            raise TypeError("K_S_0 must be either a float, int or a dict")
+
+    def get_E_K_S(self, species=None) -> float:
+        """Returns the activation energy of the solubility coefficient
+        Args:
+            species: the species we want the activation
+                energy of the solubility coefficient of. Only needed if E_K_S is a dict.
+        Returns:
+            the activation energy of the solubility coefficient
+        """
+
+        if isinstance(self.E_K_S, (float, int)):
+            return self.E_K_S
+
+        elif isinstance(self.E_K_S, dict):
+            if species is None:
+                raise ValueError("species must be provided if E_K_S is a dict")
+
+            if species in self.E_K_S:
+                return self.E_K_S[species]
+            elif species.name in self.E_K_S:
+                return self.E_K_S[species.name]
+            else:
+                raise ValueError(f"{species} is not in E_K_S keys")
+
+        else:
+            raise TypeError("E_K_S must be either a float, int or a dict")
+
     def get_diffusion_coefficient(self, mesh, temperature, species=None):
         """Defines the diffusion coefficient
         Args:
@@ -184,7 +236,7 @@ def get_solubility_coefficient(self, mesh, temperature, species=None):
         Returns:
             ufl.algebra.Product: the solubility coefficient
         """
-        # TODO use get_D_0 and get_E_D to refactore this method, something like:
+        # TODO use get_K_S0 and get_E_K_S to refactore this method, something like:
         # K_S_0 = self.get_K_S_0(species=species)
         # E_K_S = self.get_E_K_S(species=species)