Radioactive decay

[RemDelaporteMathurin]

Proposed changes

This PR implements radioactive decay (fixes #488).

The new formulation is:

$$ \frac{\partial c_m}{\partial t} = \nabla \cdot (D \nabla c_m) - \lambda c_m - \sum \left( -k_i \ c_m \ (n_i - c_{t,i}) + p_i \ c_{t,i} \right) $$

$$ \frac{\partial c_{t,i}}{\partial t}=k_i \ c_m \ (n_i - c_{t,i})- p_i \ c_{t,i} - \lambda c_{t,i} $$

where $\lambda$ is the decay constant (expressed in s^-1).

It is implemented as a volumetric source object inheriting from F.Source.
It takes a decay constant parameter.

Usage:

import festim as F
import numpy as np


my_model = F.Simulation()

my_model.traps = [
    F.Trap(0, 1, 0, 3, materials=my_model.materials.materials[0], density=100)
]

half_life = 50  # seconds
decay_constant = np.log(2) / half_life

my_model.sources = [
    F.RadioactiveDecay(decay_constant=decay_constant, volume=1, field="solute"),
    F.RadioactiveDecay(decay_constant=decay_constant, volume=1, field="1"),
]

Alternatively, a single decay source can be applied to all concentrations (default behaviour):

import festim as F
import numpy as np


my_model = F.Simulation()

my_model.traps = [
    F.Trap(0, 1, 0, 3, materials=my_model.materials.materials[0], density=100)
]

half_life = 50  # seconds
decay_constant = np.log(2) / half_life

my_model.sources = [
    F.RadioactiveDecay(decay_constant=decay_constant, volume=1, field="all"),
]

On a simple simulation with no BCs, no trapping/detrapping, but initial conditions (200 for mobile and 100 for trapped concentration), a half-life of 50 s, things work as expected:

Types of changes

What types of changes does your code introduce to FESTIM?

• Bugfix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)
• Code refactoring
• Documentation Update (if none of the other choices apply)
• New tests

Checklist

• Black formatted
• Unit tests pass locally with my changes
• I have added tests that prove my fix is effective or that my feature works
• I have added necessary documentation (if appropriate)

Further comments

@jhdark I had to modify the formulation in order to make this work. Before, we had:

$$ \frac{\partial c_m}{\partial t} = \nabla \cdot (D \nabla c_m ) + S - \sum \frac{\partial c_{t,i}}{\partial t} $$

This trapping term wouldn't work when adding decay. Since:

$$ \frac{\partial c_{t,i}}{\partial t}=k_i \ c_m \ (n_i - c_{t,i})- p_i \ c_{t,i} - \lambda c_{t,i} $$

then the decay of trapped particles would create mobile particles!

I therefore had to:

• modify the formulation
• modify the MMS tests
• modify some other units tests

I'm happy to break this PR down and open a smaller PR that just does this

Also, it might be interesting to allow several materials per source!

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (7c6ca17) 98.93% compared to head (5093be4) 98.95%.
Report is 3 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #679      +/-   ##
==========================================
+ Coverage   98.93%   98.95%   +0.01%     
==========================================
  Files          57       58       +1     
  Lines        2154     2195      +41     
==========================================
+ Hits         2131     2172      +41     
  Misses         23       23              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[jhdark]

Just some minor things from me, and it seems there is still a little coverage missing. Otherwise looks good, could be very useful!

[jhdark]

looks good!

[RemDelaporteMathurin]

@Allentro You seemed interested in this feature. Do you have anything to add about this implementation?