API reference#

Bases: object

Dirichlet boundary condition class u = value.

Parameters:

subdomain – The surface subdomain where the boundary condition is applied
value – The value of the boundary condition
enforce_weakly – Whether to enforce the boundary condition weakly using Nitsche’s method. Defaults to False.
penalty – The dimensionless penalty parameter to use if enforce_weakly is True. The Nitsche penalty term scales as penalty * D / h, so a value of order 10-100 is appropriate regardless of the material or the mesh size. Defaults to None

Variables:

subdomain (festim.subdomain.surface_subdomain.SurfaceSubdomain) – The surface subdomain where the boundary condition is applied
value (numpy.ndarray | dolfinx.fem.function.Constant | int | float | collections.abc.Callable[[numpy.ndarray], numpy.ndarray] | collections.abc.Callable[[numpy.ndarray, float], numpy.ndarray] | collections.abc.Callable[[float], float]) – The value of the boundary condition
value_fenics (dolfinx.fem.function.Function | dolfinx.fem.function.Constant | numpy.ndarray | float | None) – The value of the boundary condition in fenics format
bc_expr (dolfinx.fem.function.Expression) – The expression of the boundary condition that is used to update the value_fenics
enforce_weakly – Whether to enforce the boundary condition weakly using Nitsche’s method.
penalty – The penalty parameter to use if enforce_weakly is True.

define_surface_subdomain_dofs(facet_meshtags: MeshTags, function_space: FunctionSpace | tuple[FunctionSpace, FunctionSpace]) → ndarray[Any, dtype[int32]] | tuple[ndarray[Any, dtype[int32]], ndarray[Any, dtype[int32]]][source]#

Defines the facets and the degrees of freedom of the boundary condition.

Given the input meshtags, find all facets matching the boundary condition subdomain ID, and locate all DOFs associated with the input function space(s).

Note

For sub-spaces, a tuple of sub-spaces are expected as input, and a tuple of arrays associated to each of the function spaces are returned.

Parameters:

facet_meshtags – MeshTags describing some facets in the domain
mesh
function_space – The function space or a tuple of function spaces:
(sub
collapsed)

numerical_flux(u: Function | Indexed, D: Expr | Function | Constant, mesh) → Expr[source]#

Returns the numerical flux leaving the domain through the surface of the BC.

Nitsche’s method only enforces u = value weakly, so u - value does not vanish on the boundary and the raw gradient -D grad(u).n is not the quantity the discrete scheme conserves. The conserved quantity is this numerical flux, which includes the penalty contribution. Use it (rather than the raw gradient) whenever the flux through a weakly enforced boundary feeds another equation, so that the discrete balance is exact.

Parameters:

u – the solution function associated to the species for which the BC is applied
D – the diffusion coefficient of the species at this surface
mesh – the mesh the surface belongs to

Returns:

the numerical flux, positive when leaving the domain

property time_dependent: bool#: Returns true if the value of the boundary condition is time dependent.

update(t: float)[source]#

Updates the boundary condition value.

Parameters:: t – the time

weak_formulation(u: Function | Indexed, v: Argument | Indexed, ds: Measure, D: Expr | Function | Constant) → Expr[source]#

Returns the symmetric Nitsche weak formulation for the BC This follows the dolfinx tutorial https://jsdokken.com/dolfinx-tutorial/chapter1/nitsche.html

Parameters:

u – the solution function associated to the species for which the BC is applied
v – the test function
ds – the surface measure
D – the diffusion coefficient of the species at this surface. It multiplies the consistency and symmetry terms, without which the scheme is only consistent when D == 1.

Returns:

the weak formulation

Bases: DirichletBCBase

Parameters:

subdomain (festim.Subdomain) – the surface subdomain where the boundary condition is applied
value – The value of the boundary condition. It can be a function of
time (space and/or)
species – The name of the species
enforce_weakly – Whether to enforce the boundary condition weakly using Nitsche’s method. Defaults to False.
penalty – The dimensionless penalty parameter to use if enforce_weakly is True. The Nitsche penalty term scales as penalty * D / h, so a value of order 10-100 is appropriate regardless of the material or the mesh size. Defaults to None

Variables:

temperature_dependent (bool) – True if the value of the bc is
dependent (temperature)

Examples

FixedConcentrationBC(subdomain=my_subdomain, value=1, species="H")
FixedConcentrationBC(subdomain=my_subdomain,
                     value=lambda x: 1 + x[0], species="H")
FixedConcentrationBC(subdomain=my_subdomain,
                     value=lambda t: 1 + t, species="H")
FixedConcentrationBC(subdomain=my_subdomain,
                     value=lambda T: 1 + T, species="H")
FixedConcentrationBC(subdomain=my_subdomain,
                     value=lambda x, t: 1 + x[0] + t, species="H")

create_value(function_space: FunctionSpace, temperature: float | Constant, t: float | Constant, K_S: Function = None)[source]#

Creates the value of the boundary condition as a fenics object and sets it to self.value_fenics. If the value is a constant, it is converted to a dolfinx.fem.Constant. If the value is a function of t, it is converted to dolfinx.fem.Constant. Otherwise, it is converted to a dolfinx.fem.Function.Function and the expression of the function is stored in bc_expr.

Parameters:

function_space – the function space
temperature – The temperature
t – the time
K_S – The solubility of the species. If provided, the value of the
condition (boundary) – is divided by K_S (change of variable method).

create_value_ufl(temperature: float | Constant | Function)[source]#

Creates the value of the boundary condition as a pure ufl expression and sets it to self.value_fenics.

Unlike create_value(), this never interpolates the value into a dolfinx.fem.Function. It is required when the value depends on an unknown of the problem, such as the pressure of a festim.Enclosure, which lives in a real function space and cannot be interpolated. Such a BC can therefore only be enforced weakly.

Parameters:: temperature – the temperature
Raises:: ValueError – if the value is not a callable of the temperature only

stoichiometry: float = 1#: number of particles of species in the solid per particle of the gas species the BC is coupled to. 1 by default (Henry’s law, where the gas molecule dissolves as such); 2 for Sieverts’ law, where a diatomic molecule dissolves as two atoms. Only used when the BC is coupled to a festim.GasSpecies.

Bases: DirichletBCBase

create_value(function_space: FunctionSpace, t: Constant)[source]#

Creates the value of the boundary condition as a fenics object and sets it to self.value_fenics. If the value is a constant, it is converted to a dolfinx.fem.Constant. If the value is a function of t, it is converted to a dolfinx.fem.Constant. Otherwise, it is converted to a` dolfinx.fem.Function` and the expression of the function is stored in bc_expr.

Parameters:

function_space – the function space
t – the time

class FluxBCBase(subdomain: SurfaceSubdomain, value)[source]#

Bases: object

Flux boundary condition class.

Ensuring the gradient of the solution u at a boundary:

\[-A \nabla u \cdot \mathbf{n} = f\]

where \(A\) is some material property (diffusivity for particle flux and thermal conductivity for heat flux), \(\mathbf{n}\) is the outwards normal vector of the boundary, \(f\) is a function of space and time.

Parameters:

subdomain (festim.SurfaceSubdomain) – the surface subdomain where the boundary condition is applied
value (float, fem.Constant, callable) – the value of the boundary condition

Variables:

subdomain (festim.SurfaceSubdomain) – the surface subdomain where the boundary condition is applied
value (float, fem.Constant, callable) – the value of the boundary condition
value_fenics (fem.Function or fem.Constant) – the value of the boundary
in (condition) – fenics format
bc_expr (fem.Expression) – the expression of the boundary condition that
to (is used) – update the value_fenics

create_value_fenics(mesh, temperature, t: Constant)[source]#

Creates the value of the boundary condition as a fenics object and sets it to self.value_fenics. If the value is a constant, it is converted to a fenics.Constant. If the value is a function of t, it is converted to a fenics.Constant. Otherwise, it is converted to a ufl Expression.

Parameters:

mesh (dolfinx.mesh.Mesh) – the mesh
temperature (float) – the temperature
t (dolfinx.fem.Constant) – the time

update(t)[source]#

Updates the flux bc value.

Parameters:: t (float) – the time

class HeatFluxBC(subdomain, value)[source]#

Bases: FluxBCBase

Heat flux boundary condition class.

Ensuring the gradient of the solution T at a boundary: -lambda * grad(T) * n = f where lambda is the thermal conductivity , n is the outwards normal vector of the boundary, f is a function of space and time.

Parameters:

subdomain (festim.SurfaceSubdomain) – the surface subdomain where the heat flux is applied
value (float, callable, fem.Constant) – the value of the heat flux

Variables:

subdomain (festim.SurfaceSubdomain) – the surface subdomain where the heat flux is applied
value (float, callable, fem.Constant) – the value of the heat flux
value_fenics (fem.Function or fem.Constant) – the value of the heat flux in fenics format
bc_expr (fem.Expression) – the expression of the heat flux that is used to update the value_fenics

Examples

HeatFluxBC(subdomain=my_subdomain, value=1)
HeatFluxBC(subdomain=my_subdomain, value=lambda x: 1 + x[0])
HeatFluxBC(subdomain=my_subdomain, value=lambda t: 1 + t)
HeatFluxBC(subdomain=my_subdomain, value=lambda x, t: 1 + x[0] + t)

class HenrysBC(subdomain, H_0, E_H, pressure, species, enforce_weakly=False, penalty=None)[source]#

Bases: FixedConcentrationBC

Henrys boundary condition class.

c = H * pressure H = H_0 * exp(-E_H / k_B / T)

Parameters:

subdomain (festim.Subdomain) – the subdomain where the boundary condition is applied
species (str) – the name of the species
H_0 (float or fem.Constant) – the Henrys constant pre-exponential factor
(H/m3/Pa)
E_H (float or fem.Constant) – the Henrys constant activation energy (eV)
pressure (float or callable) – the pressure at the boundary (Pa)

Variables:

subdomain (festim.Subdomain) – the subdomain where the boundary condition is applied
value (float or fem.Constant) – the value of the boundary condition
species (festim.Species or str) – the name of the species
H_0 (float or fem.Constant) – the Henrys constant pre-exponential factor
(H/m3/Pa)
E_H (float or fem.Constant) – the Henrys constant activation energy (eV)
pressure (float or callable) – the pressure at the boundary (Pa)

Examples

HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=1e5,
species="H")
HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda
x: 1e5 + x[0], species="H")
HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda
t: 1e5 + t, species="H")
HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda
T: 1e5 + T, species="H")
HenrysBC(subdomain=my_subdomain, H_0=1e-6, E_H=0.2, pressure=lambda
x, t: 1e5 + x[0] + t, species="H")

create_new_value_function()[source]#

Creates a new value function based on the pressure attribute.

Raises:: ValueError – if the pressure function is not supported
Returns:: the value function
Return type:: callable

stoichiometry: float = 1#: one gas molecule dissolves as one particle of the solid species

class OutflowBC(subdomain: SurfaceSubdomain, species: Species)[source]#

Bases: object

Lets a drift term carry the species out through a boundary.

Adds \(+ \int_{\partial\Omega} c\, (\mathbf{v} \cdot \mathbf{n})\, w\), cancelling the boundary term the divergence form leaves behind, so the natural condition on this surface becomes zero diffusive flux – the standard “do-nothing” outflow of advection-diffusion. Without it an untagged outlet is a closed end and the species backs up against it. See Letting the species out in the user guide.

It is a no-op on a surface where no drift term acts on species.

Parameters:

subdomain – the surface subdomain the species flows out through. On a codimensional problem this may be the boundary of a manifold – the outlet of a 1D fluid, for instance
species – the species carried out

Variables:

subdomain (festim.subdomain.surface_subdomain.SurfaceSubdomain) – the surface subdomain the species flows out through
species (festim.species.Species) – the species carried out

Examples

import festim as F

F.OutflowBC(subdomain=outlet, species=H)

class ParticleFluxBC(subdomain, value, species, species_dependent_value=None)[source]#

Bases: FluxBCBase

Particle flux boundary condition class.

Ensuring the gradient of the solution c at a boundary: -D * grad(c) * n = f where D is the material diffusivity, n is the outwards normal vector of the boundary, f is a function of space and time.

Parameters:

subdomain (festim.SurfaceSubdomain) – the surface subdomain where the
flux (particle) – 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 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.

Variables:

subdomain (festim.SurfaceSubdomain) – the surface subdomain where the
flux (particle) – is applied
value (float or fem.Constant) – the value of the particle flux
species (festim.Species) – the species to which the flux is applied
value_fenics (fem.Function or fem.Constant) – the value of the particle flux in 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 mapping the argument names in a callable value to festim.Species objects

Examples

ParticleFluxBC(subdomain=my_subdomain, value=1, species="H")
ParticleFluxBC(subdomain=my_subdomain, value=lambda x: 1 + x[0],
species="H")
ParticleFluxBC(subdomain=my_subdomain, value=lambda t: 1 + t, species="H")
ParticleFluxBC(subdomain=my_subdomain, value=lambda T: 1 + T, species="H")
ParticleFluxBC(subdomain=my_subdomain, value=lambda x, t: 1 + x[0]
+ t, species="H")
ParticleFluxBC(subdomain=my_subdomain, value=lambda c1: 2 * c1**2,
species="H", species_dependent_value={"c1": species1})

create_value_fenics(mesh, temperature, t: Constant)[source]#

Creates the value of the boundary condition as a fenics object and sets it to self.value_fenics. If the value is a constant, it is converted to a fenics.Constant. If the value is a function of t, it is converted to a fenics.Constant. Otherwise, it is converted to a ufl Expression.

Parameters:

mesh (dolfinx.mesh.Mesh) – the mesh
temperature (float) – the temperature
t (dolfinx.fem.Constant) – the time

class SievertsBC(subdomain, S_0, E_S, pressure, species, enforce_weakly=False, penalty=None)[source]#

Bases: FixedConcentrationBC

Sieverts boundary condition class.

c = S * sqrt(pressure) S = S_0 * exp(-E_S / k_B / T)

Parameters:

subdomain (festim.Subdomain) – the subdomain where the boundary condition is applied
species (str) – the name of the species
S_0 (float or fem.Constant) – the Sieverts constant pre-exponential
factor (H/m3/Pa0.5)
E_S (float or fem.Constant) – the Sieverts constant activation energy (eV)
pressure (float or callable) – the pressure at the boundary (Pa)

Variables:

subdomain (festim.Subdomain) – the subdomain where the boundary condition is applied
value (float or fem.Constant) – the value of the boundary condition
species (festim.Species or str) – the name of the species
S_0 (float or fem.Constant) – the Sieverts constant pre-exponential
factor (H/m3/Pa0.5)
E_S (float or fem.Constant) – the Sieverts constant activation energy (eV)
pressure (float or callable) – the pressure at the boundary (Pa)

Examples

SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2, pressure=1e5,
species="H")
SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2,
pressure=lambda x: 1e5 + x[0], species="H")
SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2,
pressure=lambda t: 1e5 + t, species="H")
SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2,
pressure=lambda T: 1e5 + T, species="H")
SievertsBC(subdomain=my_subdomain, S_0=1e-6, E_S=0.2,
pressure=lambda x, t: 1e5 + x[0] + t, species="H")

create_new_value_function()[source]#

Creates a new value function based on the pressure attribute.

Raises:: ValueError – if the pressure function is not supported
Returns:: the value function
Return type:: callable

stoichiometry: float = 2#: a diatomic gas molecule dissolves as two particles of the solid species

class SurfaceReactionBC(reactant, gas_pressure, k_r0, E_kr, k_d0, E_kd, subdomain)[source]#

Bases: object

Boundary condition representing a surface reaction A + B <-> C.

where A, B are the reactants and C is the product the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T)) and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T)) The reaction rate is: K = K_r * C_A * C_B - K_d * P_C with C_A, C_B the concentration of species A and B, P_C the partial pressure of species C at the surface.

The flux of species A entering the surface is K. In the special case where A=B, then the flux of particle entering the surface is 2*K

Parameters:

reactant (list) – list of F.Species objects representing the reactants
gas_pressure (float, callable or F.GasSpecies) – the partial pressure of the product species. If a F.GasSpecies is given, the pressure is an unknown of the problem and the reaction feeds the mass balance of the enclosure that species belongs to.
k_r0 (float) – the pre-exponential factor of the forward reaction rate
E_kr (float) – the activation energy of the forward reaction rate (eV)
k_d0 (float) – the pre-exponential factor of the backward reaction rate
E_kd (float) – the activation energy of the backward reaction rate (eV)
subdomain (F.SurfaceSubdomain) – the surface subdomain where the reaction occurs

class SurfaceReactionBCpartial(reactant, gas_pressure, k_r0, E_kr, k_d0, E_kd, subdomain, species)[source]#

Bases: ParticleFluxBC

Boundary condition representing a surface reaction A + B <-> C.

where A, B are the reactants and C is the product the forward reaction rate is K_r = k_r0 * exp(-E_kr / (k_B * T)) and the backward reaction rate is K_d = k_d0 * exp(-E_kd / (k_B * T)) The reaction rate is: K = K_r * C_A * C_B - K_d * P_C with C_A, C_B the concentration of species A and B, P_C the partial pressure of species C at the surface.

This class is used to create the flux of a single species entering the surface Example: The flux of species A entering the surface is K.

Parameters:

reactant (list) – list of F.Species objects representing the reactants
gas_pressure (float, callable or F.GasSpecies) – the partial pressure of the product species. If a F.GasSpecies is given, the pressure is an unknown of the problem and the reaction feeds the mass balance of the enclosure that species belongs to.
k_r0 (float) – the pre-exponential factor of the forward reaction rate
E_kr (float) – the activation energy of the forward reaction rate (eV)
k_d0 (float) – the pre-exponential factor of the backward reaction rate
E_kd (float) – the activation energy of the backward reaction rate (eV)
subdomain (F.SurfaceSubdomain) – the surface subdomain where the reaction occurs
species (F.Species) – the species to which the flux is applied

create_value_fenics(mesh, temperature, t: Constant)[source]#

Creates the value of the boundary condition as a fenics object and sets it to self.value_fenics. If the value is a constant, it is converted to a fenics.Constant. If the value is a function of t, it is converted to a fenics.Constant. Otherwise, it is converted to a ufl Expression.

Parameters:

mesh (dolfinx.mesh.Mesh) – the mesh
temperature (float) – the temperature
t (dolfinx.fem.Constant) – the time

class Enclosure(volume: float, species: list[GasSpecies], temperature: float | Callable, surfaces: dict[SurfaceSubdomain, float] | list[SurfaceSubdomain] | None = None, openings: list[OpeningBase] | None = None, gas_constant: float = 1.380649e-23, name: str | None = None)[source]#

Bases: object

A gas enclosure in contact with the model through one or more surfaces.

The partial pressure of each gas species in the enclosure is an unknown of the problem, solved together with the transport problem. For each species, the pressure evolves as:

\[\frac{dP}{dt} = \frac{k_B T}{V} \left( \sum_\Gamma A_\Gamma \int_\Gamma \varphi \, dS + \sum_\text{openings} Q \right)\]

where \(\varphi\) is the rate of particles entering the gas from the solid, \(Q\) the flow rate through the openings, and \(A_\Gamma\) the physical area of each contact surface (see surfaces).

Parameters:

volume – the volume of the enclosure (m3)
species – the gas species in the enclosure
temperature – the temperature of the gas (K). Can be a callable of time. This is independent of the temperature of the transport problem.
surfaces –
the surfaces in contact with the enclosure, as a dict mapping each festim.SurfaceSubdomain to its physical area. The area is what turns the flux through the surface into a number of particles per second, and the mesh only provides it in 3D:
- 1D: a surface is a point and carries no extent, so the area is the area of the membrane facing the enclosure, in m2. Required.
- 2D: a surface is a line, so the area is the out-of-plane depth of the model, in m. Required.
- 3D: the mesh already measures the area, so pass 1.0. A plain list of surfaces is accepted in 3D and means an area of 1.0 for each.
An enclosure with no contact surfaces is allowed (it then only exchanges through its openings).
openings – the openings of the enclosure (see festim.Pump, festim.Reservoir, festim.PrescribedFlowRate, festim.EnclosureConnection)
gas_constant – the constant relating pressure to particle density (\(P = n k T\)). Defaults to festim.k_B_SI (J/K), which matches FESTIM’s convention of concentrations in particles/m3. Pass festim.R (J/mol/K) if working in mol/m3.
name – a name given to the enclosure

Variables:

volume – the volume of the enclosure (m3)
species – the gas species in the enclosure
temperature – the temperature of the gas, wrapped in a festim.helpers.Value
surfaces – a dict mapping each contact surface to its physical area
openings – the openings of the enclosure
gas_constant – the constant relating pressure to particle density
name – a name given to the enclosure

Examples

H2 = F.GasSpecies(name="H2", initial_pressure=1e5)
my_enclosure = F.Enclosure(
    volume=1e-3,
    species=[H2],
    temperature=500,
    openings=[F.Pump(pumping_speed=1e-4)],
)

convert_input_values_to_fenics_objects(function_space, t)[source]#

Converts the user input values of the enclosure and its openings to fenics objects.

Parameters:

function_space – a function space on the parent mesh
t – the time, as a fenics Constant

property thermal_energy#

The quantity relating pressure to particle density (\(P = n k T\)).

Returns gas_constant * T as a fenics object once the problem has been initialised.

update_time_dependent_values(t: float)[source]#

Updates the time-dependent values of the enclosure and its openings.

Parameters:: t – the time

class EnclosureConnection(conductance: float | Callable, species)[source]#

Bases: OpeningBase

An opening connecting two enclosures, coupling their pressures.

The flow rate into the enclosure holding species[0] is C * (P_1 - P_0) / (k * T), and the opposite for the other side.

The connection only needs to be declared in the openings of one of the two enclosures: the mirrored term is added to the partner enclosure automatically when the problem is initialised.

Parameters:

conductance – the conductance of the connection (m3/s). Can be a callable of time.
species – a tuple of the two gas species being connected, one in each enclosure

Examples

H2_a = GasSpecies(name="H2", initial_pressure=1e5)
H2_b = GasSpecies(name="H2", initial_pressure=0)
connection = EnclosureConnection(conductance=1e-4, species=(H2_a, H2_b))
enclosure_a = Enclosure(
    volume=1e-3, species=[H2_a], temperature=500, openings=[connection]
)
enclosure_b = Enclosure(volume=2e-3, species=[H2_b], temperature=500)

molar_flow_rate(gas_species, enclosure)[source]#

The flow rate of particles into the enclosure (particles/s).

A positive value means particles entering the enclosure.

Parameters:

gas_species – the gas species this rate is computed for
enclosure – the enclosure the opening belongs to

Returns:

a ufl expression for the flow rate

class GasSpecies(name: str, initial_pressure: float = 0.0)[source]#

Bases: object

A gas species living in a festim.Enclosure.

The partial pressure of the species in the enclosure is an unknown of the problem, represented by a real function space (one global degree of freedom).

Parameters:

name – a name given to the species
initial_pressure – the partial pressure at t=0 (Pa)

Variables:

name – a name given to the species
initial_pressure – the partial pressure at t=0 (Pa)
enclosure – the enclosure this species belongs to. Set by festim.Enclosure
function_space – the real function space holding the pressure
solution – the pressure at the current timestep
prev_solution – the pressure at the previous timestep
test_function – the test function of the real function space
F – the variational formulation of the pressure balance

Examples

H2 = GasSpecies(name="H2", initial_pressure=1e5)
my_enclosure = Enclosure(volume=1e-3, species=[H2], temperature=500)

property pressure#: The pressure as a fenics object, to be used in boundary conditions.

property value: float#: The current pressure as a float (Pa), collected across all MPI processes.

class OpeningBase(species=None)[source]#

Bases: object

Base class for enclosure openings.

An opening lets gas in or out of an enclosure, modifying its mass balance. Subclasses implement molar_flow_rate().

Parameters:: species – the gas species this opening applies to. If None, the opening applies to every species in the enclosure.
Variables:: species – the gas species this opening applies to, or None for all of them

applies_to(gas_species) → bool[source]#: Whether this opening acts on a given gas species.

convert_input_values_to_fenics_objects(function_space, t)[source]#

Converts the user input values to fenics objects.

Parameters:

function_space – a function space on the parent mesh. Only its mesh is used, since opening parameters are scalars.
t – the time, as a fenics Constant

molar_flow_rate(gas_species, enclosure)[source]#

The flow rate of particles into the enclosure (particles/s).

A positive value means particles entering the enclosure.

Parameters:

gas_species – the gas species this rate is computed for
enclosure – the enclosure the opening belongs to

Returns:

a ufl expression for the flow rate

class PrescribedFlowRate(flow_rate: float | Callable, species=None)[source]#

Bases: OpeningBase

An opening with a directly imposed flow rate, independent of the pressure.

Parameters:

flow_rate – the flow rate of particles into the enclosure (particles/s). A negative value removes particles. Can be a callable of time.
species – the gas species this opening applies to. If None, applies to all.

Examples

H2 = GasSpecies(name="H2")
my_enclosure = Enclosure(
    volume=1e-3,
    species=[H2],
    temperature=500,
    openings=[PrescribedFlowRate(flow_rate=1e18, species=H2)],
)

molar_flow_rate(gas_species, enclosure)[source]#

The flow rate of particles into the enclosure (particles/s).

A positive value means particles entering the enclosure.

Parameters:

gas_species – the gas species this rate is computed for
enclosure – the enclosure the opening belongs to

Returns:

a ufl expression for the flow rate

class Pump(pumping_speed: float | Callable, species=None)[source]#

Bases: OpeningBase

An opening to vacuum with a given pumping speed.

The flow rate out of the enclosure is S * P / (k * T), giving an exponential decay of the pressure P(t) = P_0 * exp(-S * t / V) for a closed enclosure.

Parameters:

pumping_speed – the pumping speed (m3/s). Can be a callable of time.
species – the gas species this pump applies to. If None, applies to all of them.

Examples

H2 = GasSpecies(name="H2", initial_pressure=1e5)
my_enclosure = Enclosure(
    volume=1e-3,
    species=[H2],
    temperature=500,
    openings=[Pump(pumping_speed=1e-4)],
)

molar_flow_rate(gas_species, enclosure)[source]#

The flow rate of particles into the enclosure (particles/s).

A positive value means particles entering the enclosure.

Parameters:

gas_species – the gas species this rate is computed for
enclosure – the enclosure the opening belongs to

Returns:

a ufl expression for the flow rate

class Reservoir(conductance: float | Callable, pressure: float | Callable, species=None)[source]#

Bases: OpeningBase

An opening to an external reservoir held at a given pressure.

The flow rate into the enclosure is C * (P_ext - P) / (k * T), giving P(t) = P_ext + (P_0 - P_ext) * exp(-C * t / V) for an enclosure with no other exchange.

Parameters:

conductance – the conductance of the opening (m3/s). Can be a callable of time.
pressure – the pressure of the reservoir (Pa). Can be a callable of time.
species – the gas species this opening applies to. If None, applies to all.

Examples

H2 = GasSpecies(name="H2", initial_pressure=1e5)
my_enclosure = Enclosure(
    volume=1e-3,
    species=[H2],
    temperature=500,
    openings=[Reservoir(conductance=1e-4, pressure=1e3)],
)

molar_flow_rate(gas_species, enclosure)[source]#

The flow rate of particles into the enclosure (particles/s).

A positive value means particles entering the enclosure.

Parameters:

gas_species – the gas species this rate is computed for
enclosure – the enclosure the opening belongs to

Returns:

a ufl expression for the flow rate

class AverageSurface(field: Species | str, surface: SurfaceSubdomain | VolumeSubdomain | int, filename: str | None = None)[source]#

Bases: SurfaceQuantity

Computes the average value of a field on a given surface.

Parameters:

field (festim.Species) – species for which the average surface is computed
surface (festim.SurfaceSubdomain) – surface subdomain
filename (str, optional) – name of the file to which the average surface is exported

Variables:

festim.SurfaceQuantity (see)

compute(u: Function | Indexed, ds: Measure, entity_maps=None, restriction: str | None = None, subdomain_id: int | None = None)[source]#

Computes the average value of the field on the defined surface subdomain, and appends it to the data list.

Parameters:

u – field for which the average value is computed
ds – surface measure of the model
entity_maps – entity maps relating parent mesh and submesh
restriction – which side of an interior facet to read the field on, when ds is an interior facet measure
subdomain_id – the id to index ds with, when it is not the surface’s own. One side of a manifold adjacent to more than two volume subdomains is integrated under an id of its own

class AverageVolume(field: Species | str, volume: VolumeSubdomain | int, filename: str | None = None)[source]#

Bases: VolumeQuantity

Computes the average value of a field in a given volume.

Parameters:

field (festim.Species) – species for which the average volume is computed
volume (festim.VolumeSubdomain) – volume subdomain
filename (str, optional) – name of the file to which the average volume
exported (is)

Variables:

festim.VolumeQuantity (see)

compute(u, dx, entity_maps=None)[source]#: Computes the average value of solution function within the defined volume subdomain, and appends it to the data list.

class CustomFieldExport(filename: str | Path, expression: Callable, species_dependent_value: dict[str, Species] | None = None, times: list[float] | list[int] | None = None, subdomain: VolumeSubdomain = None, format: str = 'vtx', backend: str | None = None, checkpoint: bool = False)[source]#

Bases: FieldExportBase

Export a custom field to a VTX file

Parameters:

filename – The name of the output file
expression – A function evaluating the custom field. Positional arguments of the function can be “t” (time), “x” (spatial coordinate), “T” (temperature), or any key from the species_dependent_value dictionary.
species_dependent_value – A dictionary mapping argument names in expression to Species objects. Defaults to None.
times – if provided, the field will be exported at these timesteps. Otherwise exports at all timesteps. Defaults to None.
subdomain – The volume subdomain on which the custom field is evaluated. Defaults to None.
format – The output format, see FieldExportBase. Defaults to "vtx".
backend – The io4dolfinx backend, for format="checkpoint".
checkpoint – Deprecated, use format="checkpoint".

Variables:

filename – The name of the output file
expression (collections.abc.Callable) – A function evaluating the custom field.
species_dependent_value (dict[str, festim.species.Species]) – A dictionary mapping argument names to Species objects.
subdomain (festim.subdomain.volume_subdomain.VolumeSubdomain) – The volume subdomain on which the custom field is evaluated.
checkpoint – True if the export is a checkpoint file.
times – if provided, the field will be exported at these timesteps. Otherwise exports at all timesteps.
function (dolfinx.fem.function.Function) – the function containing the custom field values
writer (festim.exports.writers.FieldWriter | None) – The writer object used to write the file
dolfinx_expression (dolfinx.fem.function.Expression) – the dolfinx expression used to evaluate the function

check_valid_inputs(kwargs: dict)[source]#

Check if we are in the mixed domain/discontinuous case and if the user-provided expression is valid in this case. dolfinx.fem.Expression does not support co-dim 0 submeshes and time is defined on the parent mesh, so we cannot have time-dependent custom fields in the mixed domain/discontinuous case.

When FEniCS/dolfinx#3207 is resolved we should be able to support this

property mixed_domain: bool#

Check if we are in a mixed domain/discontinuous case. This is the case if at least one of the species in species_dependent_value is defined on a subdomain or if the custom field is defined on a subdomain.

Returns:: True if we are in a mixed domain/discontinuous case, False otherwise.

set_dolfinx_expression(temperature: Constant | Function, time: Constant)[source]#

Set the dolfinx expression used to evaluate the custom field. This is done by evaluating the user-provided expression with the appropriate arguments and using the result to create a dolfinx expression.

Parameters:

temperature – The temperature field to use in the expression
time – The time to use in the expression

update()[source]#: Re-evaluate the custom field before writing.

class CustomQuantity(expr: Callable, subdomain: SurfaceSubdomain | VolumeSubdomain, title: str = 'Custom Quantity', filename: str | None = None)[source]#

Bases: DerivedQuantity

Export CustomQuantity.

Parameters:

expr – function that returns a UFL expression
subdomain – subdomain on which the quantity is evaluated
title – title of the exported quantity
filename – name of the file to which the quantity is exported

Variables:

expr – function that returns a UFL expression
subdomain – subdomain on which the quantity is evaluated
title – title of the exported quantity
filename (str | None) – name of the file to which the quantity is exported
t (list[float]) – list of time values
data (list[float]) – list of values of the quantity

Usage:

import numpy as np
import festim as F

material = F.Material(D_0=1, E_D=0)
volume = F.VolumeSubdomain(id=1, material=material)
surface = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1))

def total_concentration(**kwargs):
    return kwargs["A"] + kwargs["B"]

quantity = F.CustomQuantity(
    expr=total_concentration,
    subdomain=volume,
    title="Total quantity",
)

surface_quantity = F.CustomQuantity(
    expr=lambda **kwargs: -kwargs["D_A"] * ufl.dot(
        ufl.grad(kwargs["A"]), kwargs["n"]
    ),
    subdomain=surface,
    title="Surface flux",
)

The callable passed to expr receives keyword arguments assembled by the problem class. Common entries are:

A, B, …: Concentrations of the species present in the problem (here A and B).
n: The facet normal on the selected surface subdomain.
T: The temperature field.
D_A, D_B, …: Species-specific diffusion coefficients.
D: The diffusion coefficient data, either a single field for one species or a dictionary keyed by species name when several species are present.
x: The spatial coordinate (x[0], x[1], x[2]).

For a surface quantity, the returned UFL expression can represent a flux such as

\[q = -D\,\nabla c \cdot n\]

and FESTIM will assemble

\[Q = \int_{\Gamma} q\,\mathrm{d}\Gamma\]

over the selected surface subdomain.

The expression returned by expr is treated as an integrand and assembled over the chosen subdomain.

\[Q = \int_{\Omega} q\,\mathrm{d}\Omega\]

where q is the UFL expression returned by expr and \Omega is either a surface or a volume subdomain.

compute(measure: Measure, entity_maps: dict | None = None)[source]#

Computes the value of the custom quantity and appends it to the data list.

Parameters:

measure – volume or surface measure of the model
entity_maps – entity maps relating parent mesh and submesh

class DerivedQuantity(filename: str | None = None)[source]#

Bases: ABC

Base class for all derived quantities.

Variables:

filename (str | None) – name of the file to which the quantity is exported
t (list[float]) – list of time values
data (list[float]) – list of values of the quantity

write(t)[source]#: If the filename doesnt exist yet, create it and write the header, then append the time and value to the file.

class FieldExportBase(filename: str | Path, times: list[float] | list[int] | None = None, format: str = 'vtx', backend: str | None = None)[source]#

Bases: object

Base class for exports of fields to a file.

Parameters:

filename – The name of the output file. If its extension doesn’t match the chosen format, the correct one is substituted and a warning is issued.
times – if provided, the field will be exported at these timesteps. Otherwise exports at all timesteps. Defaults to None.
format –
the output format. One of:
- "vtx": .bp, readable by ParaView (the default)
- "vtkhdf": .vtkhdf, readable by ParaView. A single scalable HDF5 file; several exports sharing a filename become blocks of one file.
- "checkpoint": for reloading with festim.read_function_from_file(), not readable by ParaView
- "xdmf": .xdmf + .h5, readable by ParaView
backend – only used by format="checkpoint": "adios2" (default, .bp) or "h5py" (.h5).

Variables:

filename – The name of the output file
times – The timesteps to export at, or None for all timesteps
format (str) – The output format
backend (str | None) – The io4dolfinx backend, for format="checkpoint"
writer (festim.exports.writers.FieldWriter | None) – The festim.exports.writers.FieldWriter doing the writing

close() → None[source]#: Release the underlying file.

define_writer(functions: list[Function], names: list[str], mesh: Mesh, block_name: str = 'mesh', overwrite: bool = True) → None[source]#

Create the underlying writer for this export’s format.

Parameters:

functions – the functions to write at every timestep
names – the name to store each function under, one per function
mesh – the mesh the functions live on
block_name – name of the block, for formats holding several meshes per file
overwrite – False if another export already initialised this file

subdomain: VolumeSubdomain | None = None#: subdomain the export lives on, overridden by subclasses that have one

update() → None[source]#: Refresh the data to write. Called before each write; no-op by default.

write(t: float) → None[source]#: Write the fields at time t.

class GasPressure(field: GasSpecies, filename: str | None = None)[source]#

Bases: DerivedQuantity

Exports the partial pressure of a gas species in an enclosure over time.

Parameters:

field – the gas species to export the pressure of
filename – name of the file to which the pressure is exported

Variables:

field – the gas species to export the pressure of
filename (str | None) – name of the file to which the pressure is exported
t (list[float]) – list of time values
data (list[float]) – list of pressure values (Pa)
value – the pressure at the last computed timestep (Pa)

Examples

H2 = F.GasSpecies(name="H2", initial_pressure=1e5)
my_enclosure = F.Enclosure(volume=1e-3, species=[H2], temperature=500)
my_export = F.GasPressure(field=H2, filename="pressure.csv")

compute()[source]#: Computes the pressure of the gas species and appends it to data.

class MaximumSurface(field: Species | str, surface: SurfaceSubdomain | VolumeSubdomain | int, filename: str | None = None)[source]#

Bases: SurfaceQuantity

Computes the maximum value of a field on a given surface.

Parameters:

field (festim.Species) – species for which the maximum surface is computed
surface (festim.SurfaceSubdomain) – surface subdomain
filename (str, optional) – name of the file to which the maximum surface
exported (is)

Variables:

festim.SurfaceQuantity (see)
facet_meshtags (dolfinx.mesh.MeshTags | None) – the facet meshtags of the parent mesh
volume (festim.subdomain.volume_subdomain.VolumeSubdomain | None) – the volume subdomain the surface bounds. Set by the problem; None outside festim.HydrogenTransportProblemDiscontinuous

compute()[source]#: Computes the maximum value of the field on the defined surface subdomain, and appends it to the data list.

property is_submesh: bool#: Whether the field’s solution lives on the submesh of volume, as in festim.HydrogenTransportProblemDiscontinuous. See issue #1191.

property meshtags#: Facet meshtags of whichever mesh solution lives on.

class MaximumVolume(field: Species | str, volume: VolumeSubdomain | int, filename: str | None = None)[source]#

Bases: VolumeQuantity

Computes the maximum value of a field in a given volume.

Parameters:

field (festim.Species) – species for which the maximum volume is computed
volume (festim.VolumeSubdomain) – volume subdomain
filename (str, optional) – name of the file to which the maximum volume
exported (is)

Variables:

festim.VolumeQuantity (see)
volume_meshtags (dolfinx.mesh.MeshTags | None) – the cell meshtags of the mesh the field is defined on. Is None when the field is defined on a submesh that already coincides with the volume subdomain (festim.HydrogenTransportProblemDiscontinuous)

compute()[source]#: Computes the maximum value of solution function within the defined volume subdomain, and appends it to the data list.

property is_submesh: bool#: Whether the field’s solution lives on a submesh of volume itself, as in festim.HydrogenTransportProblemDiscontinuous. See issue #1191.

class MinimumSurface(field: Species | str, surface: SurfaceSubdomain | VolumeSubdomain | int, filename: str | None = None)[source]#

Bases: SurfaceQuantity

Computes the minimum value of a field on a given surface.

Parameters:

field (festim.Species) – species for which the minimum surface is computed
surface (festim.SurfaceSubdomain) – surface subdomain
filename (str, optional) – name of the file to which the minimum surface
exported (is)

Variables:

festim.SurfaceQuantity (see)
facet_meshtags (dolfinx.mesh.MeshTags | None) – the facet meshtags of the parent mesh
volume (festim.subdomain.volume_subdomain.VolumeSubdomain | None) – the volume subdomain the surface bounds. Set by the problem; None outside festim.HydrogenTransportProblemDiscontinuous

compute()[source]#: Computes the minimum value of the field on the defined surface subdomain, and appends it to the data list.

property is_submesh: bool#: Whether the field’s solution lives on the submesh of volume, as in festim.HydrogenTransportProblemDiscontinuous. See issue #1191.

property meshtags#: Facet meshtags of whichever mesh solution lives on.

class MinimumVolume(field: Species | str, volume: VolumeSubdomain | int, filename: str | None = None)[source]#

Bases: VolumeQuantity

Computes the minimum value of a field in a given volume.

Parameters:

field (festim.Species) – species for which the minimum volume is computed
volume (festim.VolumeSubdomain) – volume subdomain
filename (str, optional) – name of the file to which the minimum volume
exported (is)

Variables:

festim.VolumeQuantity (see)
volume_meshtags (dolfinx.mesh.MeshTags | None) – the cell meshtags of the mesh the field is defined on. Is None when the field is defined on a submesh that already coincides with the volume subdomain (festim.HydrogenTransportProblemDiscontinuous)

compute()[source]#: Computes the minimum value of solution function within the defined volume subdomain, and appends it to the data list.

property is_submesh: bool#: Whether the field’s solution lives on a submesh of volume itself, as in festim.HydrogenTransportProblemDiscontinuous. See issue #1191.

class Profile1DExport(field: Species, subdomain: VolumeSubdomain = None, times: list[float] | None = None)[source]#

Bases: object

Class to export 1D profiles of a field in a simulation.

Parameters:

field – the species for which the profile is computed.
subdomain – the volume subdomain to compute the profile on. If None, the profile is computed over the entire domain.
times – if provided, the profile will be exported at these timesteps. Otherwise, exports at all timesteps. Defaults to None.

Variables:

field (festim.species.Species) – the species for which the profile is computed.
subdomain (festim.subdomain.volume_subdomain.VolumeSubdomain | None) – the volume subdomain to compute the profile on. If None, the profile is computed over the entire domain.
times (list[float] | None) – if provided, the profile will be exported at these timesteps. Otherwise, exports at all timesteps.
x (numpy.ndarray) – the coordinates along which the profile is computed.
data (list) – the computed profile data.
t (list[float]) – the list of time values at which the profile is computed.

class ReactionRateExport(reaction: ArrheniusReaction, filename: str | Path, direction: str = 'both', times: list[float] | None = None, subdomain: VolumeSubdomain | None = None, format: str = 'vtx', backend: str | None = None, checkpoint: bool = False)[source]#

Bases: CustomFieldExport

Export a reaction rate to a VTX file

Parameters:

reaction – The reaction to export the rate of.
filename – The name of the output file.
direction – The direction of the reaction to export. Can be “forward”, “backward” or “both”. Defaults to “both”.
times – if provided, the field will be exported at these timesteps. Otherwise exports at all timesteps. Defaults to None.
subdomain – The volume subdomain on which the reaction rate is evaluated. Defaults to None.
format – The output format, see FieldExportBase. Defaults to "vtx".
backend – The io4dolfinx backend, for format="checkpoint".
checkpoint – Deprecated, use format="checkpoint".

override_signature(expression: Callable, reactant_names: list[str], product_names: list[str])[source]#

Override the signature of the expression function. This is needed to ensure that the expression has the correct arguments for set_dolfinx_expression().

Parameters:: expression – The user-provided expression for the reaction rate. The arguments of the expression must be T (temperature) and the names of the reactants and products.

Bases: FieldExportBase

Export species concentration fields to a file.

Parameters:

filename – The name of the output file
field – Set of species to export
subdomain – A field can be defined on multiple domains. This arguments specifies what subdomains we export on. If None we export on all domains.
format – The output format, see FieldExportBase. Defaults to "vtx".
backend – The io4dolfinx backend, for format="checkpoint".
times – if provided, the field will be exported at these timesteps. Otherwise exports at all timesteps. Defaults to None.

Variables:

filename – The name of the output file
field – Set of species to export
subdomain (festim.subdomain.volume_subdomain.VolumeSubdomain | None) – The subdomain the species are exported on
times – The timesteps to export at, or None for all timesteps
writer (festim.exports.writers.FieldWriter | None) – The writer object used to write the file

Example

# one ParaView-readable .vtkhdf file
F.SpeciesExport("results.vtkhdf", field=[H], subdomain=vol,
                format="vtkhdf")

# a checkpoint, to restart from later
F.SpeciesExport("state.bp", field=[H], subdomain=vol,
                format="checkpoint")

get_functions() → list[Function][source]#

Returns list of species for a given subdomain.

If using legacy mode, return the whole species.

class SurfaceFlux(field: Species, surface: SurfaceSubdomain | VolumeSubdomain, filename: str | None = None)[source]#

Bases: SurfaceQuantity

Computes the flux of a field on a given surface.

The total flux -D grad(c) . n + c v . n is computed, so a drift term acting on the field – advection, Soret, electromigration – contributes.

Parameters:

field – species for which the surface flux is computed
surface – surface subdomain. See festim.SurfaceQuantity for the codimensional cases – a manifold, or the boundary of one.
filename – name of the file to which the surface flux is exported

Variables:

drift_velocity – the summed velocity of the drift terms acting on field in the subdomain this surface bounds, set by the problem during initialise(). None when there is none
rest (see festim.SurfaceQuantity for the)

compute(u: Function | Indexed, ds: Measure, entity_maps=None, restriction: str | None = None, subdomain_id: int | None = None)[source]#

Computes the value of the flux at the surface.

Parameters:

u – field for which the flux is computed
ds – surface measure of the model
entity_maps – entity maps relating parent mesh and submesh
restriction – which side of an interior facet to evaluate the flux on, when ds is an interior facet measure. The whole integrand is restricted, so the normal is the one pointing out of that side and the sign convention matches an exterior surface.
subdomain_id – the id to index ds with, when it is not the surface’s own. One side of a manifold adjacent to more than two volume subdomains is integrated under an id of its own

class SurfaceQuantity(field: Species | str, surface: SurfaceSubdomain | VolumeSubdomain | int, filename: str | None = None)[source]#

Bases: DerivedQuantity

Export SurfaceQuantity.

Parameters:

field – species for which the surface flux is computed
surface – surface subdomain. A codim-1 VolumeSubdomain (a manifold) may be given instead, to compute the quantity on the facets it occupies – for a bulk species, since a manifold’s own species has no flux across it. A codim-2 SurfaceSubdomain computes it on the boundary of a manifold.
filename – name of the file to which the surface flux is exported

Variables:

field (festim.species.Species) – species for which the surface flux is computed
surface (festim.subdomain.surface_subdomain.SurfaceSubdomain | festim.subdomain.volume_subdomain.VolumeSubdomain) – surface subdomain
filename (str | None) – name of the file to which the surface flux is exported
t (list[float]) – list of time values
data (list[float]) – list of values of the surface quantity

class TemperatureExport(filename: str | Path, format: str = 'vtx', backend: str | None = None, times: list[float] | list[int] | None = None)[source]#

Bases: FieldExportBase

Export the temperature field to a file.

Parameters:

filename – The name of the output file
format – The output format, see FieldExportBase. Defaults to "vtx".
backend – The io4dolfinx backend, for format="checkpoint".
times – if provided, the field will be exported at these timesteps. Otherwise exports at all timesteps. Defaults to None.

Variables:

filename – The name of the output file
times – The timesteps to export at, or None for all timesteps
writer (festim.exports.writers.FieldWriter | None) – The writer object used to write the file

class TotalSurface(field: Species | str, surface: SurfaceSubdomain | VolumeSubdomain | int, filename: str | None = None)[source]#

Bases: SurfaceQuantity

Computes the total value of a field on a given surface.

Parameters:

field (festim.Species) – species for which the total volume is computed
surface (festim.SurfaceSubdomain) – surface subdomain
filename (str, optional) – name of the file to which the total volume is exported

Variables:

festim.SurfaceQuantity (see)

compute(u: Function | Indexed, ds: Measure, entity_maps=None, restriction: str | None = None, subdomain_id: int | None = None)[source]#

Computes the total value of the field on the defined surface subdomain, and appends it to the data list.

Parameters:

u – field for which the total value is computed
ds – surface measure of the model
entity_maps – entity maps relating parent mesh and submesh
restriction – which side of an interior facet to read the field on, when ds is an interior facet measure
subdomain_id – the id to index ds with, when it is not the surface’s own. One side of a manifold adjacent to more than two volume subdomains is integrated under an id of its own

class TotalVolume(field: Species | str, volume: VolumeSubdomain | int, filename: str | None = None)[source]#

Bases: VolumeQuantity

Computes the total value of a field in a given volume.

Parameters:

field (festim.Species) – species for which the total volume is computed
volume (festim.VolumeSubdomain) – volume subdomain
filename (str, optional) – name of the file to which the total volume is exported

Variables:

festim.VolumeQuantity (see)

compute(u, dx: Measure, entity_maps=None)[source]#

Computes the value of the total volume of the field in the volume subdomain and appends it to the data list.

Parameters:

u – field for which the total volume is computed
dx – volume measure of the model
entity_maps – entity maps relating parent mesh and submesh

class VTXSpeciesExport(filename: str | Path, field: Species | list[Species], subdomain: VolumeSubdomain = None, checkpoint: bool = False, times: list[float] | list[int] | None = None)[source]#

Bases: SpeciesExport

Export species fields to a VTX (.bp) file.

Deprecated since version Use: SpeciesExport, which supports other formats, instead. This class will be removed in a future release.

Parameters:

filename – The name of the output file
field – Set of species to export
subdomain – The subdomain to export on
checkpoint – If True, write a checkpoint instead of a VTX file. Equivalent to SpeciesExport(..., format="checkpoint").
times – The timesteps to export at, or None for all timesteps

class VTXTemperatureExport(filename: str | Path, times: list[float] | list[int] | None = None)[source]#

Bases: TemperatureExport

Export the temperature field to a VTX (.bp) file.

Deprecated since version Use: TemperatureExport, which supports other formats, instead. This class will be removed in a future release.

class VolumeQuantity(field: Species | str, volume: VolumeSubdomain | int, filename: str | None = None)[source]#

Bases: DerivedQuantity

Export VolumeQuantity.

Parameters:

field – species for which the volume quantity is computed
volume – volume subdomain
filename – name of the file to which the volume quantity is exported

Variables:

field (festim.species.Species) – species for which the volume quantity is computed
volume (festim.subdomain.volume_subdomain.VolumeSubdomain) – volume subdomain
filename (str | None) – name of the file to which the volume quantity is exported
t (list[float]) – list of time values
data (list[float]) – list of values of the volume quantity

class XDMFExport(filename: str | Path, field: list[Species] | Species, subdomain: VolumeSubdomain = None, times: list[float] | list[int] | None = None)[source]#

Bases: SpeciesExport

Export species fields to an XDMF file.

Thin wrapper over festim.SpeciesExport with format="xdmf".

Deprecated since version Use: festim.SpeciesExport with format="xdmf" instead. This class will be removed in a future release.

Parameters:

filename – The name of the output file
field – The field(s) to export
subdomain – The subdomain to export on. If None we export on all domains.
times – if provided, the field will be exported at these timesteps. Otherwise exports at all timesteps. Defaults to None.

class CoordinateSystem(value)[source]#

Bases: Enum

classmethod from_string(s: str)[source]#: Can be removed with Python 3.11+.

class Mesh(mesh: Mesh | None = None, coordinate_system: str | CoordinateSystem = CoordinateSystem.CARTESIAN)[source]#

Bases: object

Mesh class.

Parameters:

mesh – The mesh. Defaults to None.
coordinate_system – the coordinate system of the mesh (“cartesian”, “cylindrical”, “spherical”). Defaults to “cartesian”.

Variables:

mesh (mesh The)
vdim (int) – the dimension of the mesh cells
fdim (int) – the dimension of the mesh facets
n (ufl.geometry.FacetNormal) – Symbolic representation of the vector normal to the facets of the mesh.

check_mesh_dim_coords()[source]#: Checks if the used coordinates can be applied for geometry with the specified dimensions.

define_meshtags(surface_subdomains, volume_subdomains, interfaces=None)[source]#

Defines the facet and volume meshtags of the mesh.

Parameters:

surface_subdomains (list of festim.SufaceSubdomains) – the surface
model (subdomains of the)
volume_subdomains (list of festim.VolumeSubdomains) – the volume
model
interfaces (dict, optional) – the interfaces between volume subdomains {int: [VolumeSubdomain, VolumeSubdomain]}. Defaults to None.

Returns:

the facet meshtags dolfinx.mesh.MeshTags: the volume meshtags

Return type:

dolfinx.mesh.MeshTags

class Mesh1D(vertices, **kwargs)[source]#

Bases: Mesh

1D Mesh.

The vertices can be given as a flat list of x-coordinates, or as a list of lists of x-coordinates. In the latter case each sublist produces a disconnected block of cells, so that no cell is created between the last vertex of a block and the first vertex of the next one. This makes it possible to represent several solids separated by a gap (which can then be coupled through e.g. an enclosure).

Parameters:

vertices (list or np.ndarray) – the mesh x-coordinates (m), either flat or as a list of lists (one per disconnected block)

Variables:

vertices (np.ndarray) – all the mesh x-coordinates (m), sorted
vertex_blocks (list of np.ndarray) – the x-coordinates of each disconnected block of the mesh

Examples

import festim as F

# a single continuous domain
F.Mesh1D(vertices=[0, 0.1, 0.2, 0.3])

# two blocks separated by a gap between x=0.3 and x=1
F.Mesh1D(vertices=[[0, 0.1, 0.2, 0.3], [1, 1.1, 1.2]])

check_borders(volume_subdomains)[source]#

Checks that the borders of the subdomain are within the domain.

Parameters:: volume_subdomains (list of festim.VolumeSubdomain1D) – the volume subdomains
Raises:: Value error – if borders outside the domain

define_meshtags(surface_subdomains, volume_subdomains, interfaces=None)[source]#

Defines the facet and volume meshtags of the mesh.

Parameters:

surface_subdomains (list of festim.SufaceSubdomains) – the surface
model (subdomains of the)
volume_subdomains (list of festim.VolumeSubdomains) – the volume
model
interfaces (dict, optional) – the interfaces between volume subdomains {int: [VolumeSubdomain, VolumeSubdomain]}. Defaults to None.

Returns:

the facet meshtags dolfinx.mesh.MeshTags: the volume meshtags

Return type:

dolfinx.mesh.MeshTags

generate_mesh()[source]#: Generates a 1D mesh.

class MeshFromXDMF(volume_file, facet_file, mesh_name='Grid', surface_meshtags_name='Grid', volume_meshtags_name='Grid')[source]#

Bases: Mesh

Mesh read from the XDMF files.

Parameters:

volume_file (str) – path to the volume file
facet_file (str) – path to the facet file
mesh_name (str, optional) – name of the mesh in the volume XDMF file. Defaults to “Grid”.
surface_meshtags_name (str, optional) – name of the surface meshtags in the facet XDMF file. Defaults to “Grid”.
volume_meshtags_name (str, optional) – name of the volume meshtags in the volume XDMF file. Defaults to “Grid”.

Variables:

volume_file (str) – path to the volume file
facet_file (str) – path to the facet file
mesh_name (str) – name of the mesh in the volume XDMF file.
surface_meshtags_name (str) – name of the surface meshtags in the facet
file. (XDMF)
volume_meshtags_name (str) – name of the volume meshtags in the volume XDMF file
mesh (fenics.mesh.Mesh) – the fenics mesh

define_surface_meshtags()[source]#

Creates the facet meshtags.

Returns:: the facet meshtags
Return type:: dolfinx.MeshTags

define_volume_meshtags()[source]#

Creates the volume meshtags.

Returns:: the volume meshtags
Return type:: dolfinx.MeshTags

class Interface(id: int, subdomains: list[VolumeSubdomain], penalty_term: float = 10.0, method: InterfaceMethod = InterfaceMethod.penalty)[source]#

Bases: InterfaceBase

Represents an interface between two subdomains with discontinuous solutions.

This class handles the coupling of solutions across an interface between two volume subdomains using either penalty or Nitsche methods. It manages the exchange of boundary conditions and enforces continuity across the interface.

Variables:

id (int) – Tag of the interface subdomain in the parent mesh tags.
subdomains (tuple[festim.subdomain.volume_subdomain.VolumeSubdomain, festim.subdomain.volume_subdomain.VolumeSubdomain]) – The two subdomains sharing this interface.
parent_mesh (dolfinx.mesh.Mesh) – The parent mesh containing the interface.
mt (dolfinx.mesh.MeshTags) – Mesh tags for the parent mesh.
restriction (list[str, str]) – FEniCS restriction operators for each side of the interface, defaults to (“+”, “-“).
method (festim.subdomain.interface.InterfaceMethod) – The method used to enforce interface conditions (penalty or Nitsche).
penalty_term – Penalty parameter for the interface formulation.

Ds(species: Species, temperature)[source]#

Get diffusion coefficients for both sides of the interface.

Parameters:

species – The species for which to compute diffusivity.
temperature – A function that returns temperature at given restrictions.

Returns:

Diffusion coefficients (D_0, D_1) for subdomains 0 and 1.

Ks(species: Species, temperature)[source]#

Get solubility coefficients for both sides of the interface.

Computes the solubility coefficient at the interface temperature for each subdomain’s material.

Parameters:

species – The species for which to compute solubility.
temperature – A function that returns temperature at given restrictions.

Returns:

Solubility coefficients (K_0, K_1) for subdomains 0 and 1.

equality(species: Species, temperature)[source]#

The interface constraint, as a residual that vanishes at equilibrium.

Both sides are expressed in the same quantity so that the difference is meaningful: the partial pressure when the two materials obey different solubility laws (c/K for Henry, (c/K)**2 for Sievert), and plainly c/K when they obey the same one – for a matching pair the squared and unsquared constraints have the same non-negative roots, so the linear form is preferred as it keeps the coupling linear.

Note that penalty_term therefore carries different units in the two cases, and its values are not comparable across law pairs.

Parameters:

species – The species for which to compute the constraint.
temperature – A function that returns temperature at given restrictions.

Returns:

The constraint residual, zero when the two sides are in equilibrium.

Raises:

ValueError – If either material has an unsupported solubility law.

equality_scale(species: Species, temperature)[source]#

The factor that converts equality() into concentration units.

equality is written in potential units – c/K for a matching pair, a partial pressure for a mixed one – so it cannot be compared to a flux directly. Nitsche’s stabilisation and adjoint terms need it in concentration units, so that penalty_term * D / h * scale * equality is a flux and penalty_term is the dimensionless O(10) stabilisation parameter Nitsche’s theory calls for, whatever units the problem is posed in.

For a matching pair the scale is the mean solubility, which recovers the textbook jump c_0 - c_1 when the two materials share a solubility. For a Sievert/Henry pair it is the Henry coefficient, exactly dc/dP on that side, so scale * equality reads as the concentration the Henry side is missing relative to equilibrium – polynomial in both unknowns, with none of the degeneracy of the Sievert side’s dc/dP = K**2/(2c) at c = 0.

Parameters:

species – The species for which to compute the scale.
temperature – A function that returns temperature at given restrictions.

Returns:

A factor with units of concentration over equality’s units.

get_formulation(dS: Measure, species: list[Species], temperature) → tuple[Form, Form][source]#

Generate the interface formulation for all species.

Parameters:

dS – Integration measure for the interface, with correct integration data.
species – Species for which interface conditions should be applied. Must be defined in both subdomains of the interface.
temperature – Temperature field/function for temperature-dependent laws.

Returns:

Variational forms to be added to each subdomain.

Raises:

AssertionError – If the interface method is unknown or species is not defined in both subdomains.

property method: InterfaceMethod#

Get the interface coupling method.

Returns:: The current interface method (penalty or Nitsche).
Return type:: InterfaceMethod

nitsche_method(dS, species, temperature)[source]#

Generate interface formulation using the Nitsche method.

Nitsche’s method adds, on top of the penalty stabilisation, the term that makes the formulation consistent: the flux the two sides must agree on, {D grad(c) . n}, appears explicitly, so the exact solution satisfies the discrete form for any penalty_term. The penalty then only has to make the system stable rather than to enforce the interface condition on its own, which is why a value of order 10-100 reaches an accuracy the pure penalty needs orders of magnitude more for.

The symmetric (adjoint-consistent) variant is used. Like the penalty, both the consistency and the stabilisation term enter the two sides equally and oppositely, so particles are conserved exactly whatever equality is; and like the penalty it goes through equality(), so a Sievert/Henry pair is coupled through partial pressures rather than through c/K. Unlike the penalty, penalty_term is dimensionless here: the constraint is brought into concentration units by equality_scale() first, so the same value of order 10 works whatever units the problem is posed in.

Parameters:

dS – Integration measure for the interface.
species – The species for which to compute the interface form.
temperature – A function returning temperature at given restrictions.

Returns:

Variational forms for subdomains 0 and 1.

penalty_method(dS, species, temperature)[source]#

Generate interface formulation using the penalty method.

The interface is modelled as a surface exchange: the same flux penalty_term * equality leaves one side and enters the other, so the chemical potential drops across the interface by flux / penalty_term and particles are conserved exactly at any penalty_term.

Parameters:

dS – Integration measure for the interface.
species – The species for which to compute the interface form.
temperature – A function returning temperature at given restrictions.

Returns:

Variational forms for subdomains 0 and 1.

class SurfaceSubdomain(id: int, locator: Callable | None = None, dim: int | None = None)[source]#

Bases: object

Surface subdomain class.

A surface subdomain is a portion of the boundary of a volume subdomain, where boundary conditions and surface exports live.

Parameters:

id – the id of the surface subdomain
locator – a callable function that locates the boundary facets of the subdomain
dim – the topological dimension of the surface. Defaults to None, meaning the facet dimension of the mesh – the boundary of an ordinary (codim-0) volume subdomain. Set it to mesh_dim - 2 to bound a manifold volume subdomain (VolumeSubdomain(dim=mesh_dim - 1)): the endpoints of a line in a 2D mesh, the rim of a surface in a 3D mesh. Such a surface carries no meshtag – its entities are located directly on the manifold’s submesh, and which manifold that is follows from the species of the boundary condition using it.

Examples

SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[0], 1.0))
SurfaceSubdomain(id=1, locator=lambda x:
np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0)))
SurfaceSubdomain(id=1, locator=lambda x:
np.logical_and(np.isclose(x[0], 0.0), np.isclose(x[1], 1.0)))
SurfaceSubdomain(id=1, locator=lambda x:
np.logical_and(np.isclose(x[0], 0.0), x[1] <= 0.5))

codim(mesh_dim: int) → int[source]#

The codimension of the surface in a mesh of dimension mesh_dim.

Parameters:: mesh_dim – the topological dimension of the parent mesh
Returns:: 1 for the boundary of an ordinary volume subdomain, 2 for the boundary of a manifold volume subdomain
Raises:: ValueError – if the resulting codimension is not 1 or 2

locate_boundary_facet_indices(mesh: Mesh) → ndarray[source]#

Locate the boundary entities of the subdomain in mesh.

mesh is the parent mesh for an ordinary surface, and the submesh of the manifold it bounds for a codim-2 one – in both cases the entities searched are the facets of the mesh they are located in.

Parameters:: mesh – a dolfinx mesh object
Raises:: ValueError – if no locator function is provided
Returns:: the list of entities that belong to the subdomain

class SurfaceSubdomain1D(id: int, x: float)[source]#

Bases: SurfaceSubdomain

Surface subdomain class for 1D cases.

Parameters:

id – the id of the surface subdomain
x – the x coordinate of the surface subdomain

Variables:

id (int) – the id of the surface subdomain
x (float) – the x coordinate of the surface subdomain

Examples

SurfaceSubdomain1D(id=1, x=1)

class VolumeSubdomain(id, material, locator: Callable | None = None, name: str | None = None, dim=None)[source]#

Bases: object

Volume subdomain class.

Parameters:

id – the id of the volume subdomain (> 0)
submesh – the submesh of the volume subdomain
cell_map – the cell map of the volume subdomain
parent_mesh – the parent mesh of the volume subdomain
v_map – the vertex map of the volume subdomain
n_map – the normal map of the volume subdomain
ft – the facet meshtags of the volume subdomain
u – the solution function of the subdomain
u_n – the previous solution function of the subdomain
material – the material assigned to the subdomain
sub_T – the sub temperature field in the subdomain
sub_t – for a manifold (codim-1) subdomain, the current time as a constant living on its submesh. None for a codim-0 subdomain, which uses the parent-mesh constant
sub_dt – for a manifold (codim-1) subdomain, the timestep as a constant living on its submesh. None for a codim-0 subdomain
dim – the topological dimension of the subdomain. Defaults to None, meaning the dimension of the mesh. Set it to mesh_dim - 1 to solve a transport equation on a manifold embedded in the mesh (a line in a 2D mesh, a surface in a 3D mesh). Such a subdomain is tagged in the facet meshtags, and can be used wherever a surface is expected (eg. ParticleFluxBC).

codim(mesh_dim: int) → int[source]#

The codimension of the subdomain in a mesh of dimension mesh_dim.

Parameters:: mesh_dim – the topological dimension of the parent mesh
Returns:: 0 for a regular volume subdomain, 1 for a manifold subdomain
Raises:: ValueError – if the resulting codimension is not 0 or 1

create_subdomain(mesh: Mesh, marker: MeshTags)[source]#

Creates the following attributes: .parent_mesh, .submesh, .cell_map, and .v_map.

Only used in festim.HydrogenTransportProblemDiscontinuous

Parameters:

mesh (dolfinx.mesh.Mesh) – the parent mesh
marker (dolfinx.mesh.MeshTags) – the markers the subdomain is tagged in: the volume markers for a codim-0 subdomain, the facet markers for a codim-1 one

locate_subdomain_entities(mesh: Mesh) → ndarray[Any, dtype[int32]][source]#

Locates all entities of the subdomain within the domain.

These are the cells of the mesh for a regular volume subdomain, and the facets for a codim-1 (manifold) one.

Parameters:: mesh – the mesh of the model
Returns:: the entities of the subdomain
Return type:: entities

class VolumeSubdomain1D(id, borders, material)[source]#

Bases: VolumeSubdomain

Volume subdomain class for 1D cases.

Parameters:

id (int) – the id of the volume subdomain
borders (list of float) – the borders of the volume subdomain
material (festim.Material) – the material of the volume subdomain

Variables:

id (int) – the id of the volume subdomain
borders (list of float) – the borders of the volume subdomain
material (festim.Material) – the material of the volume subdomain

Examples

VolumeSubdomain1D(id=1, borders=[0, 1], material=my_mat)

compute_one_sided_interior_facet_data(cell_tags: MeshTags, facets, subdomain: VolumeSubdomain)[source]#

Integration data for the facets of a manifold that touch subdomain, ordered so that "+" is always subdomain.

compute_ordered_interior_facet_data() orders the two sides of facets that all separate the same pair of subdomains. That is not enough for a manifold adjacent to more than two volumes – a grain-boundary network in a polycrystal where every grain is its own subdomain – whose facets separate a different pair from one grain to the next. There, one integral per adjacent grain is used instead of one per manifold: this selects the facets that grain lies on and puts it on "+", so every coupling term is written against a single restriction.

A facet between grains i and j is returned by both calls, once with i on "+" and once with j, which is what lets each side carry its own exchange law. A facet with the same subdomain on both sides is returned unswapped.

Parameters:

cell_tags – the cell meshtags of the parent mesh
facets – the facets of the manifold, as returned by MeshTags.find
subdomain – the volume subdomain to place on the "+" restriction

Returns:

a flat array of (cell_plus, local_facet_plus, cell_minus, local_facet_minus) quadruples, the form accepted by ufl.Measure("dS", subdomain_data=...)

compute_ordered_interior_facet_data(cell_tags: MeshTags, facet_tags: MeshTags, tag: int, subdomain_plus: VolumeSubdomain, subdomain_minus: VolumeSubdomain)[source]#

Integration data for an interior-facet (dS) integral, with the restrictions ordered so that "+" is always subdomain_plus.

DOLFINx’s own ordering of the two cells of an interior facet is arbitrary, so without this any expression that treats the two sides differently – a solubility jump, or a codimensional coupling with different exchange rates on either side – would silently get its sides swapped.

Parameters:

cell_tags – the cell meshtags of the parent mesh, marking every cell adjacent to the tagged facets with the id of the volume subdomain it belongs to
facet_tags – the facet meshtags of the parent mesh
tag – the value identifying the facets to integrate over
subdomain_plus – the volume subdomain to place on the "+" restriction
subdomain_minus – the volume subdomain to place on the "-" restriction

Returns:

(tag, integration_data), the pair accepted by ufl.Measure("dS", subdomain_data=...). integration_data is a flat array of (cell_plus, local_facet_plus, cell_minus, local_facet_minus) quadruples.

Raises:

ValueError – if a tagged facet does not separate the two subdomains

map_manifold_to_volume_subdomains(ft: MeshTags, ct: MeshTags, facet_to_cell: AdjacencyList_int32, volume_subdomains: list[VolumeSubdomain], manifold_subdomains: list[VolumeSubdomain], comm=None) → dict[VolumeSubdomain, list[VolumeSubdomain]][source]#

Maps each codim-1 (manifold) volume subdomain to the volume subdomains it is adjacent to: one for a manifold on the boundary of the domain or buried inside a single subdomain, two for one sitting on an interface, and as many as there are grains for a boundary network threading a polycrystal in which every grain is its own subdomain.

Parameters:

ft – the facet meshtags of the parent mesh
ct – the cell meshtags of the parent mesh
facet_to_cell – the facet to cell connectivity of the parent mesh
volume_subdomains – the list of volume subdomains
manifold_subdomains – the codim-1 volume subdomains
comm – MPI communicator (required for parallel runs)

Returns:

a dictionary mapping each manifold subdomain to its adjacent volume subdomains, sorted by id

Raises:

ValueError – if a manifold is adjacent to no volume subdomain at all

map_surface_to_volume_subdomains(ft: MeshTags, ct: MeshTags, facet_to_cell: AdjacencyList_int32, volume_subdomains: list[VolumeSubdomain], surface_subdomains: list[SurfaceSubdomain], comm=None) → dict[SurfaceSubdomain, VolumeSubdomain][source]#

Maps surface subdomains to volume subdomains based on the facet and cell meshtags and the facet to cell connectivity.

Raises:

AssertionError – if a surface subdomain is connected to multiple volume subdomains

Parameters:

ft – the facet meshtags of the parent mesh
ct – the cell meshtags of the parent mesh
facet_to_cell – the facet to cell connectivity of the parent mesh
volume_subdomains – the list of volume subdomains
surface_subdomains – the list of surface subdomains
comm – MPI communicator (required for parallel runs)

Returns:

a dictionary mapping surface subdomains: to volume subdomains

Return type:

dict[SurfaceSubdomain, VolumeSubdomain]

class AdvectionTerm(velocity: Function, subdomain: VolumeSubdomain, species: Species)[source]#

Bases: DriftTermBase

Advection term class.

Transport by a velocity field given directly, as opposed to one built from a driving gradient (see festim.SoretTerm, festim.ElectromigrationTerm).

Assembled in divergence form, -div(c v), like every drift term. This differs from the v . grad(c) FESTIM assembled before: the two agree in the interior wherever div(v) == 0, which holds for an incompressible flow, but only the divergence form carries species out through a boundary that has no flux condition on it – an outlet, typically. See festim.drift.drift_form().

Parameters:

velocity – the velocity field or function
subdomain – the volume subdomain where the velocity is to be applied
species – the species to which the velocity field is acting on

Variables:

velocity (dolfinx.fem.function.Function) – the velocity field or function
subdomain (festim.subdomain.volume_subdomain.VolumeSubdomain) – the volume subdomain where the velocity is to be applied
species (festim.species.Species) – the species to which the velocity field is acting on

convert_inputs(function_space, t=None, temperature=None)[source]#

Convert the user-given coefficients of this term to fenics objects.

Called by the problem after the temperature is defined, so that a coefficient given as a function of T has one to read.

Parameters:

function_space – the function space of the problem
t – the time, optional
temperature – the temperature, optional

drift_velocity(D, temperature)[source]#

The drift velocity of this term, as a ufl expression.

Parameters:

D – the diffusion coefficient of the species in this subdomain
temperature – the temperature on the mesh the term is assembled on

Returns:

the velocity, a vector-valued ufl expression. ufl.zero when the driving gradient vanishes, in which case the term is skipped.

time_dependent_inputs()[source]#: The user-given coefficients of this term that may depend on time.

class VelocityField(input_value, species_dependent_value: dict[str, Species] | None = None)[source]#

Bases: Value

A class to handle input values of velocity fields from users and convert them to a relevent fenics object.

Parameters:

input_value – The value of the user input

Variables:

input_value (dolfinx.fem.function.Function | collections.abc.Callable) – The value of the user input
fenics_interpolation_expression – The expression of the user input that is used to update the fenics_object
fenics_object (dolfinx.fem.function.Function) – The value of the user input in fenics format
explicit_time_dependent (bool) – True if the user input value is explicitly time dependent
temperature_dependent (bool) – True if the user input value is temperature dependent
vector_function_space (dolfinx.fem.function.FunctionSpace) – the vector function space of the fenics object

convert_input_value(function_space: FunctionSpace, t: Constant | None = None)[source]#

Converts a user given value to a relevent fenics object.

Parameters:

function_space – the function space of the fenics object
t – the time, optional

update(t: Constant)[source]#

Updates the velocity field.

Parameters:: t – the time

Drift terms: transport driven by a gradient other than the species’ own.

A drift term adds a velocity to the species flux:

J = -D grad(c) + c v

Different physics differ only in what sets v: a fluid carries the species along (festim.AdvectionTerm), a temperature gradient drives thermodiffusion (SoretTerm), an electric potential drives charged defects (ElectromigrationTerm). They share the assembly in drift_form().

class DriftTermBase(subdomain: VolumeSubdomain, species: Species | list[Species])[source]#

Bases: ABC

Base class for drift terms.

Subclasses supply drift_velocity(); everything else – validation, the conversion hook and the weak form – is shared.

Parameters:

subdomain – the volume subdomain where the drift applies
species – the species the drift acts on

Variables:

subdomain (festim.subdomain.volume_subdomain.VolumeSubdomain) – the volume subdomain where the drift applies
species (list[festim.species.Species]) – the species the drift acts on

convert_inputs(function_space: FunctionSpace, t: Constant | None = None, temperature=None)[source]#

Convert the user-given coefficients of this term to fenics objects.

Called by the problem after the temperature is defined, so that a coefficient given as a function of T has one to read.

Parameters:

function_space – the function space of the problem
t – the time, optional
temperature – the temperature, optional

abstractmethod drift_velocity(D, temperature)[source]#

The drift velocity of this term, as a ufl expression.

Parameters:

D – the diffusion coefficient of the species in this subdomain
temperature – the temperature on the mesh the term is assembled on

Returns:

the velocity, a vector-valued ufl expression. ufl.zero when the driving gradient vanishes, in which case the term is skipped.

inert_reason() → str[source]#

Why this term’s velocity came out identically zero, for the user.

Reported when the term is dropped, so the message can name the input at fault rather than leaving the user to work out why the physics did nothing.

time_dependent_inputs() → list[Value][source]#: The user-given coefficients of this term that may depend on time.

update_time_dependent_inputs(t: float)[source]#

Re-evaluate the explicitly time-dependent coefficients of this term.

Called once per time step. The default updates every Value returned by time_dependent_inputs().

Parameters:: t – the time

class ElectromigrationTerm(species: Species | list[Species], charge: float, potential, subdomain: VolumeSubdomain)[source]#

Bases: DriftTermBase

Electromigration of a charged species in an electric potential.

The Nernst-Planck drift term of the flux

\[J = -D \nabla c - \frac{z D c}{k_B T} \nabla \varphi\]

with \(\varphi\) in volts and \(z\) the charge number, so that \(k_B\) in eV/K carries the elementary charge. Positive species drift down the potential gradient. The term vanishes when the potential is uniform.

Parameters:

species – the species the drift acts on
charge – the charge number \(z\) of the species (+1 for a proton, +2 for an oxygen vacancy, -1 for an electron)
potential – the electric potential in V. A float, a callable of x, t and/or T, or a fenics object
subdomain – the volume subdomain where the drift applies

Variables:

charge – the charge number of the species
potential – the potential, wrapped in a festim.helpers.Value

Examples

import festim as F

F.ElectromigrationTerm(
    species=hydroxyl,
    charge=1,
    potential=lambda x: 0.5 * (1 - x[0] / 5e-4),
    subdomain=membrane,
)

convert_inputs(function_space, t=None, temperature=None)[source]#

Convert the user-given coefficients of this term to fenics objects.

Called by the problem after the temperature is defined, so that a coefficient given as a function of T has one to read.

Parameters:

function_space – the function space of the problem
t – the time, optional
temperature – the temperature, optional

drift_velocity(D, temperature)[source]#

The drift velocity of this term, as a ufl expression.

Parameters:

D – the diffusion coefficient of the species in this subdomain
temperature – the temperature on the mesh the term is assembled on

Returns:

the velocity, a vector-valued ufl expression. ufl.zero when the driving gradient vanishes, in which case the term is skipped.

inert_reason() → str[source]#

Why this term’s velocity came out identically zero, for the user.

Reported when the term is dropped, so the message can name the input at fault rather than leaving the user to work out why the physics did nothing.

time_dependent_inputs()[source]#: The user-given coefficients of this term that may depend on time.

class SoretTerm(species: Species | list[Species], Q_star, subdomain: VolumeSubdomain)[source]#

Bases: DriftTermBase

Thermodiffusion – the Soret effect, also called thermophoresis.

Adds the temperature-gradient term of the flux

\[J = -D \nabla c - D \frac{Q^* c}{k_B T^2} \nabla T\]

so hydrogen drifts down the temperature gradient for a positive heat of transport. The term vanishes when the temperature is uniform.

Parameters:

species – the species the drift acts on
Q_star – the heat of transport \(Q^*\) in eV. A float, a callable of x, t and/or T, or a fenics object
subdomain – the volume subdomain where the drift applies

Variables:

Q_star – the heat of transport, wrapped in a festim.helpers.Value

Examples

import festim as F

F.SoretTerm(species=mobile_H, Q_star=0.2, subdomain=my_volume)

convert_inputs(function_space, t=None, temperature=None)[source]#

Convert the user-given coefficients of this term to fenics objects.

Called by the problem after the temperature is defined, so that a coefficient given as a function of T has one to read.

Parameters:

function_space – the function space of the problem
t – the time, optional
temperature – the temperature, optional

drift_velocity(D, temperature)[source]#

The drift velocity of this term, as a ufl expression.

Parameters:

D – the diffusion coefficient of the species in this subdomain
temperature – the temperature on the mesh the term is assembled on

Returns:

the velocity, a vector-valued ufl expression. ufl.zero when the driving gradient vanishes, in which case the term is skipped.

inert_reason() → str[source]#

Why this term’s velocity came out identically zero, for the user.

Reported when the term is dropped, so the message can name the input at fault rather than leaving the user to work out why the physics did nothing.

time_dependent_inputs()[source]#: The user-given coefficients of this term that may depend on time.

drift_form(concentration, test_function, velocity, dx: Measure, coordinate_system: CoordinateSystem, mesh)[source]#

The weak form contribution of a drift term, in divergence form.

\[-\int c\, \mathbf{v} \cdot \nabla w\]

which is \(\nabla \cdot (c \mathbf{v})\) integrated by parts. Because the boundary term it leaves behind is the natural boundary condition, the flux boundary conditions of the problem constrain the total flux, diffusive and drift alike, and a boundary with no condition on it is genuinely no-flux.

The metric factor matters here. FESTIM multiplies the equation by m (r in cylindrical, r^2 in spherical) and uses w / m as the test function; those do not cancel in this form, unlike in the non-conservative \(\int (\mathbf{v} \cdot \nabla c)\, w\) that FESTIM assembled before.

Parameters:

concentration – the concentration of the species
test_function – the test function of the species
velocity – the drift velocity, a vector-valued ufl expression
dx – the volume measure, already indexed by subdomain id
coordinate_system – the coordinate system of the mesh
mesh – the mesh the term is assembled on

Returns:

the ufl form to add to the problem’s formulation

Raises:

NotImplementedError – for an unknown coordinate system

is_zero_velocity(velocity) → bool[source]#

Whether a drift velocity is identically zero.

ufl.grad of a spatially constant field is a Zero at construction, and that propagates through the arithmetic of a drift velocity. This is a structural zero, so it is zero at every time – a velocity that merely happens to vanish at t=0 is a Function and is never reported here.

Such a term contributes nothing to a form it is added to: UFL sums the integrands sharing a measure, so the zero collapses before FFCx sees it. It is used to warn the user, and to keep a form that would consist of nothing but zero integrands from being built at all – such a form compiles to one with no arguments, which cannot be assembled.

warn_if_no_effect(term: DriftTermBase, species: Species, velocity) → bool[source]#

Warn that term contributes nothing, and say why.

Writing a drift term is a statement of intent, so a term whose driving gradient vanishes is almost always an input that is uniform when the user meant it to vary. The term is still assembled – it costs nothing, and quietly discarding what the user asked for is worse than carrying a term that contributes zero – but it is worth saying so.

Callers that would otherwise build a form out of nothing but this term use the return value to skip it; see is_zero_velocity().

Parameters:

term – the drift term
species – the species it would have acted on
velocity – the velocity it produced

Returns:

True if the term has no effect

class CoupledTransientHeatTransferHydrogenTransport(heat_problem: HeatTransferProblem, hydrogen_problem: HydrogenTransportProblem)[source]#

Bases: object

Coupled heat transfer and hydrogen transport transient problem.

Parameters:

heat_problem – the heat transfer problem
hydrogen_problem – the hydrogen transport problem

Variables:

heat_problem (festim.heat_transfer_problem.HeatTransferProblem) – the heat transfer problem
hydrogen_problem (festim.hydrogen_transport_problem.HydrogenTransportProblem) – the hydrogen transport problem
non_matching_meshes (bool) – True if the meshes in the heat_problem and hydrogen_problem are not matching

Examples

import festim as F

my_heat_transfer_model = F.HeatTransferProblem(...)

my_h_transport_model = F.HydrogenTransportProblem(...)

coupled_problem = F.CoupledTransientHeatTransferHydrogenTransport(
    heat_problem=my_heat_transfer_model,
    hydrogen_problem=my_h_transport_model,
)

class HeatTransferProblem(mesh=None, subdomains=None, initial_condition=None, boundary_conditions=None, sources=None, exports=None, settings=None)[source]#

Bases: ProblemBase

create_dirichletbc_form(bc)[source]#

Creates a dirichlet boundary condition form.

Parameters:

bc (festim.FixedTemperatureBC) – the boundary condition

Returns:

A representation of: the boundary condition for modifying linear systems.

Return type:

dolfinx.fem.bcs.DirichletBC

create_flux_values_fenics()[source]#: For each heat flux create the value_fenics.

create_formulation()[source]#: Creates the formulation of the model.

create_initial_conditions()[source]#: For each initial condition, create the value_fenics and assign it to the previous solution of the condition’s species.

create_source_values_fenics()[source]#: For each source create the value_fenics.

define_function_space()[source]#

Creates the function space of the model, creates a mixed element if model is multispecies.

Creates the main solution and previous solution function u and u_n. Create global DG function spaces of degree 0 and 1 for the global diffusion coefficient

initialise_exports()[source]#: Defines the export writers of the model, if field is given as a string, find species object in self.species.

post_processing()[source]#: Post processes the model.

class Value(input_value, species_dependent_value: dict[str, Species] | None = None)[source]#

Bases: object

A class to handle input values from users and convert them to a relevent fenics object.

Parameters:

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.

Variables:

input_value (float | int | dolfinx.fem.function.Constant | dolfinx.fem.function.Expression | ufl.core.expr.Expr | dolfinx.fem.function.Function) – The value of the user input
species_dependent_value (dict[str, Species] | None) – 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 (dolfinx.fem.function.Expression) – The expression of the user input that is used to update the fenics_object
fenics_object (dolfinx.fem.function.Function | dolfinx.fem.function.Constant | ufl.core.expr.Expr) – The value of the user input in fenics format
explicit_time_dependent (bool) – True if the user input value is explicitly time dependent
temperature_dependent (bool) – True if the user input value is temperature dependent
species_dependent (bool) – True if the user input value depends on the concentration of other species

Converts a user given value to a relevent fenics object depending on the type of the value provided.

Parameters:

function_space – the function space of the fenics object, optional
t – the time, optional
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
mesh – the mesh the integral carrying this value is assembled over. Defaults to function_space.mesh. Pass it explicitly when the two differ – a term coupling a codim-1 subdomain to the bulk is a parent-mesh facet integral although its unknowns live on submeshes – so that ufl.SpatialCoordinate and any fem.Constant are built on the integration domain, not on the space of the unknown. Only meaningful with up_to_ufl_expr=True.
foreign_subdomain – where to read a species that does not live on subdomain. A term coupling a codim-1 subdomain to the bulk reaches into one particular side of it, and a bulk species defined on several subdomains has a solution on each; see solution_on(). Optional.

property explicit_time_dependent: bool#: Returns true if the value given is time dependent.

property species_dependent: bool#: Returns true if the value given depends on the concentration of other species.

property temperature_dependent: bool#: Returns true if the value given is temperature dependent.

update(t: float)[source]#

Updates the value.

Parameters:: t – the time

as_fenics_constant(value: float | int | Constant, mesh: Mesh) → Constant[source]#

Converts a value to a dolfinx.Constant.

Parameters:

value – the value to convert
mesh – the mesh of the domiain

Returns:

The converted value

Raises:

TypeError – if the value is not a float, an int or a dolfinx.Constant

as_fenics_interp_expr_and_function(value: Callable, function_space: FunctionSpace, t: Constant | None = None, temperature: Function | Constant | Expr | None = None, species_dependent_value: dict[str, Species] | None = None, subdomain: VolumeSubdomain | None = None) → tuple[Expression, Function][source]#

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.

Parameters:

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

Maps a user given callable function to the mesh, time, temperature or the concentration of other species within festim as needed.

Parameters:

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

is_it_time_to_export(times: list | None, current_time: float, atol=0, rtol=1e-05) → bool[source]#

Checks if the exported field should be written to a file or not based on the current time and the times in export.times

After a successful match, the corresponding time is removed from the list to prevent multiple exports for the same target time.

Parameters:

current_time – the current simulation time
atol – absolute tolerance for time comparison
rtol – relative tolerance for time comparison
times – the times at which the field should be exported, if None, returns True

Returns:

True if the exported field should be written to a file, else False

Return type:

bool

nmm_interpolate(f_out: Function, f_in: Function, cells: MeshTags | None = None, padding: float | None = 1e-11)[source]#

Non Matching Mesh Interpolate: interpolate one function (f_in) from one mesh into another function (f_out) with a mismatching mesh

Parameters:

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

restrict(expression, restriction: str | None)[source]#

Apply restriction to a whole expression.

Interior facet integrals need every discontinuous terminal restricted to one side of the facet. Restricting the expression as a whole rather than each terminal keeps a compound cross-mesh expression – k * (c_bulk - c_manifold), or a flux -D * grad(c) . n – on a single side, which also picks the outward normal of that side.

Parameters:

expression – the ufl expression to restrict
restriction – "+", "-", or None for an exterior facet integral, where nothing needs restricting

Returns:

The restricted expression, or expression unchanged when restriction is None

solution_on(species: Species, subdomain: VolumeSubdomain, foreign_subdomain: VolumeSubdomain | None = None)[source]#

The solution of species to use in an expression assembled for subdomain.

Usually the species lives on subdomain and this is just a dictionary lookup. In a codimensional coupling the expression deliberately reaches across meshes – a source on a manifold subdomain depending on the bulk concentration, say – and the solution to read is the one on the side of the manifold the term belongs to. The caller knows which side that is and passes it as foreign_subdomain; failing that, a species with a single solution has only one it could mean.

Parameters:

species – the species whose solution is needed
subdomain – the volume subdomain the expression is assembled for
foreign_subdomain – where to read a species that does not live on subdomain, optional

Returns:

The ufl expression of the species solution

Raises:

ValueError – if the species lives on several subdomains, none of which is subdomain or foreign_subdomain, so the choice would be arbitrary

class HydrogenTransportProblem(mesh: Mesh | None = None, subdomains: list[VolumeSubdomain | SurfaceSubdomain] | None = None, species: list[Species] | None = None, reactions: list[ReactionBase] | None = None, temperature: float | int | Constant | Function | Callable[[ndarray[Any, dtype[float64]]], ndarray[Any, dtype[float64]]] | Callable[[ndarray[Any, dtype[float64]], Constant], ndarray[Any, dtype[float64]]] | None = None, sources=None, initial_conditions=None, boundary_conditions=None, settings=None, exports=None, traps=None, advection_terms=None, drift_terms=None, petsc_options=None, element_immobile: str = 'CG')[source]#

Bases: ProblemBase

Hydrogen Transport Problem.

Parameters:

mesh – The mesh
subdomains – list containing the subdomains
species – list containing the species
reactions – list containing the reactions
temperature – The temperature or a function describing the temperature as a model of either space or space and time. Unit (K)
sources – The hydrogen sources
initial_conditions – The initial conditions
boundary_conditions – The boundary conditions
exports (list of festim.Export) – the exports of the model
traps (list of F.Trap) – the traps of the model
advection_terms – deprecated, use drift_terms. Appended to it
drift_terms – the drift terms of the model – advection, Soret, electromigration

Variables:

mesh (festim.mesh.mesh.Mesh) – The mesh
subdomains (list[festim.subdomain.volume_subdomain.VolumeSubdomain]) – The subdomains
species – The species
reactions – the reaction
temperature – The temperature in unit K
sources (list[festim.source.SourceBase]) – The hydrogen sources
initial_conditions – The initial conditions
boundary_conditions – list of Dirichlet boundary conditions
exports (list of festim.Export) – the export
traps (list of F.Trap) – the traps of the model
advection_terms – deprecated, use drift_terms. Appended to it
drift_terms (list[festim.drift.DriftTermBase]) – the drift terms of the model – advection, Soret, electromigration
dx (dolfinx.fem.dx) – the volume measure of the model
ds (dolfinx.fem.ds) – the surface measure of the model
function_space (dolfinx.fem.FunctionSpaceBase) – the function space of the model
facet_meshtags (dolfinx.mesh.MeshTags) – the facet meshtags of the model
volume_meshtags (dolfinx.mesh.MeshTags) – the volume meshtags of the model
formulation (ufl.form.Form) – the formulation of the model
solver (dolfinx.nls.newton.NewtonSolver) – the solver of the model
temperature_fenics (fem.Constant or fem.Function) – the temperature of the model as a fenics object (fem.Constant or fem.Function).
temperature_expr (fem.Expression) – the expression of the temperature that is used to update the temperature_fenics
temperature_time_dependent (bool) – True if the temperature is time dependent
V_DG_0 (dolfinx.fem.FunctionSpaceBase) – A DG function space of degree 0 over domain
V_DG_1 (dolfinx.fem.FunctionSpaceBase) – A DG function space of degree 1 over domain
volume_subdomains (list of festim.VolumeSubdomain) – the volume subdomains of the model
surface_subdomains (list of festim.SurfaceSubdomain) – the surface subdomains of the model

Examples

Can be used as either

import festim as F
my_model = F.HydrogenTransportProblem()
my_model.mesh = F.Mesh(...)
my_model.subdomains = [F.Subdomain(...)]
my_model.species = [F.Species(name="H"), F.Species(name="Trap")]
my_model.temperature = 500
my_model.sources = [F.ParticleSource(...)]
my_model.boundary_conditions = [F.BoundaryCondition(...)]
my_model.initialise()

or

my_model = F.HydrogenTransportProblem(
    mesh=F.Mesh(...),
    subdomains=[F.Subdomain(...)],
    species=[F.Species(name="H"), F.Species(name="Trap")],
)
my_model.initialise()

property advection_terms#

Deprecated. The advection terms among drift_terms.

An advection term is one kind of drift term, alongside festim.SoretTerm and festim.ElectromigrationTerm, and they all live in drift_terms now. This reads back only the festim.AdvectionTerm entries, so it is not an alias of the whole list.

assign_functions_to_species()[source]#: Creates the solution, prev solution, test function and post-processing solution for each species, as well as a collapsed function space for each species.

convert_drift_terms_to_fenics_objects()[source]#

For each drift term convert its user-given coefficients.

Runs after define_temperature so that a coefficient given as a function of T, or a drift velocity built from the temperature gradient, has one to read.

convert_reaction_rates_to_fenics_objects()[source]#: For each reaction convert its rate coefficients to fenics objects.

convert_source_input_values_to_fenics_objects()[source]#: For each source create the value_fenics.

create_dirichletbc_form(bc)[source]#

Creates a dirichlet boundary condition form.

Parameters:

bc (festim.DirichletBC) – the boundary condition

Returns:

A representation of: the boundary condition for modifying linear systems.

Return type:

dolfinx.fem.bcs.DirichletBC

create_flux_values_fenics()[source]#: For each particle flux create the value_fenics.

create_formulation()[source]#: Creates the formulation of the model.

create_implicit_species_value_fenics()[source]#: For each implicit species, create the value_fenics.

create_initial_conditions()[source]#: For each initial condition, create the value_fenics and assign it to the previous solution of the condition’s species.

create_sources_from_reactions()[source]#: Populate _unpacked_sources with the user-provided sources plus one volumetric particle source per species participating in each reaction.

create_species_from_traps()[source]#: Generate a species and reaction per trap defined in self.traps.

define_D_global(species)[source]#

Defines the global diffusion coefficient for a given species.

Parameters:

species (F.Species) – the species

Returns:

the global diffusion: coefficient and the expression of the global diffusion coefficient for a given species

Return type:

dolfinx.fem.Function, dolfinx.fem.Expression

define_boundary_conditions()[source]#: Defines the boundary conditions of the model.

define_function_spaces(element_degree: int = 1)[source]#

Creates the function space of the modelw with a mixed element. Creates the main solution and previous solution function u and u_n. Create global DG function spaces of degree 0 and 1 for the global diffusion coefficient.

Parameters:: element_degree – Degree order for finite element. Defaults to 1.

define_temperature()[source]#

Sets the value of temperature_fenics_value.

The type depends on self.temperature. If self.temperature is a function on t only, create a fem.Constant. Else, create an dolfinx.fem.Expression (stored in self.temperature_expr) to be updated, a dolfinx.fem.Function object is created from the Expression (stored in self.temperature_fenics_value). Raise a ValueError if temperature is None.

drift_velocity_in(field: Species, volume: VolumeSubdomain, temperature=None, mesh=None)[source]#

The summed drift velocity acting on field in volume.

Used in two places, both to do with the boundary term the divergence form leaves behind: festim.SurfaceFlux has to report it, and festim.OutflowBC cancels it.

Parameters:

field – the species whose flux is being computed
volume – the volume subdomain the flux leaves
temperature – the temperature on the mesh the integral is taken on, defaults to the problem’s
mesh – the mesh the integral is taken on, defaults to the problem’s

Returns:

the summed drift velocity as a ufl expression, or None when no drift term acts on field in volume

initialise_exports()[source]#: Defines the export writers of the model, if field is given as a string, find species object in self.species.

post_processing()[source]#: Post processes the model.

update_post_processing_solutions()[source]#: Updates the post-processing solutions of each species.

volume_subdomain_of_surface(surface: SurfaceSubdomain) → VolumeSubdomain[source]#

Returns the volume subdomain a surface subdomain belongs to.

The mapping is deduced from the mesh connectivity and the meshtags, and is computed once then cached.

Parameters:: surface – the surface subdomain
Returns:: the volume subdomain the surface belongs to
Raises:: ValueError – if the surface cannot be mapped to a volume subdomain

class HydrogenTransportProblemDiscontinuous(mesh=None, subdomains=None, species=None, reactions=None, temperature=None, sources=None, initial_conditions=None, boundary_conditions=None, settings=None, exports=None, traps=None, advection_terms=None, drift_terms=None, interfaces: list[Interface] | None = None, enclosures: list[Enclosure] | None = None, petsc_options: dict | None = None)[source]#

Bases: HydrogenTransportProblem

build_interior_facet_measure()[source]#

Builds the single dS measure that every interior-facet integral of the parent mesh shares, and stores it in interior_facet_measure.

The coupling terms of an interior manifold and the continuity terms of an festim.Interface are both dS integrals of the parent mesh, and both need integration data of their own so that "+" and "-" land on the sides they are meant to. They cannot each carry their own measure: UFL collects one subdomain_data entry per integral of a form, and dolfinx.fem.form asserts that they are all the same object before using the first of them for every id. A measure per manifold therefore breaks as soon as one volume subdomain touches two interior manifolds, or an interface and an interior manifold – both integrals end up in that subdomain’s F. So the data of every interior-facet integral goes into one list, and the measure built from it is handed to all of them, each indexing it by its own id.

convert_drift_terms_to_fenics_objects()[source]#

As the base class, but on the function space of the term’s own subdomain.

Every coefficient of a submesh integral has to be built on that submesh – FFCx cannot tabulate a parent-mesh coefficient on submesh cells – so the temperature is the subdomain’s own, as it is for reaction rates.

convert_enclosure_input_values_to_fenics_objects()[source]#: Converts the user input values of the enclosures and their openings to fenics objects.

convert_reaction_rates_to_fenics_objects()[source]#

As the base class, but with the temperature of the reaction’s own mesh.

A reaction on a manifold becomes a source integrated over that manifold’s submesh, so its rate cannot be built from the parent-mesh temperature: FFCx cannot tabulate a parent-mesh coefficient on submesh cells.

convert_source_input_values_to_fenics_objects()[source]#: For each source create the value_fenics.

coupling_measure(manifold: VolumeSubdomain, volume: VolumeSubdomain)[source]#: The measure, indexed, that the terms coupling manifold to volume are integrated over. Pair it with restriction_of().

coupling_measure_id(manifold: VolumeSubdomain, volume: VolumeSubdomain) → int[source]#

The id the terms coupling manifold to volume index the measure of facet_measure() with.

The manifold’s own id while it is adjacent to at most two volume subdomains – one integral carries both sides, told apart by their restriction. Beyond that each side has an integral of its own (see _allocate_manifold_side_ids()).

coupling_side(manifold: VolumeSubdomain, species: Species) → VolumeSubdomain[source]#

The volume subdomain on the side of manifold that species lives on.

This is how a coupling term declares which side of an interface it applies to: a flux boundary condition or a source names a bulk species, and that species identifies the side. Nothing else has to be specified by the user.

Raises:: ValueError – if the species does not live on exactly one of the manifold’s adjacent volume subdomains

create_dirichletbc_form(bc: FixedConcentrationBC)[source]#

Creates the value_fenics attribute for a given festim.FixedConcentrationBC and returns the appropriate dolfinx.fem.DirichletBC object.

Parameters:

bc (festim.FixedConcentrationBC) – the dirichlet BC

Returns:

the appropriate dolfinx representation: generated from dolfinx.fem.dirichletbc()

Return type:

dolfinx.fem.DirichletBC

create_enclosure_formulation(gas_species: GasSpecies)[source]#

Creates the variational formulation of the pressure balance of a gas species and stores it in gas_species.F.

Parameters:: gas_species – the gas species

create_flux_values_fenics()[source]#: For each particle flux create the value_fenics attribute.

create_formulation()[source]#

Takes all the formulations for each subdomain and adds the interface conditions.

Finally compute the jacobian matrix and store it in the J attribute, adds the entity_maps to the forms and store them in the forms attribute

create_implicit_species_value_fenics()[source]#

For each implicit species, create the value_fenics.

The density of an implicit species consumed by a reaction on a manifold subdomain appears in an integral over that manifold’s submesh, so like every other coefficient of such an integral it has to be built there rather than on the parent mesh (see create_submesh_time_constants()).

create_initial_conditions()[source]#: For each intial condition, create the value_fenics and assign it to the previous solution of the condition’s species.

create_solver()[source]#: Creates the solver of the model.

create_subdomain_formulation(subdomain: VolumeSubdomain)[source]#

Creates the variational formulation for each subdomain and stores it in subdomain.F and, for a manifold subdomain, subdomain.F_submesh.

Terms are split by whether they reach across meshes. Self terms – the time derivative, diffusion, advection, reactions and sources involving only fields of subdomain – are integrated over subdomain_measure(), which for a manifold subdomain is its own submesh. Coupling terms, which mix a manifold field with a bulk field, cannot live there (a bulk function cannot be resolved inside a codim-1 integral) and use facet_measure() on the parent mesh.

For a regular volume subdomain both measures are on the parent mesh and the whole formulation ends up in subdomain.F. For a manifold subdomain they are measures on two different meshes, which DOLFINx cannot compile into a single form, so the self terms are kept apart in subdomain.F_submesh and assembled as a separate group.

Parameters:: subdomain (F.VolumeSubdomain) – a subdomain of the geometry

create_submesh_time_constants()[source]#

Mirrors t and dt onto the submesh of every manifold subdomain.

The self terms of a manifold are integrated over its own submesh (see subdomain_measure()) and every coefficient of such an integral must live on that submesh: a fem.Constant bound to the parent mesh makes the integral mixed-dimensional and fails inside FFCx with an undiagnosable UnboundLocalError. The mirrors are kept in sync by update_submesh_time_constants().

define_boundary_conditions()[source]#: Defines the boundary conditions of the model.

define_enclosure_function_spaces()[source]#

Creates a real function space, a solution and a previous solution for each gas species of each enclosure.

The real function spaces live on the parent mesh: every form of the blocked system is integrated over the parent mesh, with submesh functions pulled in via entity_maps.

define_function_spaces(subdomain: VolumeSubdomain, element_degree=1)[source]#

Creates appropriate function space and functions for a given subdomain (submesh) based on the number of species existing in this subdomain. Then stores the functionspace, the current solution (u) and the previous solution (u_n) functions. It also populates the correspondance dicts attributes of the species (eg. species.subdomain_to_solution, species.subdomain_to_test_function, etc) for easy access to the right subfunctions, sub-testfunctions etc.

Parameters:

subdomain (F.VolumeSubdomain) – a subdomain of the geometry
element_degree (int, optional) – Degree order for finite element. Defaults to 1.

define_temperature()[source]#

Sets the value of temperature_fenics_value.

The type depends on self.temperature. If self.temperature is a function on t only, create a fem.Constant. Else, create an dolfinx.fem.Expression (stored in self.temperature_expr) to be updated, a dolfinx.fem.Function object is created from the Expression (stored in self.temperature_fenics_value). Raise a ValueError if temperature is None.

diffusion_coefficient(subdomain: VolumeSubdomain, species)[source]#: The diffusion coefficient of species for the gradient terms of subdomain, defined on the mesh those terms are integrated over.

export_surface_context(export: SurfaceQuantity)[source]#

Everything a surface quantity needs to know about where it is computed.

Three cases, and the ordinary one is unchanged:

an ordinary SurfaceSubdomain: the parent ds, and the volume subdomain it bounds;
a manifold (codim-1 VolumeSubdomain): the facets it occupies, read from the bulk side that export.field lives on. Interior manifolds use the dS coupling measure with that side’s restriction – the parent ds integrates to exactly zero over interior facets, which would be a silent zero rather than an error;
a codim-2 ``SurfaceSubdomain``: the boundary of the manifold that export.field lives on, integrated on that manifold’s submesh.

Returns:: (volume, mesh, measure, subdomain_id, restriction, entity_maps) where volume is the subdomain whose solution and material the quantity reads, mesh the one the integral is taken on – the parent mesh except on the boundary of a manifold – and subdomain_id the id to index measure with, which is the surface’s own except on a side of a manifold adjacent to more than two volume subdomains (see coupling_measure_id()).
Raises:: ValueError – if the volume subdomain given as a surface is not a manifold, or if a manifold’s own species is asked for on that manifold’s facets, which is not a flux across anything

export_volume_measure(volume: VolumeSubdomain)[source]#

The measure a volume quantity over volume is integrated with.

A manifold’s fields live on its own submesh, so its integrals are taken there – the parent dx carries no cell tagged with a manifold’s id, and assembling a submesh field against a parent measure fails inside FFCx anyway. The measure is given a meshtag covering the whole submesh so that the export can index it by the subdomain id, exactly as it does with the parent dx.

facet_measure(manifold: VolumeSubdomain)[source]#

The parent-mesh measure the facets of manifold are integrated in: the terms coupling it to the bulk, and the derived quantities exported on it.

A manifold on the boundary of the mesh is integrated with ds, one inside the mesh with the shared dS (see build_interior_facet_measure()), whose entities are ordered so that "+" is the first of the two volume subdomains it separates (see restriction_of()).

The measure is not restricted to the manifold: index it with manifold.id, as the derived quantities do with the parent ds. Left unindexed in a form it integrates over everything the measure carries data for – every tagged facet of the mesh for ds, every interior manifold and interface for the shared dS.

flux_bc_target(bc)[source]#

Where a ParticleFluxBC contributes: (volume subdomain, measure, restriction).

On an ordinary surface this is the single volume the surface belongs to. On a manifold subdomain it is the side bc.species lives on, which is how a user declares one flux per side of an interface without naming the side explicitly.

foreign_species(source, manifold: VolumeSubdomain) → list[Species][source]#

The species a source on manifold reads that do not live on it.

A non-empty result means the source is one half of a codimensional coupling rather than a source of the manifold’s own equation, and must therefore be integrated on the parent mesh.

Raises:: ValueError – if the source reads species from more than one side of the manifold, which cannot be expressed as a single restricted integral

gas_production_rates(surface, gas_species: GasSpecies)[source]#

The rates at which particles of a gas species are produced at a surface, in particles/s/m2, positive when entering the gas.

Parameters:

surface – the surface subdomain
gas_species – the gas species

Yields:

ufl expressions for each contribution at that surface

property gas_species#: All the gas species across all the enclosures. Defines the ordering of the pressure blocks in the solver.

initialise_exports()[source]#: Defines the export writers of the model, if field is given as a string, find species object in self.species.

interface_species(interface: Interface)[source]#

The mobile species whose continuity interface enforces.

An interface condition relates the two solutions of one species across the facets, so it applies to a species that has a solution on both of its volume subdomains. In a model with codim-1 subdomains most species do not: a manifold’s own species lives on a subdomain that is not a volume of any interface, and a bulk species may be confined to one side of a manifold. Those are simply not part of this interface’s condition.

A species present on exactly one of the two sides is the ambiguous case. It is either deliberately absent from the neighbouring material or a subdomain missing from its subdomains, and only the user can tell which, so it is skipped with a warning rather than silently.

is_manifold_self_source(source) → bool[source]#: Whether source belongs to a manifold’s own equation, as opposed to being one half of a codimensional coupling or an ordinary volumetric source.

iterate()[source]#: Iterates the model for a given time step.

link_enclosures()[source]#

Validates the enclosures and resolves the links between them, their surfaces and the boundary conditions coupled to them.

This runs before any fenics object is created.

manifold_is_interior(manifold: VolumeSubdomain) → bool[source]#

Whether manifold sits on interior facets of the mesh.

This is a question about the mesh, not about how many volume subdomains the manifold happens to separate: a grain boundary network inside a single-phase polycrystal is interior even though the same subdomain lies on both sides of it. Deciding from the subdomain count instead would pick ds for such a manifold, and ds integrates to exactly zero over interior facets – a coupling that silently does nothing.

Raises:: ValueError – if the manifold mixes interior and exterior facets, which would need two measures at once

manifold_of(surface: SurfaceSubdomain, species: Species) → VolumeSubdomain[source]#

The manifold subdomain that a codim-2 surface bounds.

A codim-2 surface is not in any meshtag, so it cannot be mapped to its volume topologically the way an ordinary surface is. It is instead resolved from the species of the term using it – the same rule that gives a flux on an interior manifold its side. That also means one surface object may be reused on several manifolds, one species each.

Raises:: ValueError – if the species does not live on exactly one manifold subdomain

outflow_form(bc, subdomain: VolumeSubdomain)[source]#

The contribution of an festim.OutflowBC to subdomain’s form.

Two positions. On an ordinary surface the integral is the parent ds of the volume the surface bounds. On the boundary of a manifold – the outlet of a 1D fluid, the case this exists for – everything lives on that manifold’s submesh: the measure, the normal, and the velocity, which was already built there when the drift term’s inputs were converted.

Parameters:

bc – the outflow boundary condition
subdomain – the volume subdomain whose form is being built

Returns:

the ufl form, or None if this bc does not belong to subdomain or no drift acts on its species there

post_processing()[source]#: Post processes the model.

restrict(expression, restriction: str | None)[source]#: Apply restriction to a whole expression, or return it unchanged when the coupling is on exterior facets.

restriction_of(manifold: VolumeSubdomain, volume: VolumeSubdomain) → str | None[source]#

Which side of an interior manifold the subdomain volume is on.

None when the manifold is on the boundary of the mesh, where the coupling is an exterior-facet integral and nothing needs restricting. "+" when the same subdomain lies on both sides: the bulk field is continuous across the facet, so both restrictions read the same value and the exchange is applied once. "+" again when the manifold is adjacent to more than two subdomains, where each side has an integral of its own on which it is the "+" one – so this must always be read together with coupling_measure_id().

run()[source]#: Runs the model.

source_coupling_side(source) → VolumeSubdomain | None[source]#

The bulk subdomain a coupling source on a manifold reads, or None if source is not one half of a codimensional coupling.

Which side of the manifold the term belongs to is decided by the bulk species the source names (see coupling_side()), and the expression has to read that species’ solution there. It cannot work that out for itself: a bulk species defined on several subdomains – because it is continuous across an interface, say – has a solution on each of them.

source_integration_mesh(source) → Mesh[source]#

The mesh the integral carrying source is assembled over.

A source involving only fields of a manifold subdomain is integrated over that manifold’s submesh (subdomain_measure()); one coupling the manifold to the bulk is a facet integral of the parent mesh (facet_measure()).

subdomain_dt(subdomain: VolumeSubdomain)[source]#: The timestep constant to use in the self terms of subdomain.

subdomain_measure(subdomain: VolumeSubdomain)[source]#

The measure the self terms of subdomain are integrated over.

Self terms are the ones that only involve fields living on subdomain: the time derivative, diffusion, advection, reactions and sources that do not reach across meshes. For a manifold (codim-1) subdomain they are integrated over its own submesh, which is what makes ufl.grad mean the tangential gradient (see the docs and issue #1208) and what makes the formulation identical whether the manifold sits on the boundary of the domain or between two volume subdomains.

Terms coupling a manifold to the bulk cannot use this measure – a bulk field cannot be resolved inside a codim-1 integral – and use facet_measure() instead.

subdomain_temperature(subdomain: VolumeSubdomain)[source]#: The temperature to use in the self terms of subdomain.

subdomain_time(subdomain: VolumeSubdomain)[source]#: The time constant to use in the self terms of subdomain.

update_submesh_time_constants()[source]#: Copies the current t and dt into the submesh mirrors created by create_submesh_time_constants(), so that an adaptive stepsize and explicitly time-dependent values on a manifold see the same values as the rest of the problem.

class InitialConcentration(value, volume, species: Species)[source]#

Bases: InitialConditionBase

Initial concentration class.

Parameters:

value – the value of the initial concentration of a given species.
species – the species to which the condition is applied
volume – the volume subdomain where the initial condition is applied

Variables:

value (float | int | dolfinx.fem.function.Constant | numpy.ndarray | dolfinx.fem.function.Expression | ufl.core.expr.Expr | dolfinx.fem.function.Function) – the value of the initial concentration of a given species.
species (festim.species.Species) – the species to which the condition is applied
volume (festim.subdomain.volume_subdomain.VolumeSubdomain) – the volume subdomain where the initial condition is applied
expr_fenics (collections.abc.Callable | dolfinx.fem.function.Expression) – the value of the initial condition in fenics expr format

Examples

InitialConcentration(value=1, species=my_species, volume=my_vol)
InitialConcentration(
    value=lambda x: 1 + x[0],
    species=my_species,
    volume=my_vol
)
InitialConcentration(
    value=lambda T: 1 + T,
    species=my_species,
    volume=my_vol
)
InitialConcentration(
    value=lambda x, T: 1 + x[0] + T,
    species=my_species,
    volume=my_vol
)

create_expr_fenics(mesh: Mesh, temperature: Function | Constant, function_space: functionspace)[source]#

Creates the expr_fenics of the initial condition.

If the value is a float or int, a function is created with an array with the shape of the mesh and all set to the value. Otherwise, it is converted to a fem.Expression.

Parameters:

mesh – the mesh
temperature – the temperature
function_space – the function space of the species

Bases: object

Base initial condition class.

Parameters:

value – the value of the initial condition.
volume – the volume subdomain where the initial condition is applied

Variables:

value (float | int | dolfinx.fem.function.Constant | numpy.ndarray | dolfinx.fem.function.Expression | ufl.core.expr.Expr | dolfinx.fem.function.Function) – the value of the initial condition.
volume (festim.subdomain.volume_subdomain.VolumeSubdomain) – the volume subdomain where the initial condition is applied

class InitialTemperature(value, volume)[source]#

Bases: InitialConditionBase

Initial temperature class.

Parameters:

value – the value of the initial temperature
volume – the volume subdomain where the initial condition is applied

Variables:

value (float | int | dolfinx.fem.function.Constant | numpy.ndarray | dolfinx.fem.function.Expression | ufl.core.expr.Expr | dolfinx.fem.function.Function) – the value of the initial temperature
volume (festim.subdomain.volume_subdomain.VolumeSubdomain) – the volume subdomain where the initial condition is applied
expr_fenics – the value of the initial condition in fenics expr format

Examples

InitialTemperature(value=1, volume=my_vol)
InitialTemperature(value=lambda x: 1 + x[0], volume=my_vol)
InitialTemperature(value=lambda x, t: 1 + x[0] + t, volume=my_vol)

create_expr_fenics(mesh: Mesh, function_space: functionspace)[source]#

Creates the expr_fenics of the initial condition.

If the value is a float or int, a function is created with an array with the shape of the mesh and all set to the value. Otherwise, it is converted to a fem.Expression.

Parameters:

mesh – the mesh
function_space – the function space of the species

read_function_from_file(filename: str, name: str, timestamp: int | float, family='P', order: int = 1, mesh: Mesh | None = None, backend: str = 'adios2') → Function[source]#

Read a function from a checkpoint file.

Reads a checkpoint written by an export with format="checkpoint" (or directly by io4dolfinx.write_function()). The visualisation formats ("vtx", "vtkhdf", "xdmf") store values interpolated onto the mesh nodes rather than the degrees of freedom, and cannot be read back with this function.

The function space the checkpoint is read into is built from family and order, so these must match the space the function was written from.

note::: The function is read from a file using io4dolfinx. For more information see the [io4dolfinx documentation]( scientificcomputing.github.io/io4dolfinx/README.html).

Parameters:

filename – the filename
name – the name of the function
timestamp – the timestamp of the function
family – the family of the function space
order – the order of the function space
mesh – Mesh to create input space on.
backend – the io4dolfinx backend the checkpoint was written with, "adios2" (.bp) or "h5py" (.h5). Must match the backend given to the export that wrote the file.

Returns:

the function

Raises:

ValueError – if backend is not a backend that can store a checkpoint

Bases: object

Material class.

Parameters:

D_0 – the pre-exponential factor of the diffusion coefficient (m2/s)
E_D – the activation energy of the diffusion coeficient (eV)
K_S_0 – the pre-exponential factor of the solubility coefficient (H/m3/Pa0.5)
E_K_S – the activation energy of the solubility coeficient (eV)
name – the name of the material
thermal_conductivity – the thermal conductivity of the material (W/m/K)
density – the density of the material (kg/m3)
heat_capacity – the heat capacity of the material (J/kg/K)
solubility_law – the solubility law of the material (SIEVERT, HENRY or NONE). For single material problems one can use NONE. This does not work for multi-material problems
D – the diffusion coefficient of the material (m2/s)

Variables:

D_0 – the pre-exponential factor of the diffusion coefficient (m2/s)
E_D – the activation energy of the diffusion coeficient (eV)
K_S_0 – the pre-exponential factor of the solubility coefficient (H/m3/Pa0.5)
E_K_S – the activation energy of the solubility coeficient (eV)
name – the name of the material
thermal_conductivity – the thermal conductivity of the material (W/m/K)
density – the density of the material (kg/m3)
heat_capacity – the heat capacity of the material (J/kg/K)
solubility_law – the solubility law of the material (SIEVERT, HENRY or NONE). For single material problems one can use NONE. This does not work for multi-material problems

Examples

# if only one species:
Material(D_0=1.9e-7, E_D=0.2, name="my_mat")

# if several species:
Material(
    D_0={"Species_1": 1.9e-7, "Species_2": 2.0e-7},
    E_D={"Species_1": 0.2, "Species_2": 0.3},
    name="my_mat"
)

get_D_0(species=None)[source]#

Returns the pre-exponential factor of the diffusion coefficient.

Parameters:: species (festim.Species or str, optional) – the species we want the pre-exponential factor of the diffusion coefficient of. Only needed if D_0 is a dict.
Returns:: the pre-exponential factor of the diffusion coefficient
Return type:: float

get_E_D(species=None)[source]#

Returns the activation energy of the diffusion coefficient.

Parameters:: species (festim.Species or str, optional) – the species we want the activation energy of the diffusion coefficient of. Only needed if E_D is a dict.
Returns:: the activation energy of the diffusion coefficient
Return type:: float

get_E_K_S(species=None) → float[source]#

Returns the activation energy of the solubility coefficient.

Parameters:: 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

get_K_S_0(species=None) → float[source]#

Returns the pre-exponential factor of the solubility coefficient.

Parameters:: 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

get_diffusion_coefficient(mesh=None, temperature=None, species=None)[source]#

Defines the diffusion coefficient.

Parameters:

mesh (dolfinx.mesh.Mesh) – the domain mesh
temperature (dolfinx.fem.Constant) – the temperature
species (festim.Species, optional) – the species we want the diffusion coefficient of. Only needed if D_0 and E_D are dicts.

Returns:

the diffusion coefficient

Return type:

ufl.algebra.Product

get_solubility_coefficient(mesh, temperature, species=None)[source]#

Defines the solubility coefficient.

Parameters:

mesh (dolfinx.mesh.Mesh) – the domain mesh
temperature (dolfinx.fem.Constant) – the temperature
species (festim.Species, optional) – the species we want the solubility coefficient of. Only needed if K_S_0 and E_K_S are dicts.

Returns:

the solubility coefficient

Return type:

ufl.algebra.Product

class SolubilityLaw(value)[source]#

Bases: Enum

classmethod from_string(s: str)[source]#: Can be removed with Python 3.11+.

class ProblemBase(mesh: Mesh = None, sources=None, exports=None, subdomains=None, boundary_conditions=None, settings=None, petsc_options=None)[source]#

Bases: object

Base class for HeatTransferProblem and HydrogenTransportProblem.

Variables:

show_progress_bar (bool) – If True a progress bar is displayed during the simulation
progress_bar (tqdm.asyncio.tqdm_asyncio | None) – the progress bar

create_solver()[source]#: Creates the solver of the model.

define_boundary_conditions()[source]#: Defines the dirichlet boundary conditions of the model.

define_meshtags_and_measures()[source]#: Defines the facet and volume meshtags of the model which are used to define the measures fo the model, dx and ds.

property facet_surface_subdomains#: The surface subdomains that bound an ordinary volume subdomain, ie. the ones living in the facet meshtags. These are all of them unless a surface bounds a manifold.

get_petsc_options() → dict[str, Any][source]#

Gets the PETSc options to pass to the NewtonProblem solver. Default options are updated with user-provided options, if any.

Returns:: the petsc options to pass to the NewtonProblem solver.

iterate()[source]#: Iterates the model for a given time step.

property manifold_boundary_subdomains#

the boundary of a manifold subdomain (the endpoints of a line, the rim of a surface).

They carry no meshtag – their entities are located on the manifold’s submesh when the boundary condition using them is created – so they take no part in the facet-tag id namespace and in surface_to_volume.

Type:: The codim-2 surface subdomains

property manifold_subdomains#

manifolds embedded in the mesh carrying their own transport equation. They are tagged in the facet meshtags and may be used wherever a surface subdomain is expected.

Type:: The codim-1 volume subdomains

run()[source]#: Runs the model.

Bases: GenericReaction

A reaction between species, with forward and backward rate coefficients built from Arrhenius laws. This is typically used to model trapping/detrapping.

Parameters:

reactant – The reactant(s).
product – The product(s). None for an irreversible reaction with no product.
k_0 – The forward rate constant pre-exponential factor.
E_k – The forward rate constant activation energy.
volume – The volume subdomain where the reaction takes place.
p_0 – The backward rate constant pre-exponential factor. Must be None when there is no product.
E_p – The backward rate constant activation energy. Must be None when there is no product.

Variables:

reactant (list[festim.species.Species | festim.species.ImplicitSpecies]) – The reactant(s).
product (list[festim.species.Species]) – The product(s), as a list (empty for an irreversible reaction).
k_0 (float) – The forward rate constant pre-exponential factor.
E_k (float) – The forward rate constant activation energy.
p_0 (float | None) – The backward rate constant pre-exponential factor.
E_p (float | None) – The backward rate constant activation energy.
volume (festim.subdomain.volume_subdomain.VolumeSubdomain) – The volume subdomain where the reaction takes place.

Examples

import festim as F

# create a volume subdomain for the reaction to take place in
material = F.Material(D_0=1, E_D=0)
volume = F.VolumeSubdomain(id=1, material=material)

# create two reactant species and a product species
reactant = [F.Species("A"), F.Species("B")]
product = F.Species("C")

# create a reaction between the reactants and the product
reaction = F.ArrheniusReaction(
    reactant=reactant,
    product=product,
    k_0=1.0,
    E_k=0.2,
    p_0=0.1,
    E_p=0.3,
    volume=volume,
)
print(reaction)

A + B <--> C

class DecayReaction(reactant: Species | ImplicitSpecies, half_life: float, volume: VolumeSubdomain, product: Species | list[Species] | None = None)[source]#

Bases: GenericReaction

A first-order radioactive decay reaction, consuming a single reactant species at a rate set by its decay half-life.

The net rate follows the mass-action form of GenericReaction

\[R = \lambda c\]

where the decay constant \(\lambda = \ln(2) / t_{1/2}\) is built from the half_life \(t_{1/2}\). The reaction is irreversible (no backward rate): the reactant is consumed at rate R and each product, if any, is produced at rate R (e.g. helium from the decay of tritium).

Parameters:

reactant – The decaying reactant species. Exactly one species decays: a decay is first order, so a list of more than one reactant is rejected.
half_life – The decay half-life, in the simulation’s time unit. Must be a positive float.
volume – The volume subdomain where the decay takes place.
product – The decay product(s). None if the products are not tracked.

Variables:

reactant (list[festim.species.Species | festim.species.ImplicitSpecies]) – The reactant, as a list of one species.
half_life (float) – The decay half-life.
product (festim.species.Species | list[festim.species.Species] | None) – The product(s), as a list (empty if none are tracked).
forward_rate (float | int | collections.abc.Callable | dolfinx.fem.function.Constant | dolfinx.fem.function.Expression | ufl.core.expr.Expr | dolfinx.fem.function.Function | festim.helpers.Value) – The decay constant \(\lambda\), as a festim.Value.
volume (festim.subdomain.volume_subdomain.VolumeSubdomain) – The volume subdomain where the decay takes place.

Examples

import festim as F

material = F.Material(D_0=1, E_D=0)
volume = F.VolumeSubdomain(id=1, material=material)

T = F.Species("T")  # tritium
He = F.Species("He")  # helium-3 produced by the decay

# tritium decays into helium with a half-life of ~12.3 years (in seconds)
reaction = F.DecayReaction(
    reactant=T,
    half_life=3.888e8,
    volume=volume,
    product=He,
)
print(reaction)

T --> He

Bases: ReactionBase

A reaction between one or more reactant species and zero or more product species, following a mass-action net rate

\[R = k_1 \prod_i c_i^{\text{reactant}} - k_2 \prod_j c_j^{\text{product}}\]

where \(k_1\) and \(k_2\) are the forward_rate and backward_rate coefficients. Both coefficients are festim.Value objects, so they can be a float, a ufl expression, or a callable of the temperature (argument T) and/or of other species concentrations (referenced by name through arg_to_species). This paves the way for building other reaction types (eg. hydride formation) on top of this class.

A reaction does not enter the formulation directly: it is expanded into one volumetric ParticleSource per participating Species (see create_sources()), each appearance in reactant consuming the species at rate R and each appearance in product producing it at rate R.

Parameters:

reactant – The reactant(s). A species listed twice appears squared in the mass-action rate and is consumed at rate 2 R (eg. 2 A --> B).
product – The product(s), following the same rule. None for an irreversible reaction with no product.
forward_rate – The forward reaction rate coefficient \(k_1\).
volume – The volume subdomain where the reaction takes place.
backward_rate – The backward reaction rate coefficient \(k_2\). If None, the reaction is irreversible.
arg_to_species – A dictionary mapping argument names in a callable rate coefficient to festim.Species objects, allowing a coefficient to depend on species concentrations. Every argument of the forward and backward coefficients other than the reserved t/x/T must appear as a key; extra keys that no coefficient uses are ignored with a warning. Every value must be a festim.Species. Alternatively the mapping may be attached directly to a rate coefficient passed as a festim.Value (its species_dependent_value), but the two ways are mutually exclusive: giving a mapping both here and on a rate Value raises a ValueError.

Variables:

reactant (list[festim.species.Species | festim.species.ImplicitSpecies]) – The reactant(s).
product (list[festim.species.Species]) – The product(s), as a list (empty for an irreversible reaction).
forward_rate (float | int | collections.abc.Callable | dolfinx.fem.function.Constant | dolfinx.fem.function.Expression | ufl.core.expr.Expr | dolfinx.fem.function.Function | festim.helpers.Value) – The forward reaction rate coefficient, as a festim.Value.
backward_rate (float | int | collections.abc.Callable | dolfinx.fem.function.Constant | dolfinx.fem.function.Expression | ufl.core.expr.Expr | dolfinx.fem.function.Function | festim.helpers.Value) – The backward reaction rate coefficient, as a festim.Value.
volume (festim.subdomain.volume_subdomain.VolumeSubdomain) – The volume subdomain where the reaction takes place.
arg_to_species (dict[str, festim.species.Species]) – The mapping used to resolve concentration arguments in a callable rate coefficient.

Examples

import festim as F

material = F.Material(D_0=1, E_D=0)
volume = F.VolumeSubdomain(id=1, material=material)

A = F.Species("A")
B = F.Species("B")
C = F.Species("C")

# net rate R = k1 * c_A * c_B - k2 * c_C, with k1 depending on T
reaction = F.GenericReaction(
    reactant=[A, B],
    product=C,
    forward_rate=lambda T: 1e-3 * T,
    backward_rate=2.0,
    volume=volume,
)
print(reaction)

A + B <--> C

property rate_coefficients: list[Value]#: The festim.Value rate coefficients that must be converted to fenics objects and updated in time. Subclasses with a different rate structure (e.g. GenericReaction) override this.

reaction_term(reactant_concentrations: list | None = None, product_concentrations: list | None = None) → Expr[source]#

Compute the net reaction rate R as a ufl expression.

The rate coefficients must already be converted to fenics objects (done by HydrogenTransportProblem.convert_reaction_rates_to_fenics_objects).

Parameters:

reactant_concentrations – The concentration to use for each reactant, same length as the reactants. Where an entry is None the reactant’s own concentration is used. If None, all reactant concentrations are used.
product_concentrations – The concentrations of the products, following the same rules as reactant_concentrations.

Returns:

The net reaction rate to be used in a formulation.

class Reaction(*args, **kwargs)[source]#

Bases: ArrheniusReaction

Deprecated alias for ArrheniusReaction.

Bases: object

A reaction that produces or consumes species at an arbitrary net rate, taking place within a volume.

The net rate is a single coefficient

\[R = f(c_i, x, T, t)\]

given by reaction_rate, a festim.Value that can be a float, a ufl expression, or a callable of the temperature (argument T), the spatial coordinate (x), the time (t) and/or other species concentrations (referenced by name through arg_to_species). R is whatever reaction_rate returns, so rates such as \(R = k (c_1 - c_2)\) are expressible.

A reaction does not enter the formulation directly: it is expanded into volumetric ParticleSource objects (see create_sources()), each reactant getting a sink -R and each product a source +R. The stoichiometry is carried by the reactant and product lists themselves: a species listed twice as a reactant is consumed at rate 2 R.

Parameters:

reaction_rate – The net reaction rate coefficient \(R\).
volume – The volume subdomain where the reaction takes place.
reactant – The reactant(s), consumed at rate R each. Implicit species have no governing equation and receive no source.
product – The product(s), produced at rate R each. None for a reaction with no product.
arg_to_species – A dictionary mapping argument names in a callable reaction_rate to festim.Species objects, allowing the rate to depend on species concentrations. Every argument of reaction_rate other than the reserved t/x/T must appear as a key; extra keys that the rate does not use are ignored with a warning. Every value must be a festim.Species. Alternatively the mapping may be attached directly to a reaction_rate passed as a festim.Value (its species_dependent_value), but the two ways are mutually exclusive: giving a mapping both here and on the rate Value raises a ValueError.

Variables:

reaction_rate (float | int | collections.abc.Callable | dolfinx.fem.function.Constant | dolfinx.fem.function.Expression | ufl.core.expr.Expr | dolfinx.fem.function.Function | festim.helpers.Value) – The net reaction rate coefficient, as a festim.Value.
volume (festim.subdomain.volume_subdomain.VolumeSubdomain) – The volume subdomain where the reaction takes place.
reactant (list[festim.species.Species | festim.species.ImplicitSpecies]) – The reactant(s), as a list.
product (list[festim.species.Species]) – The product(s), as a list (empty for a reaction with no product).
arg_to_species (dict[str, festim.species.Species]) – The mapping used to resolve concentration arguments in a callable reaction rate.

Examples

import festim as F

material = F.Material(D_0=1, E_D=0)
volume = F.VolumeSubdomain(id=1, material=material)

A = F.Species("A")
B = F.Species("B")

# an arbitrary rate R = 2 * (c_A - c_B): A is consumed and B produced
reaction = F.ReactionBase(
    reaction_rate=lambda c_A, c_B: 2.0 * (c_A - c_B),
    volume=volume,
    reactant=A,
    product=B,
    arg_to_species={"c_A": A, "c_B": B},
)
print(reaction)

A --> B

create_sources() → list[ParticleSource][source]#

Express the reaction as a list of volumetric particle sources, one per appearance of a Species in reactant or product.

Each reactant is given a sink -R and each product a source +R, where R is reaction_term(). The stoichiometry therefore comes from the lists themselves: a species listed twice gets two terms, which add up in the formulation. Implicit species (which have no governing equation) are skipped. The rate coefficients must already be converted to fenics objects.

Returns:: A list of festim.ParticleSource objects.

property rate_coefficients: list[Value]#: The festim.Value rate coefficients that must be converted to fenics objects and updated in time. Subclasses with a different rate structure (e.g. GenericReaction) override this.

reaction_term() → Expr[source]#

The net reaction rate R as a ufl expression.

reaction_rate must already be converted to a fenics object (done by HydrogenTransportProblem.convert_reaction_rates_to_fenics_objects).

Returns:: The net reaction rate to be unpacked into particle sources.

property species: list[Species | ImplicitSpecies]#

those receiving a source term and those the rate depends on. Used e.g. to update implicit-species densities.

Type:: All species involved in the reaction

class Settings(atol, rtol, max_iterations=30, transient=True, final_time=None, element_degree=1, stepsize=None, convergence_criterion: Literal['residual', 'incremental'] = 'residual')[source]#

Bases: object

Settings for a festim simulation.

Parameters:

atol (float or callable) – Absolute tolerance for the solver.
rtol (float or callable) – Relative tolerance for the solver.
max_iterations (int, optional) – Maximum number of iterations for the solver. Defaults to 30.
transient (bool, optional) – Whether the simulation is transient or not.
final_time (float, optional) – Final time for a transient simulation. Defaults to None
element_degree (int, optional) – Degree order for finite element. Defaults to 1.
stepsize (festim.Stepsize, optional) – stepsize for a transient simulation. Defaults to None
convergence_criterion – resiudal or incremental (for Newton solver)

Variables:

atol (float or callable) – Absolute tolerance for the solver.
rtol (float or callable) – Relative tolerance for the solver.
max_iterations (int) – Maximum number of iterations for the solver.
transient (bool) – Whether the simulation is transient or not.
final_time (float) – Final time for a transient simulation.
element_degree (int) – Degree order for finite element.
stepsize (festim.Stepsize) – stepsize for a transient simulation.
convergence_criterion – resiudal or incremental (for Newton solver)

class HeatSource(value, volume)[source]#

Bases: SourceBase

Heat source class.

Parameters:

value – the value of the source
volume – the volume subdomains where the source is applied

Variables:

value (float | int | dolfinx.fem.function.Constant | dolfinx.fem.function.Expression | ufl.core.expr.Expr | dolfinx.fem.function.Function) – the value of the source
volume (festim.subdomain.volume_subdomain.VolumeSubdomain) – the volume subdomains where the source is applied

Examples

from festim import HeatSource

HeatSource(volume=my_vol, value=1)
HeatSource(volume=my_vol, value=lambda x: 1 + x[0])
HeatSource(volume=my_vol, value=lambda t: 1 + t)
HeatSource(volume=my_vol, value=lambda x, t: 1 + x[0] + t)

class ParticleSource(value, volume, species: Species, species_dependent_value: dict[str, Species] | None = None)[source]#

Bases: SourceBase

Particle source class.

Parameters:

value – the value of the source
volume – the volume subdomains where the source is applied
species – the species to which the source is applied
species_dependent_value – a dictionary mapping the argument names in a callable value to festim.Species objects, allowing the source 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.

Variables:

value (float | int | dolfinx.fem.function.Constant | dolfinx.fem.function.Expression | ufl.core.expr.Expr | dolfinx.fem.function.Function) – the value of the source, as a festim.Value object
volume (festim.subdomain.volume_subdomain.VolumeSubdomain) – the volume subdomains where the source is applied
species (festim.species.Species) – the species to which the source is applied
species_dependent_value (dict[str, festim.species.Species] | None) – a dictionary mapping the argument names in a callable value to festim.Species objects, allowing the source 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.

Examples

from festim import ParticleSource

ParticleSource(volume=my_vol, value=1, species="H")
ParticleSource(volume=my_vol, value=lambda x: 1 + x[0], species="H")
ParticleSource(volume=my_vol, value=lambda t: 1 + t, species="H")
ParticleSource(volume=my_vol, value=lambda T: 1 + T, species="H")
ParticleSource(volume=my_vol, value=lambda x, t: 1 + x[0] + t, species="H")
ParticleSource(volume=my_vol, value=lambda c1: 2 * c1**2, species="H",
species_dependent_value={"c1": species1})

Bases: object

Source base class.

Parameters:

value – the value of the source
volume – the volume subdomains where the source is applied

Variables:

value (float | int | dolfinx.fem.function.Constant | dolfinx.fem.function.Expression | ufl.core.expr.Expr | dolfinx.fem.function.Function) – the value of the source
volume (festim.subdomain.volume_subdomain.VolumeSubdomain) – the volume subdomains where the source is applied

class ImplicitSpecies(n: float | callable, others: list[Species] | None = None, name: str | None = None)[source]#

Bases: object

Implicit species class for H transport simulation. c = n - others

Parameters:

n – the total concentration of the species
others – the list of species from which the implicit species concentration is computed (c = n - others)
name – a name given to the species. Defaults to None.

Variables:

n (float | callable) – the total concentration of the species
others (list[festim.species.Species] | None) – the list of species from which the implicit species concentration is computed (c = n - others)
name (str | None) – a name given to the species. Defaults to None.
concentration (ufl.form.Form) – the concentration of the species
value_fenics (dolfinx.fem.function.Constant | ufl.core.expr.Expr) – the total concentration as a fenics object

create_value_fenics(mesh, t: Constant)[source]#

Creates the value of the density as a fenics object and sets it to self.value_fenics. If the value is a constant, it is converted to a fenics.Constant. If the value is a function of t, it is converted to a fenics.Constant. Otherwise, it is converted to a ufl Expression.

Parameters:

mesh (dolfinx.mesh.Mesh) – the mesh
t (dolfinx.fem.Constant) – the time

update_density(t)[source]#

Updates the density value (only if the density is a function of time only)

Parameters:: t (float) – the time

class Species(name: str | None = None, mobile: bool = True, subdomains: list[VolumeSubdomain] | VolumeSubdomain | None = None)[source]#

Bases: object

Hydrogen species class for H transport simulation.

Parameters:

name – a name given to the species. Defaults to None.
mobile – whether the species is mobile or not. Defaults to True.
subdomain – the volume subdomain where the species is. Defaults to None.

Variables:

name (str | None) – a name given to the species.
mobile (bool) – whether the species is mobile or not.
solution (dolfinx.fem.function.Function | None) – the solution for the current timestep
prev_solution (dolfinx.fem.function.Function | None) – the solution for the previous timestep
test_function (ufl.argument.Argument | None) – the testfunction associated with this species
sub_function – the sub function of the species in case of multiple species in the same function space
sub_function_space (dolfinx.fem.function.FunctionSpace | None) – the subspace of the function space
collapsed_function_space (dolfinx.fem.function.FunctionSpace | None) – the collapsed function space for a species in the function space. In case single species case, this is None.
map_sub_to_main_solution (list | None) – the mapping from the sub solution dofs to the main solution dofs
post_processing_solution (dolfinx.fem.function.Function | None) – the solution for post processing
concentration (dolfinx.fem.function.Function | None) – the concentration of the species
subdomains (list[festim.subdomain.volume_subdomain.VolumeSubdomain] | festim.subdomain.volume_subdomain.VolumeSubdomain | None) – the volume subdomains where the species is
subdomain_to_solution (dict) – a dictionary mapping subdomains to solutions
subdomain_to_prev_solution (dict) – a dictionary mapping subdomains to previous solutions
subdomain_to_test_function (dict) – a dictionary mapping subdomains to test functions
subdomain_to_post_processing_solution (dict) – a dictionary mapping subdomains to post processing solutions
subdomain_to_collapsed_function_space (dict) – a dictionary mapping subdomains to collapsed function spaces
subdomain_to_function_space (dict) – a dictionary mapping subdomains to function spaces

Examples

:: testsetup:: Species

from festim import Species

:: testcode:: Species

Species(name=”H”) Species(name=”Trap”, mobile=False)

property legacy: bool#: Check if we are using FESTIM 1.0 implementation or FESTIM 2.0.

find_species_from_name(name: str, species: list)[source]#

Returns the correct species object from a list of species based on a string.

Parameters:

name (str) – the name of the species
species (list) – the list of species

Returns:

the species object with the correct name

Return type:

species (festim.Species)

Raises:

ValueError – if the species name is not found in the list of species

class Stepsize(initial_value, growth_factor=None, cutback_factor=None, target_nb_iterations=None, max_stepsize=None, milestones=None, milestone_tolerance=1e-05)[source]#

Bases: object

A class for evaluating the stepsize of transient simulations.

Parameters:

initial_value (float, int) – initial stepsize.
growth_factor (float, optional) – factor by which the stepsize is increased when adapting
cutback_factor (float, optional) – factor by which the stepsize is decreased when adapting
target_nb_iterations (int, optional) – number of Newton iterations over (resp. under) which the stepsize is increased (resp. decreased)
max_stepsize (float or callable, optional) – Maximum stepsize. If callable, has to be a function of t. Defaults to None.
milestones (list, optional) – list of times by which the simulation must pass. Defaults to an empty list.
milestone_tolerance (float, optional) – relative tolerance passed to numpy.isclose (rtol). Defaults to 1e-5.

Variables:

initial_value (float, int) – initial stepsize.
growth_factor (float) – factor by which the stepsize is increased when adapting
cutback_factor (float) – factor by which the stepsize is decreased when adapting
target_nb_iterations (int) – number of Newton iterations over (resp. under) which the stepsize is increased (resp. decreased)
adaptive (bool) – True if the stepsize is adaptive, False otherwise.
max_stepsize (float, callable) – Maximum stepsize.
milestones (list) – list of times by which the simulation must pass.
milestone_tolerance (float) – relative tolerance for how closely a time must align with a milestone to be triggered.

get_max_stepsize(t)[source]#

Returns the maximum stepsize at time t.

Parameters:: t (float) – the current time
Returns:: the maximum stepsize at time t
Return type:: float or None

is_adapt(t)[source]#

Methods that defines if the stepsize should be adapted or not.

Parameters:: t (float) – the current time
Returns:: True if needs to adapt, False otherwise.
Return type:: bool

next_milestone(current_time: float)[source]#

Returns the next milestone that the simulation must pass. Returns None if there are no more milestones.

Parameters:: current_time (float) – current time.
Returns:: next milestone.
Return type:: float

class Trap(name: str, mobile_species, k_0, E_k, p_0, E_p, n, volume)[source]#

Bases: Species

Trap species class for H transport simulation.

This class only works for 1 mobile species and 1 trapping level and is for convenience, for more details see notes.

Parameters:

name (str, optional) – a name given to the trap. Defaults to None.
mobile_species (_Species) – the mobile species to be trapped
k_0 (float) – the trapping rate constant pre-exponential factor (m3 s-1)
E_k (float) – the trapping rate constant activation energy (eV)
p_0 (float) – the detrapping rate constant pre-exponential factor (s-1)
E_p (float) – the detrapping rate constant activation energy (eV)
volume (F.VolumeSubdomain) – The volume subdomain where the trap is.

Variables:

name (str, optional) – a name given to the trap. Defaults to None.
mobile_species (_Species) – the mobile species to be trapped
k_0 (float) – the trapping rate constant pre-exponential factor (m3 s-1)
E_k (float) – the trapping rate constant activation energy (eV)
p_0 (float) – the detrapping rate constant pre-exponential factor (s-1)
E_p (float) – the detrapping rate constant activation energy (eV)
volume (F.VolumeSubdomain) – The volume subdomain where the trap is.
trapped_concentration (_Species) – The immobile trapped concentration
trap_reaction (_ArrheniusReaction) – The reaction for trapping the mobile conc.
empty_trap_sites (F.ImplicitSpecies) – The implicit species for the
sites (empty trap)

Examples

trap = Trap(name="Trap", mobile_species=H, k_0=1.0, E_k=0.2,
p_0=0.1, E_p=0.3, n=100, volume=my_vol)

my_model = HydrogenTransportProblem()
my_model.traps = [trap]

Notes

This convenience class replaces the need to specify an implicit species and the associated reaction, thus:

cm = _Species("mobile")
my_trap = F.Trap(
    name="trapped",
    mobile_species=cm,
    k_0=1,
    E_k=1,
    p_0=1,
    E_p=1,
    n=1,
    volume=my_vol,
)
my_model.species = [cm]
my_model.traps = [my_trap]

is equivalent to:

cm = _Species("mobile")
ct = _Species("trapped")
trap_sites = F.ImplicitSpecies(n=1, others=[ct])
trap_reaction = _ArrheniusReaction(
    reactant=[cm, trap_sites],
    product=ct,
    k_0=1,
    E_k=1,
    p_0=1,
    E_p=1,
    volume=my_vol,
)
my_model.species = [cm, ct]
my_model.reactions = [trap_reaction]

create_species_and_reaction()[source]#: Create the immobile trapped species object and the reaction for trapping.