WIP Change of var method for interface discontinuities

src/festim/__init__.py

@@ -38,6 +38,7 @@
     HTransportProblemDiscontinuous,
     HydrogenTransportProblem,
 )
+from .problem_change_of_var import HydrogenTransportProblemDiscontinuousChangeVar
 from .initial_condition import InitialCondition, InitialTemperature
 from .material import Material
 from .mesh.mesh import Mesh
@@ -47,7 +48,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/hydrogen_transport_problem.py

@@ -149,6 +149,7 @@
         settings=None,
         exports=None,
         traps=None,
+        petsc_options=None,
     ):
         super().__init__(
             mesh=mesh,
@@ -167,6 +168,9 @@
         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
@@ -446,14 +450,25 @@
             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)
@@ -690,7 +705,8 @@
                         # 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:

src/festim/problem_change_of_var.py

@@ -0,0 +1,175 @@
+import festim.boundary_conditions
+from festim.hydrogen_transport_problem import HydrogenTransportProblem
+from festim import boundary_conditions, as_fenics_constant
+import festim
+import festim.species as _species
+import festim.exports as exports
+import ufl
+from dolfinx import fem
+
+from typing import List
+
+
+class HydrogenTransportProblemDiscontinuousChangeVar(HydrogenTransportProblem):
+    species: List[_species.Species]
+
+    def initialise(self):
+        self.create_species_from_traps()
+        self.define_function_spaces()
+        self.define_meshtags_and_measures()
+        self.assign_functions_to_species()
+
+        self.t = fem.Constant(self.mesh.mesh, 0.0)
+        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.settings.stepsize.initial_value, self.mesh.mesh
+            )
+
+        self.define_temperature()
+        self.define_boundary_conditions()
+        self.create_source_values_fenics()
+        self.create_flux_values_fenics()
+        self.create_initial_conditions()
+        self.create_formulation()
+        self.create_solver()
+        self.override_post_processing_solution()
+        self.initialise_exports()
+
+    def create_formulation(self):
+        """Creates the formulation of the model"""
+
+        self.formulation = 0
+
+        # add diffusion and time derivative for each species
+        for spe in self.species:
+            u = spe.solution
+            u_n = spe.prev_solution
+            v = spe.test_function
+
+            for vol in self.volume_subdomains:
+                D = vol.material.get_diffusion_coefficient(
+                    self.mesh.mesh, self.temperature_fenics, spe
+                )
+                K_S = vol.material.get_solubility_coefficient(
+                    self.mesh.mesh, self.temperature_fenics, spe
+                )
+                if spe.mobile:
+                    c = u * K_S
+                    c_n = u_n * K_S
+                else:
+                    c = u
+                    c_n = u_n
+                if spe.mobile:
+                    self.formulation += ufl.dot(D * ufl.grad(c), ufl.grad(v)) * self.dx(
+                        vol.id
+                    )
+
+                if self.settings.transient:
+                    self.formulation += ((c - c_n) / self.dt) * v * self.dx(vol.id)
+
+        for reaction in self.reactions:
+            if isinstance(reaction.product, list):
+                products = reaction.product
+            else:
+                products = [reaction.product]
+
+            # hack enforce the concentration attribute of the species for all species
+            # to be used in reaction.reaction_term
+            for spe in self.species:
+                K_S = reaction.volume.material.get_solubility_coefficient(
+                    self.mesh.mesh, self.temperature_fenics, spe
+                )
+                if spe.mobile:
+                    spe.concentration = spe.solution * K_S
+
+            # reactant
+            for reactant in reaction.reactant:
+                if isinstance(reactant, festim.species.Species):
+                    self.formulation += (
+                        reaction.reaction_term(self.temperature_fenics)
+                        * reactant.test_function
+                        * self.dx(reaction.volume.id)
+                    )
+
+            # product
+            for product in products:
+                self.formulation += (
+                    -reaction.reaction_term(self.temperature_fenics)
+                    * product.test_function
+                    * self.dx(reaction.volume.id)
+                )
+        # add sources
+        for source in self.sources:
+            self.formulation -= (
+                source.value_fenics
+                * source.species.test_function
+                * self.dx(source.volume.id)
+            )
+
+        # add fluxes
+        for bc in self.boundary_conditions:
+            if isinstance(bc, boundary_conditions.ParticleFluxBC):
+                self.formulation -= (
+                    bc.value_fenics
+                    * bc.species.test_function
+                    * self.ds(bc.subdomain.id)
+                )
+
+        # check if each species is defined in all volumes
+        if not self.settings.transient:
+            for spe in self.species:
+                # if species mobile, already defined in diffusion term
+                if not spe.mobile:
+                    not_defined_in_volume = self.volume_subdomains.copy()
+                    for vol in self.volume_subdomains:
+                        # check reactions
+                        for reaction in self.reactions:
+                            if vol == reaction.volume:
+                                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:
+                        self.formulation += (
+                            spe.solution * spe.test_function * self.dx(vol.id)
+                        )
+
+    def override_post_processing_solution(self):
+        # override the post-processing solution c = theta * K_S
+        Q0 = fem.functionspace(self.mesh.mesh, ("DG", 0))
+        Q1 = fem.functionspace(self.mesh.mesh, ("DG", 1))
+        K_S0 = fem.Function(Q0)
+        E_KS = fem.Function(Q0)
+        for subdomain in self.volume_subdomains:
+            entities = subdomain.locate_subdomain_entities(self.mesh.mesh)
+            K_S0.x.array[entities] = subdomain.material.K_S_0
+            E_KS.x.array[entities] = subdomain.material.E_K_S
+
+        K_S = K_S0 * ufl.exp(-E_KS / (festim.k_B * self.temperature_fenics))
+
+        for spe in self.species:
+            theta = spe.solution
+
+            spe.dg_expr = fem.Expression(theta * K_S, Q1.element.interpolation_points())
+            spe.post_processing_solution = fem.Function(Q1)
+            spe.post_processing_solution.interpolate(spe.dg_expr)
+
+    def post_processing(self):
+
+        # need to compute c = theta * K_S
+        # this expression is stored in species.dg_expr
+        for spe in self.species:
+            spe.post_processing_solution.interpolate(spe.dg_expr)
+
+        super().post_processing()
+
+    # def define_boundary_conditions(self):
+    #     for bc in self.boundary_conditions:
+    #         if isinstance(bc.species, str):
+    #             # if name of species is given then replace with species object
+    #             bc.species = _species.find_species_from_name(bc.species, self.species)
+
+    #     for bc in self.boundary_conditions:
+    #         pass

src/festim/species.py

@@ -166,3 +166,13 @@
         if spe.name == name:
             return spe
     raise ValueError(f"Species {name} not found in list of species")
+
+
+class SpeciesChangeVar(Species):
+    @property
+    def concentration(self):
+        return self._concentration
+
+    @concentration.setter
+    def concentration(self, value):
+        self._concentration = value

src/festim/subdomain/volume_subdomain.py

@@ -2,6 +2,7 @@
 import numpy as np
 
 from festim.helpers_discontinuity import transfer_meshtags_to_submesh
+from festim.material import Material
 
 
 class VolumeSubdomain:
@@ -24,6 +25,7 @@ class VolumeSubdomain:
     padded: bool
     u: dolfinx.fem.Function
     u_n: dolfinx.fem.Function
+    material: Material
 
     def __init__(self, id, material):
         self.id = id