Time dependent transport rates

libra_toolbox/tritium/model.py

@@ -5,6 +5,8 @@
 
 from libra_toolbox.tritium import ureg
 
+from typing import Callable
+
 SPECIFIC_ACT = 3.57e14 * ureg.Bq * ureg.g**-1
 MOLAR_MASS = 6.032 / 2 * ureg.g * ureg.mol**-1
 
@@ -70,8 +72,8 @@ def __init__(
         height: pint.Quantity,
         TBR: pint.Quantity,
         neutron_rate: pint.Quantity,
-        k_wall: pint.Quantity,
-        k_top: pint.Quantity,
+        k_wall: pint.Quantity | Callable[[pint.Quantity], pint.Quantity],
+        k_top: pint.Quantity | Callable[[pint.Quantity], pint.Quantity],
         irradiations: list,
     ) -> None:
 
@@ -154,7 +156,7 @@ def source(self, t):
                 return self.TBR * self.neutron_rate
         return 0 * self.TBR * self.neutron_rate
 
-    def Q_wall(self, c_salt):
+    def Q_wall(self, t, c_salt):
         """
         Calculate the release rate of tritium through the wall.
 
@@ -167,10 +169,13 @@ def Q_wall(self, c_salt):
         Returns:
             pint.Quantity: The release rate of tritium through the wall.
         """
+        if callable(self.k_wall):
+            k_wall = self.k_wall(t)
+        else:
+            k_wall = self.k_wall
+        return self.A_wall * k_wall * c_salt
 
-        return self.A_wall * self.k_wall * c_salt
-
-    def Q_top(self, c_salt):
+    def Q_top(self, t, c_salt):
         """
         Calculate the release rate of tritium through the top surface of the salt.
 
@@ -183,7 +188,11 @@ def Q_top(self, c_salt):
         Returns:
             pint.Quantity: The release rate of tritium through the top.
         """
-        return self.A_top * self.k_top * c_salt
+        if callable(self.k_top):
+            k_top = self.k_top(t)
+        else:
+            k_top = self.k_top
+        return self.A_top * k_top * c_salt
 
     def rhs(self, t, c):
         """
@@ -201,8 +210,8 @@ def rhs(self, t, c):
 
         return self.volume.to(ureg.m**3) ** -1 * (
             self.source(t).to(ureg.particle * ureg.s**-1)
-            - self.Q_wall(c).to(ureg.particle * ureg.s**-1)
-            - self.Q_top(c).to(ureg.particle * ureg.s**-1)
+            - self.Q_wall(t, c).to(ureg.particle * ureg.s**-1)
+            - self.Q_top(t, c).to(ureg.particle * ureg.s**-1)
         )
 
     def _generate_time_intervals(self, t_final):
@@ -222,7 +231,7 @@ def _generate_time_intervals(self, t_final):
 
         for irr in self.irradiations:
             t0 = irr[0]
-            tf = irr[1]
+            tf = min(irr[1], t_final)
 
             if previous_tf is not None:
                 time_intervals.append((previous_tf, t0))
@@ -287,7 +296,7 @@ def integrated_release_top(self):
             ndarray: array with units same size as the number of time steps,
             the integrated release of tritium through the top surface.
         """
-        top_release = self.Q_top(self.concentrations)
+        top_release = self.Q_top(t=self.times, c_salt=self.concentrations)
         integrated_top = cumulative_trapezoid(
             top_release.to(ureg.particle * ureg.h**-1).magnitude,
             self.times.to(ureg.h).magnitude,
@@ -304,7 +313,7 @@ def integrated_release_wall(self):
             ndarray: array with units same size as the number of time steps,
             the integrated release of tritium through the walls.
         """
-        wall_release = self.Q_wall(self.concentrations)
+        wall_release = self.Q_wall(t=self.times, c_salt=self.concentrations)
         integrated_wall = cumulative_trapezoid(
             wall_release.to(ureg.particle * ureg.h**-1).magnitude,
             self.times.to(ureg.h).magnitude,

libra_toolbox/tritium/plotting.py

@@ -285,14 +285,14 @@ def plot_salt_inventory(model: Model, **kwargs):
 
 
 def plot_top_release(model: Model, **kwargs):
-    top_release = model.Q_top(model.concentrations)
+    top_release = model.Q_top(t=model.times, c_salt=model.concentrations)
     top_release = quantity_to_activity(top_release).to(ureg.Bq * ureg.h**-1)
     l = plt.plot(model.times.to(ureg.day), top_release, **kwargs)
     return l
 
 
 def plot_wall_release(model: Model, **kwargs):
-    wall_release = model.Q_wall(model.concentrations)
+    wall_release = model.Q_wall(t=model.times, c_salt=model.concentrations)
     wall_release = quantity_to_activity(wall_release).to(ureg.Bq * ureg.h**-1)
     l = plt.plot(model.times.to(ureg.day), wall_release, **kwargs)
     return l

test/tritium/test_mass_conservation.py → test/tritium/test_model.py

@@ -32,3 +32,25 @@ def test_simple_case(TBR):
 
     assert np.isclose(model.concentrations[-1], 0)
     assert np.isclose(expected_total_production, computed_total_production, rtol=1e-2)
+
+
+def test_mass_transport_coeffs_time_dependent():
+
+    def k_top(t):
+        return 2 / (1 + t.magnitude) * ureg.m * ureg.s**-1
+
+    def k_wall(t):
+        return 3 / (1 + t.magnitude) * ureg.m * ureg.s**-1
+
+    model = Model(
+        radius=2 * ureg.m,
+        height=4 * ureg.m,
+        TBR=1 * ureg.particle * ureg.neutron**-1,
+        k_top=k_top,
+        k_wall=k_wall,
+        irradiations=[(0 * ureg.s, 10 * ureg.s), (60 * ureg.s, 70 * ureg.s)],
+        neutron_rate=30 * ureg.neutron * ureg.s**-1,
+    )
+
+    # run
+    model.run(100 * ureg.s)