diff --git a/src/festim/__init__.py b/src/festim/__init__.py
index ecbcee40d..a3e85a9ca 100644
--- a/src/festim/__init__.py
+++ b/src/festim/__init__.py
@@ -44,22 +44,22 @@
Value,
get_interpolation_points,
)
+from .problem import ProblemBase
from .hydrogen_transport_problem import (
HTransportProblemDiscontinuous,
HTransportProblemPenalty,
HydrogenTransportProblem,
+ 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, SpeciesChangeVar, find_species_from_name
+from .species import ImplicitSpecies, Species, find_species_from_name
from .stepsize import Stepsize
from .subdomain.interface import Interface
from .subdomain.surface_subdomain import SurfaceSubdomain, find_surface_from_id
diff --git a/src/festim/coupled_heat_hydrogen_problem.py b/src/festim/coupled_heat_hydrogen_problem.py
index e288fb217..325c122b1 100644
--- a/src/festim/coupled_heat_hydrogen_problem.py
+++ b/src/festim/coupled_heat_hydrogen_problem.py
@@ -7,8 +7,8 @@
HTransportProblemDiscontinuous,
HTransportProblemPenalty,
HydrogenTransportProblem,
+ HydrogenTransportProblemDiscontinuousChangeVar,
)
-from festim.problem_change_of_var import HydrogenTransportProblemDiscontinuousChangeVar
class CoupledTransientHeatTransferHydrogenTransport:
diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 2b7326997..59c3dfd60 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -1,4 +1,5 @@
from collections.abc import Callable
+from typing import List
from mpi4py import MPI
@@ -1471,3 +1472,259 @@ def create_formulation(self):
entity_maps=entity_maps,
)
self.J = dolfinx.fem.form(J, entity_maps=entity_maps)
+
+
+class HydrogenTransportProblemDiscontinuousChangeVar(HydrogenTransportProblem):
+ species: List[_species.Species]
+
+ def initialise(self):
+ # check if a SurfaceReactionBC is given
+ for bc in self.boundary_conditions:
+ if isinstance(bc, (boundary_conditions.SurfaceReactionBC)):
+ raise ValueError(
+ f"{type(bc)} not implemented for HydrogenTransportProblemDiscontinuousChangeVar"
+ )
+ if isinstance(bc, boundary_conditions.ParticleFluxBC):
+ if bc.species_dependent_value:
+ raise ValueError(
+ f"{type(bc)} concentration-dependent not implemented for HydrogenTransportProblemDiscontinuousChangeVar"
+ )
+
+ super().initialise()
+
+ 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
+ )
+ if spe.mobile:
+ K_S = vol.material.get_solubility_coefficient(
+ self.mesh.mesh, self.temperature_fenics, spe
+ )
+ 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:
+ self.add_reaction_term(reaction)
+
+ # add sources
+ for source in self.sources:
+ self.formulation -= (
+ source.value.fenics_object
+ * 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 (
+ spe in reaction.product
+ ): # TODO we probably need this in HydrogenTransportProblem too no?
+ 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 add_reaction_term(self, reaction: _reaction.Reaction):
+ """Adds the reaction term to the formulation"""
+
+ products = (
+ reaction.product
+ if isinstance(reaction.product, list)
+ else [reaction.product]
+ )
+
+ # we cannot use the `concentration` attribute of the mobile species and need to use u * K_S instead
+
+ def get_concentrations(species_list) -> List:
+ concentrations = []
+ for spe in species_list:
+ if isinstance(spe, _species.ImplicitSpecies):
+ concentrations.append(None)
+ elif spe.mobile:
+ K_S = reaction.volume.material.get_solubility_coefficient(
+ self.mesh.mesh, self.temperature_fenics, spe
+ )
+ concentrations.append(spe.solution * K_S)
+ else:
+ concentrations.append(None)
+ return concentrations
+
+ reactant_concentrations = get_concentrations(reaction.reactant)
+ product_concentrations = get_concentrations(products)
+
+ # get the reaction term from the reaction
+ reaction_term = reaction.reaction_term(
+ temperature=self.temperature_fenics,
+ reactant_concentrations=reactant_concentrations,
+ product_concentrations=product_concentrations,
+ )
+
+ # add reaction term to formulation
+ # reactant
+ for reactant in reaction.reactant:
+ if isinstance(reactant, festim.species.Species):
+ self.formulation += (
+ reaction_term * reactant.test_function * self.dx(reaction.volume.id)
+ )
+
+ # product
+ for product in products:
+ self.formulation += (
+ -reaction_term * product.test_function * self.dx(reaction.volume.id)
+ )
+
+ def initialise_exports(self):
+ self.override_post_processing_solution()
+ super().initialise_exports()
+
+ 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))
+
+ for spe in self.species:
+ if not spe.mobile:
+ continue
+ K_S0 = fem.Function(Q0)
+ E_KS = fem.Function(Q0)
+ for subdomain in self.volume_subdomains:
+ entities = subdomain.locate_subdomain_entities_correct(
+ self.volume_meshtags
+ )
+ K_S0.x.array[entities] = subdomain.material.get_K_S_0(spe)
+ E_KS.x.array[entities] = subdomain.material.get_E_K_S(spe)
+
+ K_S = K_S0 * ufl.exp(-E_KS / (festim.k_B * self.temperature_fenics))
+
+ theta = spe.solution
+
+ spe.dg_expr = fem.Expression(
+ theta * K_S, get_interpolation_points(Q1.element)
+ )
+ spe.post_processing_solution = fem.Function(Q1)
+ spe.post_processing_solution.interpolate(
+ spe.dg_expr
+ ) # NOTE: do we need this line since it's in initialise?
+
+ def post_processing(self):
+ # need to compute c = theta * K_S
+ # this expression is stored in species.dg_expr
+ for spe in self.species:
+ if not spe.mobile:
+ continue
+ spe.post_processing_solution.interpolate(spe.dg_expr)
+
+ super().post_processing()
+
+ def create_dirichletbc_form(self, bc: festim.FixedConcentrationBC):
+ """Creates a dirichlet boundary condition form
+
+ Args:
+ bc (festim.DirichletBC): the boundary condition
+
+ Returns:
+ dolfinx.fem.bcs.DirichletBC: A representation of
+ the boundary condition for modifying linear systems.
+ """
+ # create value_fenics
+ if not self.multispecies:
+ function_space_value = bc.species.sub_function_space
+ else:
+ function_space_value = bc.species.collapsed_function_space
+
+ # create K_S function
+ Q0 = fem.functionspace(self.mesh.mesh, ("DG", 0))
+ K_S0 = fem.Function(Q0)
+ E_KS = fem.Function(Q0)
+ for subdomain in self.volume_subdomains:
+ entities = subdomain.locate_subdomain_entities_correct(self.volume_meshtags)
+ K_S0.x.array[entities] = subdomain.material.get_K_S_0(bc.species)
+ E_KS.x.array[entities] = subdomain.material.get_E_K_S(bc.species)
+
+ K_S = K_S0 * ufl.exp(-E_KS / (festim.k_B * self.temperature_fenics))
+
+ bc.create_value(
+ temperature=self.temperature_fenics,
+ function_space=function_space_value,
+ t=self.t,
+ K_S=K_S,
+ )
+
+ # get dofs
+ if self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
+ function_space_dofs = (
+ bc.species.sub_function_space,
+ bc.species.collapsed_function_space,
+ )
+ else:
+ function_space_dofs = bc.species.sub_function_space
+
+ bc_dofs = bc.define_surface_subdomain_dofs(
+ facet_meshtags=self.facet_meshtags,
+ function_space=function_space_dofs,
+ )
+
+ # create form
+ if not self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
+ # no need to pass the functionspace since value_fenics is already a Function
+ function_space_form = None
+ else:
+ function_space_form = bc.species.sub_function_space
+
+ form = fem.dirichletbc(
+ value=bc.value_fenics,
+ dofs=bc_dofs,
+ V=function_space_form,
+ )
+
+ return form
+
+ def update_time_dependent_values(self):
+ super().update_time_dependent_values()
+
+ if self.temperature_time_dependent:
+ for bc in self.boundary_conditions:
+ if isinstance(bc, boundary_conditions.FixedConcentrationBC):
+ bc.update(self.t)
diff --git a/src/festim/problem_change_of_var.py b/src/festim/problem_change_of_var.py
deleted file mode 100644
index 65ee6db74..000000000
--- a/src/festim/problem_change_of_var.py
+++ /dev/null
@@ -1,254 +0,0 @@
-from typing import List
-
-import ufl
-from dolfinx import fem
-
-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):
- 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.create_implicit_species_value_fenics()
-
- self.define_temperature()
- self.define_boundary_conditions()
- self.convert_source_input_values_to_fenics_objects()
- self.create_flux_values_fenics()
- self.create_initial_conditions()
- self.create_formulation()
- self.create_solver()
- self.override_post_processing_solution() # NOTE this is the only difference with parent class
- 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
- )
- if spe.mobile:
- K_S = vol.material.get_solubility_coefficient(
- self.mesh.mesh, self.temperature_fenics, spe
- )
- 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:
- if spe.mobile:
- K_S = reaction.volume.material.get_solubility_coefficient(
- self.mesh.mesh, self.temperature_fenics, spe
- )
- 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_object
- * 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 (
- spe in reaction.product
- ): # TODO we probably need this in HydrogenTransportProblem too no?
- 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))
-
- for spe in self.species:
- if not spe.mobile:
- continue
- K_S0 = fem.Function(Q0)
- E_KS = fem.Function(Q0)
- for subdomain in self.volume_subdomains:
- entities = subdomain.locate_subdomain_entities_correct(
- self.volume_meshtags
- )
- K_S0.x.array[entities] = subdomain.material.get_K_S_0(spe)
- E_KS.x.array[entities] = subdomain.material.get_E_K_S(spe)
-
- K_S = K_S0 * ufl.exp(-E_KS / (festim.k_B * self.temperature_fenics))
-
- theta = spe.solution
-
- spe.dg_expr = fem.Expression(
- theta * K_S, get_interpolation_points(Q1.element)
- )
- spe.post_processing_solution = fem.Function(Q1)
- spe.post_processing_solution.interpolate(
- spe.dg_expr
- ) # NOTE: do we need this line since it's in initialise?
-
- def post_processing(self):
- # need to compute c = theta * K_S
- # this expression is stored in species.dg_expr
- for spe in self.species:
- if not spe.mobile:
- continue
- spe.post_processing_solution.interpolate(spe.dg_expr)
-
- super().post_processing()
-
- def create_dirichletbc_form(self, bc: festim.FixedConcentrationBC):
- """Creates a dirichlet boundary condition form
-
- Args:
- bc (festim.DirichletBC): the boundary condition
-
- Returns:
- dolfinx.fem.bcs.DirichletBC: A representation of
- the boundary condition for modifying linear systems.
- """
- # create value_fenics
- if not self.multispecies:
- function_space_value = bc.species.sub_function_space
- else:
- function_space_value = bc.species.collapsed_function_space
-
- # create K_S function
- Q0 = fem.functionspace(self.mesh.mesh, ("DG", 0))
- K_S0 = fem.Function(Q0)
- E_KS = fem.Function(Q0)
- for subdomain in self.volume_subdomains:
- entities = subdomain.locate_subdomain_entities_correct(self.volume_meshtags)
- K_S0.x.array[entities] = subdomain.material.get_K_S_0(bc.species)
- E_KS.x.array[entities] = subdomain.material.get_E_K_S(bc.species)
-
- K_S = K_S0 * ufl.exp(-E_KS / (festim.k_B * self.temperature_fenics))
-
- bc.create_value(
- temperature=self.temperature_fenics,
- function_space=function_space_value,
- t=self.t,
- K_S=K_S,
- )
-
- # get dofs
- if self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
- function_space_dofs = (
- bc.species.sub_function_space,
- bc.species.collapsed_function_space,
- )
- else:
- function_space_dofs = bc.species.sub_function_space
-
- bc_dofs = bc.define_surface_subdomain_dofs(
- facet_meshtags=self.facet_meshtags,
- function_space=function_space_dofs,
- )
-
- # create form
- if not self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
- # no need to pass the functionspace since value_fenics is already a Function
- function_space_form = None
- else:
- function_space_form = bc.species.sub_function_space
-
- form = fem.dirichletbc(
- value=bc.value_fenics,
- dofs=bc_dofs,
- V=function_space_form,
- )
-
- return form
-
- def update_time_dependent_values(self):
- super().update_time_dependent_values()
-
- if self.temperature_time_dependent:
- for bc in self.boundary_conditions:
- if isinstance(bc, boundary_conditions.FixedConcentrationBC):
- bc.update(self.t)
diff --git a/src/festim/reaction.py b/src/festim/reaction.py
index 64a4effea..3b485d99e 100644
--- a/src/festim/reaction.py
+++ b/src/festim/reaction.py
@@ -1,6 +1,7 @@
-from typing import Optional, Union
+from typing import Optional, Union, List
from ufl import exp
+from ufl.core.expr import Expr
from festim import k_B as _k_B
from festim.species import ImplicitSpecies as _ImplicitSpecies
@@ -102,11 +103,27 @@ def __str__(self) -> str:
products = self.product
return f"{reactants} <--> {products}"
- def reaction_term(self, temperature):
+ def reaction_term(
+ self,
+ temperature,
+ reactant_concentrations: List = None,
+ product_concentrations: List = None,
+ ) -> Expr:
"""Compute the reaction term at a given temperature.
Arguments:
- temperature (): The temperature at which the reaction term is computed.
+ temperature: The temperature at which the reaction term is computed.
+ reactant_concentrations: The concentrations of the reactants. Must
+ be the same length as the reactants. If None, the ``concentration``
+ attribute of each reactant is used. If an element is None, the
+ ``concentration`` attribute of the reactant is used.
+ product_concentrations: The concentrations of the products. Must
+ be the same length as the products. If None, the ``concentration``
+ attribute of each product is used. If an element is None, the
+ ``concentration`` attribute of the product is used.
+
+ Returns:
+ The reaction term to be used in a formulation.
"""
if self.product == []:
@@ -130,6 +147,9 @@ def reaction_term(self, temperature):
"E_p cannot be None when reaction products are present."
)
+ products = self.product if isinstance(self.product, list) else [self.product]
+
+ # reaction rates
k = self.k_0 * exp(-self.E_k / (_k_B * temperature))
if self.p_0 and self.E_p:
@@ -139,20 +159,35 @@ def reaction_term(self, temperature):
else:
p = 0
+ # if reactant_concentrations is provided, use these concentrations
reactants = self.reactant
- product_of_reactants = reactants[0].concentration
- for reactant in reactants[1:]:
- product_of_reactants *= reactant.concentration
-
- if isinstance(self.product, list):
- products = self.product
+ if reactant_concentrations is not None:
+ assert len(reactant_concentrations) == len(reactants)
+ for i, reactant in enumerate(reactants):
+ if reactant_concentrations[i] is None:
+ reactant_concentrations[i] = reactant.concentration
+ else:
+ reactant_concentrations = [reactant.concentration for reactant in reactants]
+
+ # if product_concentrations is provided, use these concentrations
+ if product_concentrations is not None:
+ assert len(product_concentrations) == len(products)
+ for i, product in enumerate(products):
+ if product_concentrations[i] is None:
+ product_concentrations[i] = product.concentration
else:
- products = [self.product]
+ product_concentrations = [product.concentration for product in products]
- if len(products) > 0:
- product_of_products = products[0].solution
- for product in products[1:]:
- product_of_products *= product.solution
+ # multiply all concentrations to be used in the term
+ product_of_reactants = reactant_concentrations[0]
+ for reactant_conc in reactant_concentrations[1:]:
+ product_of_reactants *= reactant_conc
+
+ if products:
+ product_of_products = product_concentrations[0]
+ for product_conc in product_concentrations[1:]:
+ product_of_products *= product_conc
else:
product_of_products = 0
+
return k * product_of_reactants - (p * product_of_products)
diff --git a/src/festim/species.py b/src/festim/species.py
index 437e109d8..1a9e6b373 100644
--- a/src/festim/species.py
+++ b/src/festim/species.py
@@ -225,13 +225,3 @@ def find_species_from_name(name: str, species: list):
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
diff --git a/test/system_tests/test_multi_material_change_of_var.py b/test/system_tests/test_multi_material_change_of_var.py
index b41d2ded9..efc7137e2 100644
--- a/test/system_tests/test_multi_material_change_of_var.py
+++ b/test/system_tests/test_multi_material_change_of_var.py
@@ -103,8 +103,8 @@ def test_run():
right_surface,
]
- H = F.SpeciesChangeVar("H", mobile=True)
- trapped_H = F.SpeciesChangeVar("H_trapped", mobile=False)
+ H = F.Species("H", mobile=True)
+ trapped_H = F.Species("H_trapped", mobile=False)
empty_trap = F.ImplicitSpecies(n=0.5, others=[trapped_H])
my_model.species = [H, trapped_H]
@@ -193,7 +193,7 @@ def c_exact_bot_np(x):
bottom_surface,
]
- H = F.SpeciesChangeVar("H", mobile=True)
+ H = F.Species("H", mobile=True)
my_model.species = [H]