Here’s a working example. I’ll try to walk you through the methodology:
Step 1: Write a custom trap class
We write a custom class called Hydride that is inheriting from Trap and we override the __init__ method (construction) and the create_trapping_form method with the correct trapping term.
import festim as F
import fenics as f
class Hydride(F.Trap):
def __init__(self, TSSd, materials):
self.TSSd = TSSd
super().__init__(k_0=0, E_k=0, p_0=0, E_p=0, materials=materials, density=0, id=None)
def create_trapping_form(self, mobile, materials, T, dx, dt=None):
"""
Args:
mobile (festim.Mobile): the mobile concentration of the simulation
materials (festim.Materials): the materials of the simulation
T (festim.Temperature): the temperature of the simulation
dx (fenics.Measure): the dx measure of the sim
dt (festim.Stepsize, optional): If None assuming steady state.
Defaults to None.
"""
solution = self.solution
prev_solution = self.previous_solution
test_function = self.test_function
if not all(isinstance(mat, F.Material) for mat in self.materials):
self.make_materials(materials)
expressions_trap = []
F_trapping = 0 # initialise the form
if dt is not None:
# d(c_t)/dt in trapping equation
F_trapping += ((solution - prev_solution) / dt.value) * test_function * dx
else:
# if the sim is steady state and
# if a trap is not defined in one subdomain
# add c_t = 0 to the form in this subdomain
for mat in materials:
if mat not in self.materials:
F_trapping += solution * test_function * dx(mat.id)
for i, mat in enumerate(self.materials):
if isinstance(mobile, F.Theta) and mat.solubility_law == "henry":
raise NotImplementedError(
"Henry law of solubility is not implemented with traps"
)
c_0, c_0_n = mobile.get_concentration_for_a_given_material(mat, T)
# ------- WE OVERRIDE THIS PART ---- >>>>>
F_trapping += -self.coupling_term(c_0) * test_function * dx(mat.id)
# ----------------------------------- <<<<<
self.F_trapping = F_trapping
self.F += self.F_trapping
self.sub_expressions += expressions_trap
def coupling_term(self, c_mobile):
# dc_hydride/dt = k_F * c_mobile - k_B
k_F = 1e16 / self.TSSd
k_B = 1e16 * f.exp(-1e-10/(self.solution + f.DOLFIN_EPS)) # added epsilon to avoid division by zero
return k_F * c_mobile - k_B
Step 2 monkey patching festim.Mobile.create_diffusion_form
A current limitation of FESTIM 1 is that users don’t explicitly define the mobile concentration instance + only the “classical” trapping form is considered.
Therefore we need to “mokey patch” the Mobile class to account for this new expression.
# Monkey patch the F.Mobile class (the create_diffusion_form method)
def create_diffusion_form(
self, materials, mesh, T, dt=None, traps=None, soret=False
):
"""Creates the variational formulation for the diffusive part.
Args:
materials (festim.Materials): the materials
mesh (festim.Mesh): the mesh
T (festim.Temperature): the temperature
dt (festim.Stepsize, optional): the stepsize. Defaults to None.
traps (festim.Traps, optional): the traps. Defaults to None.
chemical_pot (bool, optional): if True, conservation of chemical
potential is assumed. Defaults to False.
soret (bool, optional): If True, Soret effect is assumed. Defaults
to False.
"""
from fenics import dot, grad, exp, SpatialCoordinate
from festim import k_B
F = 0
for material in materials:
D_0 = material.D_0
E_D = material.E_D
c_0, c_0_n = self.get_concentration_for_a_given_material(material, T)
subdomains = material.id # list of subdomains with this material
if not isinstance(subdomains, list):
subdomains = [subdomains] # make sure subdomains is a list
# add to the formulation F for every subdomain
for subdomain in subdomains:
dx = mesh.dx(subdomain)
# transient form
if dt is not None:
F += ((c_0 - c_0_n) / dt.value) * self.test_function * dx
D = D_0 * exp(-E_D / k_B / T.T)
if mesh.type == "cartesian":
F += dot(D * grad(c_0), grad(self.test_function)) * dx
if soret:
Q = material.Q
if callable(Q):
Q = Q(T.T)
F += (
dot(
D * Q * c_0 / (k_B * T.T**2) * grad(T.T),
grad(self.test_function),
)
* dx
)
# see https://fenicsproject.discourse.group/t/method-of-manufactured-solution-cylindrical/7963
elif mesh.type == "cylindrical":
r = SpatialCoordinate(mesh.mesh)[0]
F += r * dot(D * grad(c_0), grad(self.test_function / r)) * dx
if soret:
Q = material.Q
if callable(Q):
Q = Q(T.T)
F += (
r
* dot(
D * Q * c_0 / (k_B * T.T**2) * grad(T.T),
grad(self.test_function / r),
)
* dx
)
elif mesh.type == "spherical":
r = SpatialCoordinate(mesh.mesh)[0]
F += (
D
* r
* r
* dot(grad(c_0), grad(self.test_function / r / r))
* dx
)
if soret:
Q = material.Q
if callable(Q):
Q = Q(T.T)
F += (
D
* r
* r
* dot(
Q * c_0 / (k_B * T.T**2) * grad(T.T),
grad(self.test_function / r / r),
)
* dx
)
# add the trapping terms
F_trapping = 0
if traps is not None:
for trap in traps:
for i, mat in enumerate(trap.materials):
if isinstance(trap, Hydride):
c_m, _ = self.get_concentration_for_a_given_material(mat, T)
F_trapping += -trap.coupling_term(c_m) * self.test_function * dx(mat.id)
else:
if type(trap.k_0) is list:
k_0 = trap.k_0[i]
E_k = trap.E_k[i]
p_0 = trap.p_0[i]
E_p = trap.E_p[i]
density = trap.density[i]
else:
k_0 = trap.k_0
E_k = trap.E_k
p_0 = trap.p_0
E_p = trap.E_p
density = trap.density[0]
c_m, _ = self.get_concentration_for_a_given_material(mat, T)
F_trapping += (
-k_0
* exp(-E_k / k_B / T.T)
* c_m
* (density - trap.solution)
* self.test_function
* dx(mat.id)
)
F_trapping += (
p_0
* exp(-E_p / k_B / T.T)
* trap.solution
* self.test_function
* dx(mat.id)
)
F += -F_trapping
self.F_diffusion = F
self.F += F
F.Mobile.create_diffusion_form = create_diffusion_form
Step 3 Integrate in a FESTIM model and plot results
Here, we use the newly created Hydride class in the traps of a festim model.
This is a simple 1D model with insulated boundaries and an initial concentration of mobile particles (equivalent to a 0D model).
import numpy as np
my_model = F.Simulation()
my_model.mesh = F.MeshFromVertices(np.linspace(0, 1e-3, 500))
my_model.materials = F.Material(id=1, D_0=1, E_D=0)
my_model.traps = [Hydride(TSSd=1e15, materials=my_model.materials)]
my_model.initial_conditions = [
F.InitialCondition(field="solute", value=1e16)
]
my_model.T = 600
my_model.settings = F.Settings(
absolute_tolerance=1e6,
relative_tolerance=1e-10,
final_time=1,
)
my_model.dt = F.Stepsize(initial_value=0.01)
c_mobile_export = F.AverageVolume(field="solute", volume=1)
c_hydride_export = F.AverageVolume(field=1, volume=1)
my_model.exports = [F.DerivedQuantities([c_mobile_export, c_hydride_export])]
# my_model.log_level = 20
my_model.initialise()
my_model.run()
# plot the results
import matplotlib.pyplot as plt
plt.plot(c_mobile_export.t, c_mobile_export.data, label="c_mobile")
plt.plot(c_hydride_export.t, c_hydride_export.data, label="c_hydride")
plt.legend()
plt.yscale("log")
plt.xlabel("Time (s)")
plt.ylabel("Concentration (m^-3)")
plt.show()
