Merge pull request #679 from festim-dev/radioactive-decay

Radioactive decay

docs/source/api/sources.rst

@@ -9,4 +9,8 @@ Sources
 
 .. autoclass:: ImplantationFlux
     :members:
-    :show-inheritance:
+    :show-inheritance:
+
+.. autoclass:: RadioactiveDecay
+    :members:
+    :show-inheritance:

festim/__init__.py

@@ -90,6 +90,7 @@
 
 from .sources.source import Source
 from .sources.source_implantation_flux import ImplantationFlux
+from .sources.radioactive_decay import RadioactiveDecay
 
 from .materials.material import Material
 from .materials.materials import Materials

festim/concentration/mobile.py

@@ -1,6 +1,5 @@
-from festim import Concentration, FluxBC, k_B, R
+from festim import Concentration, FluxBC, k_B, R, RadioactiveDecay
 from fenics import *
-import sympy as sp
 
 
 class Mobile(Concentration):
@@ -105,15 +104,41 @@ def create_diffusion_form(
                         * dx
                     )
 
-        # add the traps transient terms
-        if dt is not None:
-            if traps is not None:
-                for trap in traps.traps:
-                    F += (
-                        ((trap.solution - trap.previous_solution) / dt.value)
+        # add the trapping terms
+        F_trapping = 0
+        if traps is not None:
+            for trap in traps.traps:
+                for i, mat in enumerate(trap.materials):
+                    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
-                        * mesh.dx
+                        * dx(mat.id)
                     )
+        F += -F_trapping
+
         self.F_diffusion = F
         self.F += F
 
@@ -132,6 +157,9 @@ def create_source_form(self, dx):
                 volumes = source.volume
             else:
                 volumes = [source.volume]
+            if isinstance(source, RadioactiveDecay):
+                source.value = source.form(self.mobile_concentration())
+
             for volume in volumes:
                 F_source += -source.value * self.test_function * dx(volume)
             if isinstance(source.value, (Expression, UserExpression)):

festim/concentration/traps/trap.py

@@ -1,4 +1,4 @@
-from festim import Concentration, k_B, Material, Theta
+from festim import Concentration, k_B, Material, Theta, RadioactiveDecay
 from fenics import *
 import sympy as sp
 import numpy as np
@@ -204,6 +204,9 @@ def create_source_form(self, dx):
             dx (fenics.Measure): the dx measure of the sim
         """
         for source in self.sources:
+            if isinstance(source, RadioactiveDecay):
+                source.value = source.form(self.solution)
             self.F_source = -source.value * self.test_function * dx(source.volume)
             self.F += self.F_source
-            self.sub_expressions.append(source.value)
+            if isinstance(source.value, (Expression, UserExpression)):
+                self.sub_expressions.append(source.value)

festim/generic_simulation.py

@@ -188,7 +188,14 @@ def attribute_source_terms(self):
 
         # set sources
         for source in self.sources:
-            field_to_object[source.field].sources.append(source)
+            if isinstance(source, festim.RadioactiveDecay) and source.field == "all":
+                # assign source to each of the unique festim.Concentration
+                # objects in field_to_object
+                for obj in set(field_to_object.values()):
+                    if isinstance(obj, festim.Concentration):
+                        obj.sources.append(source)
+            else:
+                field_to_object[source.field].sources.append(source)
 
     def attribute_boundary_conditions(self):
         """Assigns boundary_conditions to mobile and T"""

festim/sources/radioactive_decay.py

@@ -0,0 +1,35 @@
+from festim import Source
+
+
+class RadioactiveDecay(Source):
+    """
+    A radioactive decay source.
+    dc/dt = ... - lambda * c
+    where lambda is the decay constant in 1/s.
+
+    Args:
+        decay_constant (float): The decay constant of the source in 1/s.
+        volume (int): The volume of the source.
+        field (str, optional): The field to which the source is
+            applied. If "all" the decay will be applied to all
+            concentrations. Defaults to "all".
+    """
+
+    def __init__(self, decay_constant, volume, field="all") -> None:
+        self.decay_constant = decay_constant
+        super().__init__(value=None, volume=volume, field=field)
+
+    @property
+    def decay_constant(self):
+        return self._decay_constant
+
+    @decay_constant.setter
+    def decay_constant(self, value):
+        if not isinstance(value, (float, int)):
+            raise TypeError("decay_constant must be a float or an int")
+        if value <= 0:
+            raise ValueError("decay_constant must be positive")
+        self._decay_constant = value
+
+    def form(self, concentration):
+        return -self.decay_constant * concentration

test/system/test_system.py

@@ -181,7 +181,8 @@ def test_run_MMS(tmpdir):
 
     f = (
         sp.diff(u, festim.t)
-        + sp.diff(v, festim.t)
+        - p * v
+        + k * u * (n_trap - v)
         - D * sp.diff(u, festim.x, 2)
         - sp.diff(D, festim.x) * sp.diff(u, festim.x)
     )
@@ -303,7 +304,8 @@ def test_run_MMS_chemical_pot(tmpdir):
 
     f = (
         sp.diff(u, festim.t)
-        + sp.diff(v, festim.t)
+        - p * v
+        + k * u * (n_trap - v)
         - D * sp.diff(u, festim.x, 2)
         - sp.diff(D, festim.x) * sp.diff(u, festim.x)
     )
@@ -589,12 +591,11 @@ def test_run_MMS_steady_state(tmpdir):
     k = k_0 * sp.exp(-E_k / k_B / T)
 
     f = (
-        sp.diff(u, festim.t)
-        + sp.diff(v, festim.t)
-        - D * sp.diff(u, festim.x, 2)
+        -D * sp.diff(u, festim.x, 2)
         - sp.diff(D, festim.x) * sp.diff(u, festim.x)
+        - (p * v - k * u * (n_trap - v))
     )
-    g = sp.diff(v, festim.t) + p * v - k * u * (n_trap - v)
+    g = p * v - k * u * (n_trap - v)
 
     def run(h):
         my_materials = festim.Materials(
@@ -1031,3 +1032,171 @@ def test_completion_tone():
     )
     my_model.initialise()
     my_model.run(completion_tone=True)
+
+
+def test_mms_radioactive_decay():
+    """MMS test for radioactive decay
+    Steady state, only solute
+    """
+    u = 1 + sp.sin(2 * fenics.pi * festim.x)
+    size = 1
+    T = 700 + 30 * festim.x
+    E_D = 0.1
+    D_0 = 2
+    decay_constant = 0.1
+    k_B = festim.k_B
+    D = D_0 * sp.exp(-E_D / k_B / T)
+
+    f = (
+        -D * sp.diff(u, festim.x, 2)
+        - sp.diff(D, festim.x) * sp.diff(u, festim.x)
+        + decay_constant * u
+    )
+
+    my_materials = festim.Material(name="mat", id=1, D_0=D_0, E_D=E_D)
+
+    my_initial_conditions = [
+        festim.InitialCondition(field=0, value=u),
+    ]
+
+    my_mesh = festim.MeshFromVertices(np.linspace(0, size, 1000))
+
+    my_bcs = [
+        festim.DirichletBC(surfaces=[1, 2], value=u, field=0),
+    ]
+
+    my_temp = festim.Temperature(T)
+
+    my_sources = [
+        festim.Source(f, 1, "0"),  # MMS source term
+        festim.RadioactiveDecay(decay_constant=decay_constant, volume=1),
+    ]
+
+    my_settings = festim.Settings(
+        absolute_tolerance=1e-10,
+        relative_tolerance=1e-9,
+        maximum_iterations=50,
+        transient=False,
+    )
+
+    my_sim = festim.Simulation(
+        mesh=my_mesh,
+        materials=my_materials,
+        initial_conditions=my_initial_conditions,
+        boundary_conditions=my_bcs,
+        temperature=my_temp,
+        sources=my_sources,
+        settings=my_settings,
+    )
+
+    my_sim.initialise()
+    my_sim.run()
+    error_max_u = compute_error(
+        u,
+        computed=my_sim.mobile.post_processing_solution,
+        t=0,
+        norm="error_max",
+    )
+
+    tol_u = 1e-7
+    msg = f"Maximum error on u is: {error_max_u}"
+    print(msg)
+    assert error_max_u < tol_u
+
+
+def test_MMS_decay_with_trap():
+    """MMS test for radioactive decay
+    Steady state, solute and trap
+    """
+    u = 1 + festim.x
+    v = 1 + festim.x * 2
+    size = 1
+    k_0 = 2
+    E_k = 1.5
+    p_0 = 0.2
+    E_p = 0.1
+    T = 700 + 30 * festim.x
+    decay_constant = 0.1
+    n_trap = 1
+    E_D = 0.1
+    D_0 = 2
+    k_B = festim.k_B
+    D = D_0 * sp.exp(-E_D / k_B / T)
+    p = p_0 * sp.exp(-E_p / k_B / T)
+    k = k_0 * sp.exp(-E_k / k_B / T)
+
+    f = (
+        -p * v
+        + k * u * (n_trap - v)
+        - D * sp.diff(u, festim.x, 2)
+        - sp.diff(D, festim.x) * sp.diff(u, festim.x)
+        + decay_constant * u
+    )
+    g = p * v - k * u * (n_trap - v) + decay_constant * v
+
+    my_materials = festim.Materials(
+        [festim.Material(name="mat", id=1, D_0=D_0, E_D=E_D)]
+    )
+
+    my_trap = festim.Trap(k_0, E_k, p_0, E_p, ["mat"], n_trap)
+
+    size = 0.1
+    my_mesh = festim.MeshFromVertices(np.linspace(0, size, 1600))
+
+    my_sources = [
+        festim.Source(f, 1, "solute"),
+        festim.Source(g, 1, "1"),
+        festim.RadioactiveDecay(decay_constant=decay_constant, volume=1),
+    ]
+
+    my_temp = festim.Temperature(T)
+
+    my_bcs = [
+        festim.DirichletBC(surfaces=[1, 2], value=u, field=0),
+        festim.DirichletBC(surfaces=[1, 2], value=v, field=1),
+    ]
+
+    my_settings = festim.Settings(
+        absolute_tolerance=1e-10,
+        relative_tolerance=1e-9,
+        maximum_iterations=50,
+        transient=False,
+        traps_element_type="DG",
+    )
+
+    my_sim = festim.Simulation(
+        mesh=my_mesh,
+        materials=my_materials,
+        traps=my_trap,
+        boundary_conditions=my_bcs,
+        sources=my_sources,
+        temperature=my_temp,
+        settings=my_settings,
+    )
+
+    my_sim.initialise()
+    my_sim.run()
+    error_max_u = compute_error(
+        u,
+        computed=my_sim.mobile.post_processing_solution,
+        t=my_sim.t,
+        norm="error_max",
+    )
+    error_max_v = compute_error(
+        v,
+        computed=my_sim.traps.traps[0].post_processing_solution,
+        t=my_sim.t,
+        norm="error_max",
+    )
+
+    tol_u = 1e-10
+    tol_v = 1e-7
+    msg = (
+        "Maximum error on u is:"
+        + str(error_max_u)
+        + "\n \
+        Maximum error on v is:"
+        + str(error_max_v)
+    )
+    print(msg)
+    assert error_max_u < tol_u and error_max_v < tol_v