Species and reactions

[RemDelaporteMathurin]

This is a proposal for a species & reaction architecture for the fenicsx rewrite.

The main idea is to define explicitely the species (unknowns of the problem) and reactions ( A + B <-> C). We would make use of the Species class and implement a Reaction class that takes reactants, product, and rate constants as arguments.

Basic case

A very simple example would be to do:

my_model = F.HTransportProblem()

cm = F.Species("cm")
ct = F.Species("ct")
empty_trap = F.Species("empty_trap")
my_model.species = [cm, ct, empty_trap]

my_model.reactions = [
    F.Reaction(
        reactant1=cm, reactant2=empty_trap, product=[ct], k_0=1, E_k=0.1, p_0=2, E_p=0.2
    ),
]

In the example above, cm reacts with empty_trap to create ct. A reaction is created and passed to the model.

The system of equations would be:

$$\frac{\partial c{m}}{\partial t} = \nabla \cdot D \nabla c{m} - k c{m} c{empty} + p c_{t}$$

$$\frac{\partial c{empty}}{\partial t} = -k c{m} c{empty} + p c{t}$$

$$\frac{\partial c{t}}{\partial t} = k c{m} c{empty} - p c{t}$$

This opens up the possibility to have diffusing trapped concentrations and/or diffusing empty traps (#376).

Implicit species

Now what if we don't want to explicitly define the empty trap but rather define it as $c_{empty} = n - c_t$ where $n=10$? Then we have to come with a new class ImplicitSpecies (we can find better naming...) that would, under the hood, do the c_empty = n - c_t thing.

my_model = F.HTransportProblem()

cm = F.Species("cm")
ct = F.Species("ct")

empty_trap = F.ImplicitSpecies("empty_trap", n=10, others=[ct])

my_model.species = [cm, ct]

my_model.reactions = [
    F.Reaction(
        reactant1=cm, reactant2=empty_trap, product=[ct], k_0=1, E_k=0.1, p_0=2, E_p=0.2
    ),
]

The system of equations is now:

$$ \frac{\partial c{m}}{\partial t} = \nabla \cdot D \nabla c{m} - k c{m} (n - c{t}) + p c_{t}$$

$$ \frac{\partial c{t}}{\partial t} = k c{m} (n - c{t}) - p c{t}$$

Multiple traps with implicit density

Same principle but extended to two traps:

my_model = F.HTransportProblem()

cm = F.Species("cm")
ct1 = F.Species("ct1")
ct2 = F.Species("ct2")

empty_trap1 = F.ImplicitSpecies("empty_trap1", n=10, others=[ct1])  # c_empty = n1 - ct1
empty_trap2 = F.ImplicitSpecies("empty_trap2", n=10, others=[ct2]) # c_empty = n2 - ct2

my_model.species = [cm, ct1, ct2]

my_model.reactions = [
    F.Reaction(
        reactant1=cm,
        reactant2=empty_trap1,
        product=[ct1],
        k_0=1,
        E_k=0.1,
        p_0=2,
        E_p=0.2,
    ),
    F.Reaction(
        reactant1=cm,
        reactant2=empty_trap2,
        product=[ct2],
        k_0=1,
        E_k=0.1,
        p_0=2,
        E_p=0.2,
    ),
]

The system of equations is:

$$\frac{\partial c{m}}{\partial t} = \nabla \cdot D \nabla c{m} - k{1} c{m} (n{1} - c{t1}) - k{2} c{m} (n{2} - c{t2}) + p{1} c{t1} + p{2} c{t2}$$

$$ \frac{\partial c{t1}}{\partial t} = k{1} c{m} (n{1} - c{t1}) - p{1} c_{t1}$$

$$ \frac{\partial c{t2}}{\partial t} = k{2} c{m} (n{2} - c{t2}) - p{2} c_{t2}$$

One trap, multiple levels

Now things start to get spicy... 🌶️♨️ With this architecture, it is possible to have multiple level trapping!

my_model = F.HTransportProblem()

cm = F.Species("cm")
ct_lev1 = F.Species("ct_lev1")
ct_lev2 = F.Species("ct_lev2")
ct_lev3 = F.Species("ct_lev3")
empty_trap = F.ImplicitSpecies("empty_trap", n=10, others=[ct_lev1, ct_lev2, ct_lev3])  # c_empty = n - ct_lev1 - ct_lev2 - ct_lev3

my_model.species = [cm, ct_lev1, ct_lev2, ct_lev3]

my_model.reactions = [
    F.Reaction(
        reactant1=cm,
        reactant2=empty_trap,
        product=[ct_lev1],
        k_0=1,
        E_k=0.1,
        p_0=2,
        E_p=0.2,
    ),
    F.Reaction(
        reactant1=cm,
        reactant2=ct_lev1,
        product=[ct_lev2],
        k_0=1,
        E_k=0.1,
        p_0=2,
        E_p=0.2,
    ),
    F.Reaction(
        reactant1=cm,
        reactant2=ct_lev2,
        product=[ct_lev3],
        k_0=1,
        E_k=0.1,
        p_0=2,
        E_p=0.2,
    ),
]

The system of equations is:

$$ c{empty} = n - c{t1} - c{t2} - c{t3} $$

$$\frac{\partial c{m}}{\partial t} = \nabla \cdot D \nabla c{m} - k{1} c{m} c{empty} + p{1} c{t1} - k{2} c{m} c{t1} + p{2} c{t2} - k{3} c{m} c{t2} + p{3} c_{t3}$$

$$\frac{\partial c{t1}}{\partial t} = k{1} c{m} c{empty} - p{1} c{t1} - k{2} c{m} c{t1} + p{2} c_{t2}$$

$$\frac{\partial c{t2}}{\partial t} = k{2} c{m} c{t1} - p{2} c{t2} - k{3} c{m} c{t2} + p{3} c_{t3}$$

$$\frac{\partial c{t3}}{\partial t} = k{3} c{m} c{t2} - p{3} c{t3}$$

2 diffusing species, 1 trap

Say we have two diffusing species that can react with one trap. #520 This architecture allows it very easily!

my_model = F.HTransportProblem()

cm1 = F.Species("cm1")  # ex: mobile H
cm2 = F.Species("cm2")  # ex: mobile D
ct1 = F.Species("ct1")  # ex: H in trap
ct2 = F.Species("ct2")  # ex: D in trap
empty_trap = F.ImplicitSpecies("empty_trap", n=10, others=[ct1, ct2])  # c_empty = n - ct1 - ct2

my_model.species = [cm1, cm2, ct1, ct2]

my_model.reactions = [
    F.Reaction(
        reactant1=cm1,
        reactant2=empty_trap,
        product=[ct1],
        k_0=1,
        E_k=0.1,
        p_0=2,
        E_p=0.2,
    ),
    F.Reaction(
        reactant1=cm2,
        reactant2=empty_trap,
        product=[ct2],
        k_0=1,
        E_k=0.1,
        p_0=2,
        E_p=0.2,
    ),
]

$$\frac{\partial c_{m1}}{\partial t} = \nabla \cdot D1 \nabla c{m1} - k{1} c{m1} (n - c{t1} - c{t2}) + p{1} c{t1}$$

$$ \frac{\partial c_{m2}}{\partial t} = \nabla \cdot D2 \nabla c{m2} - k{2} c{m2} (n - c{t1} - c{t2}) + p{2} c{t2}$$

$$ \frac{\partial c{t1}}{\partial t} = k{1} c{m1} (n - c{t1} - c{t2}) - p{1} c_{t1}$$

$$ \frac{\partial c{t2}}{\partial t} = k{2} c{m2} (n - c{t1} - c{t2}) - p{2} c_{t2}$$

Trap class for convenience

We want the simple use case to be simple. 1 trap, 1 mobile species.

my_model = F.HTransportProblem()

ct = F.Trap("ct", k_0=1, E_k=0.1, p_0=2, E_p=0.2, density=10)

my_model.species = [ct]

We would detect that:

• no mobile species was given, therefore create one under the hood
• a Trap object has been given so generate the reactions and trapped species from that.

The Trap class could be something like this:

class Trap:
    def __init__(self, name, k_0, E_k, p_0, E_p, density) -> None:
        self.name = name
        self.k_0 = k_0
        self.E_k = E_k
        self.p_0 = p_0
        self.E_p = E_p
        self.density = density

    def generate_reactions(self, mobile_species):
        self.trapped_conc = F.Species(f"c_{self.name}")

        self.implicit_density = F.ImplicitSpecies(
            f"empty_trap_{self.name}", n=self.density, others=[self.trapped_conc]
        )
        self.reaction = F.Reaction(
            reactant1=mobile_species,
            reactant2=self.implicit_density,
            product=[self.trapped_species[0]],
            k_0=self.k_0,
            E_k=self.E_k,
            p_0=self.p_0,
            E_p=self.E_p,
        )

We can also extend it to work for multi-level trapping and multiple mobile species, just a bit of logic.

Let's discuss!

@jhdark @stephen-dixon @gabriele-ferrero @SamueleMeschini @mougenj @ehodille

[mougenj]

Great programme ! I can't wait to test (and converge) it all.

A. Concerning the "One trap, multiple levels part", it's based on the formalism of a discrete variation of the trap energy as a function of filling. With 42 levels there would be almost as many reactions. Some authors, such as Ebihara et al [1], propose a continuous variation of the energy as a function of the filling rate (equation 5 in the paper), which makes it possible to model the multiple level without multiplying the number of variables. Would it be possible to have E_p not as a discrete value but as a function dependent on ct, for example?

B. For part "2 diffusing species". Here n is the trap density, or more excatly the number of diffusing objects that can enter the trap multiplie by the trap density; i.e. for the same trap density, n1 towels or n2 tea towels can come in. The filling capacity of the trap therefore depends on the size of the objects. If species 1 and 2 share the same trap, we could imagine that they do not contribute to filling the trap in the same way. Could we weight ct1 and ct2: in other words, instead of having "[ct1, ct2]" we could use "[0.5*ct1, ct2]"?

C. In this same section, if cm1 and cm2 are likely to be in the same trap, they are likely to use the same diffusion path. We could therefore imagine that after an exposure of cm1, which would fill the diffusion sites, the diffusion of cm2 would be modified. So could we foresee a development allowing us to write "grad(cm1+cm2)" instead of "grad(cm1)" in the first Fickian term?

D. Would all these developments make it possible to have a spatially distribution of empty_trap? i.e. with n=n(x)?

[1] Ebihara, K., Sugiyama, Y., Matsumoto, R., Takai, K. & Suzudo, T. Numerical Interpretation of Hydrogen Thermal Desorption Spectra for Iron with Hydrogen-Enhanced Strain-Induced Vacancies. Metallurgical and Materials Transactions A 52, 257–269 (2021).

[RemDelaporteMathurin]

Glad you like it 🚀

A. Concerning the "One trap, multiple levels part", it's based on the formalism of a discrete variation of the trap energy as a function of filling. With 42 levels there would be almost as many reactions. Some authors, such as Ebihara et al [1], propose a continuous variation of the energy as a function of the filling rate (equation 5 in the paper), which makes it possible to model the multiple level without multiplying the number of variables. Would it be possible to have E_p not as a discrete value but as a function dependent on ct, for example?

Yes the number of equations can quickly increase. I understand the idea of having a continuum for the detrapping energy. If you could write what the final system of equation would look like then we can try and have a look at this. I must admit I have a hard time seeing how c_t would give an information on the filling level? But surely writing the equations will help!

B. For part "2 diffusing species". Here n is the trap density, or more excatly the number of diffusing objects that can enter the trap multiplie by the trap density; i.e. for the same trap density, n1 towels or n2 tea towels can come in. The filling capacity of the trap therefore depends on the size of the objects. If species 1 and 2 share the same trap, we could imagine that they do not contribute to filling the trap in the same way. Could we weight ct1 and ct2: in other words, instead of having "[ct1, ct2]" we could use "[0.5*ct1, ct2]"?

n1 towels and n2 tea towels? 🤔🤔

In this example, the trap has only one filling level meaning that either mobile species 1 or mobile species 2 can be in the trap. Again, writing the system of equation would help understanding the behaviour

C. In this same section, if cm1 and cm2 are likely to be in the same trap, they are likely to use the same diffusion path. We could therefore imagine that after an exposure of cm1, which would fill the diffusion sites, the diffusion of cm2 would be modified. So could we foresee a development allowing us to write "grad(cm1+cm2)" instead of "grad(cm1)" in the first Fickian term?

@ehodille could talk to this but I'm fairly sure that one of the assumptions of this overall model is the "dilute approximation" meaning that the density of free diffusion site is much greater than the concentration. So very unlikely to "fill the diffusion sites".

I don't see right now what having grad(cm1+cm2) would do, although it should be possible to simply modify the formulation. We could also have concentration dependent diffusion coefficient.

D. Would all these developments make it possible to have a spatially distribution of empty_trap? i.e. with n=n(x)?

Yes. We haven't implemented spatially (and time) dependent trap densities but I guess we could:

Pass a function of x and t to ImplicitSpecies:

empty_trap = F.ImplicitSpecies("empty_traps", n=lambda x: x[0] + 2 * x[1], others=[ct])

or give an initial condition to Species if the empty traps are explicitly defined

or add a volumetric source to the empty trap species to emulate generation of empty traps (cc @jhdark )

[stephen-dixon]

Delighted to be included in the discussion :). Overall looks like it does everything required. Not explicitly listed in your examples but I think it would also work for the case of multiple isotopes in combination with multi-level trapping. The automation you mention at the end would be particularly important if that was a use-case of interest as the number of species/equations will grow quickly beyond what I'd personally have the patience to type out by hand...

I also think trap density as a function of space/time would be useful (and this was previously a feature in FESTIM-1 I thought?). Similarly, for the explicit trap-site case, is it easy to add source/sink terms for empty traps ? You could then inform trap distribution/evolution from e.g. a separate irradiation damage study for either implicit or explicit trap species. (possibly annealing too? I'm less informed on that case).

Another thought was on the Reaction class constructor. Could take a separate rate object (default case being Arhhenius) instead of the arhhenius parameters directly (is temperature missing or captured elsewhere?) in case you wanted to play with different rate models. Not sure if it's worth the hassle, just a thought.

[RemDelaporteMathurin]

Not explicitly listed in your examples but I think it would also work for the case of multiple isotopes in combination with multi-level trapping.

Yes! Here's the code for it

Code

```python my_model = F.HTransportProblem() c_H = F.Species("H") c_D = F.Species("D") ct_H1D0 = F.Species("ct_H1D0") ct_H2D0 = F.Species("ct_H2D0") ct_H0D1 = F.Species("ct_H0D1") ct_H0D2 = F.Species("ct_H0D2") ct_H1D1 = F.Species("ct_H1D1") trapped_species = [ ct_H1D0, ct_H2D0, ct_H0D1, ct_H0D2, ct_H1D1, ] empty_trap = F.ImplicitSpecies("empty_trap", n=10, others=trapped_species) my_model.species = [c_H, c_D] + trapped_species my_model.reactions = [ F.Reaction( reactant1=c_H, reactant2=empty_trap, product=[ct_H1D0], k_0=1, E_k=0.1, p_0=2, E_p=0.2, ), F.Reaction( reactant1=c_H, reactant2=ct_H1D0, product=[ct_H2D0], ... ), F.Reaction( reactant1=c_D, reactant2=empty_trap, product=[ct_H0D1], ... ), F.Reaction( reactant1=c_D, reactant2=ct_H0D1, product=[ct_H0D2], ... ), F.Reaction( reactant1=c_H, reactant2=ct_H0D1, product=[ct_H1D1], ... ), F.Reaction( reactant1=c_D, reactant2=ct_H1D0, product=[ct_H1D1], ... ), ] ```

And the equivalent reaction scheme:

$$ \frac{\partial c\mathrm{H}}{\partial t} = \nabla \cdot (D\mathrm{H} \nabla c_\mathrm{H}) \textcolor{red}{- k1 c\mathrm{H} c_{\mathrm{empty}} + p1 c{t1\mathrm{H}0\mathrm{D}}} \textcolor{blue}{- k2 c\mathrm{H} c_{t1\mathrm{H}0\mathrm{D}} + p2 c{t2\mathrm{H}0\mathrm{D}}} \textcolor{green}{- k5 c\mathrm{H} c_{t0\mathrm{H}1\mathrm{D}} + p5 c{t1\mathrm{H}1\mathrm{D}}}$$

$$\frac{\partial c\mathrm{D}}{\partial t} = \nabla \cdot (D\mathrm{D} \nabla c_\mathrm{D}) \textcolor{orange}{- k3 c\mathrm{D} c_{\mathrm{empty}} + p3 c{t0\mathrm{H}1\mathrm{D}}} \textcolor{purple}{- k4 c\mathrm{D} c_{t0\mathrm{H}1\mathrm{D}} + p4 c{t0\mathrm{H}2\mathrm{D}}} \textcolor{pink}{- k6 c\mathrm{D} c_{t1\mathrm{H}0\mathrm{D}} + p6 c{t1\mathrm{H}1\mathrm{D}}} $$

$$\frac{\partial c_{t1\mathrm{H}0\mathrm{D}}}{\partial t} = \textcolor{red}{k1 c\mathrm{H} c_{\mathrm{empty}} - p1 c{t1\mathrm{H}0\mathrm{D}}} \textcolor{blue}{- k2 c\mathrm{H} c_{t1\mathrm{H}0\mathrm{D}} + p2 c{t2\mathrm{H}0\mathrm{D}}} \textcolor{pink}{- k6 c\mathrm{D} c_{t1\mathrm{H}0\mathrm{D}} + p6 c{t1\mathrm{H}1\mathrm{D}}} $$

$$\frac{\partial c_{t2\mathrm{H}0\mathrm{D}}}{\partial t} = \textcolor{blue}{k2 c\mathrm{H} c_{t1\mathrm{H}0\mathrm{D}} - p2 c{t2\mathrm{H}0\mathrm{D}}} $$

$$\frac{\partial c_{t0\mathrm{H}1\mathrm{D}}}{\partial t} = \textcolor{orange}{k3 c\mathrm{D} c_{\mathrm{empty}} - p3 c{t0\mathrm{H}1\mathrm{D}}} \textcolor{purple}{- k4 c\mathrm{D} c_{t0\mathrm{H}1\mathrm{D}} + p4 c{t0\mathrm{H}2\mathrm{D}}} \textcolor{green}{- k5 c\mathrm{H} c_{t0\mathrm{H}1\mathrm{D}} + p5 c{t1\mathrm{H}1\mathrm{D}}} $$

$$\frac{\partial c_{t0\mathrm{H}2\mathrm{D}}}{\partial t} = \textcolor{purple}{k4 c\mathrm{D} c_{t0\mathrm{H}1\mathrm{D}} - p4 c{t0\mathrm{H}2\mathrm{D}}} \$$

$$\frac{\partial c_{t1\mathrm{H}1\mathrm{D}}}{\partial t} = \textcolor{pink}{k6 c\mathrm{D} c_{t1\mathrm{H}0\mathrm{D}} - p6 c{t1\mathrm{H}1\mathrm{D}}} \textcolor{green}{+ k5 c\mathrm{H} c_{t0\mathrm{H}1\mathrm{D}} - p5 c{t1\mathrm{H}1\mathrm{D}}} $$

the number of species/equations will grow quickly beyond what I'd personally have the patience to type out by hand

Yes x2.........

[RemDelaporteMathurin]

I also think trap density as a function of space/time would be useful (and this was previously a feature in FESTIM-1 I thought?).

We haven't implemented it yet but it will be needed for sure! And yes it is possible to do this in the current version of FESTIM.

Similarly, for the explicit trap-site case, is it easy to add source/sink terms for empty traps ?

We haven't implemented sources yet. But i reckon the usage would be:

empty_trap = F.Species("empty traps")
mobile = F.Species("H")
trapped = F.Species("trapped H")

my_model.species = [empty_trap, mobile, trapped]

trap_creation = F.Source(species=empty_trap, value=10, subdomain=...)
annealing = F.Source(species=empty_trap, value=..., subdomain=...)

my_model.sources = [trap_creation, annealing]

[RemDelaporteMathurin]

Another thought was on the Reaction class constructor. Could take a separate rate object (default case being Arhhenius) instead of the arhhenius parameters directly (is temperature missing or captured elsewhere?) in case you wanted to play with different rate models. Not sure if it's worth the hassle, just a thought.

Good shout, it has also been suggested at today's FESTIM meeting. Maybe we could have something like:

reaction = F.Reaction(...., k=lambda T: 2 * T, p=lambda T: 20 + T)

or

reaction = F.Reaction(...., k=2, p=4)

We would maybe keep the p_0, E_p behaviour as it's 99% of the use cases. Just need some logic to prevent users from providing both p_0, E_p and p.

[RemDelaporteMathurin]

This has been implemented. I'm closing this in favour of #653 to track progress for convenience classes