begin
    using SciMLBase: ODEProblem, solve
    using OrdinaryDiffEqLowOrderRK: DP5
    using OrdinaryDiffEqRosenbrock: Rosenbrock23
    using OrdinaryDiffEqTsit5:Tsit5
    using Catalyst
    using VoronoiFVM: VoronoiFVM, enable_species!, enable_boundary_species!
    using VoronoiFVM: ramp, boundary_dirichlet!
    using ExtendableGrids: simplexgrid
    using GridVisualize: GridVisualize, GridVisualizer, reveal, scalarplot!, gridplot, available_kwargs
    if doplots # defined in the Appendix
      using Plots: Plots, plot, theme
      using PlotThemes
      Plots.gr()
      Plots.theme(:dark)
      GridVisualize.default_plotter!(Plots)
    end
    import PlutoUI
    using Test
    PlutoUI.TableOfContents(; depth = 4)
end

Towards Heterogeneous Catalysis

How to model and simulate heterogeneous catalysis with Catalyst.jl and VoronoiFVM.jl.

Mass action kinetics

Sources:

Assume $j=1 … M$ reversible reactions of educt (substrate) species to product species

$$ α_1^j S_1 + α_2^j S_2 + … + α_n^jS_n \underset{k_j^+}{\stackrel{k_j^-}{\longrightleftharpoons}} β_1^jS_1 + β_2^j S_2 + … + β_n^j S_n$$

or equivalently,

$$∑_{i=1}^n α_{i}^j S_i \underset{k_j^+}{\stackrel{k_j^-}{\longrightleftharpoons}} ∑_{i=1}^n β_i^j S_i $$

The rate of these reactions depend on the concentrations $[S_i]$ of the species. Within $Catalyst.jl$, due to consistency with the derivation from stochastic approaches, the default "combinatoric rate law" is

$$ r_j=k_j^+ ∏_{i=1}^n \frac{[S_i]^{α_i^j}}{α_i^j!} - k_j^- ∏_{i=1}^n \frac{[S_i]^{β_i^j}}{β_i^j!} $$

while it appears that in most textbooks, the rate law is

$$ r_j=k_j^+ ∏_{i=1}^n [S_i]^{α_i^j} - k_j^- ∏_{i=1}^n [S_i]^{β_i^j}.$$

See the corresponding remark in the Catalyst.jl docs and the github issue. We will stick to the secobd version which can be achieved by setting combinatoric_ratelaws to false at appropriate places.

Later in the project we will see that instead of the concentrations, we need to work with so called activities.

The numbers $σ_i^j=α_i^j-β_i^j$ are the net stoichiometric coefficients of the system.

The rate differential equations then are (TODO: check this)

$$ ∂_t [S_i] + \sum_{j=1}^M \sigma_i^jr_j = 0$$

These assume that the reactions take place in a continuously stirred tank reactor (CSTR) which means that we have complete mixing, and species concentrations are spatially constant, and we have just one concentration value for each species at a given point of time.

Example 1: A $\longrightleftharpoons$ B

$$\begin{aligned} A& \underset{k_1^+}{\stackrel{k_1^-}{\longrightleftharpoons}} B\\ r_1&= k_1^+ [A] - k_1^- [B]\\ \partial_t[A] &= -r_1\\ \partial_t[B] &= r_1 \end{aligned} $$

Solution via plain ODE problem using OrdinaryDiffEq.jl:

Set the parameters such that the forward reaction is faster then the backward reaction:

p1 = (k_p = 1, k_m = 0.1)

(k_p = 1, k_m = 0.1)

Define an ODE function describing the right hand side of the ODE System:

function example1(du, u, p, t)
    (; k_p, k_m) = p
    r1 = k_p * u[1] - k_m * u[2]
    du[1] = -r1
    du[2] = r1
end

example1 (generic function with 1 method)

Define some initial value:

u1_ini = [1.0, 0.0]

2-element Vector{Float64}:
 1.0
 0.0

Create an solve the corresponding ODEProblem

prob1 = ODEProblem(example1, u1_ini, (0, 10), p1)

ODEProblem with uType Vector{Float64} and tType Int64. In-place: true
timespan: (0, 10)
u0: 2-element Vector{Float64}:
 1.0
 0.0

sol1 = solve(prob1, DP5())

	timestamp	value1	value2
1	0.0	1.0	0.0
2	0.000999001	0.999002	0.000998452
3	0.010989	0.989077	0.0109229
4	0.110889	0.895607	0.104393
5	0.418818	0.664402	0.335598
6	0.859981	0.443912	0.556088
7	1.47125	0.271123	0.728877
8	2.18013	0.173557	0.826443
9	2.97759	0.125302	0.874698
10	3.85357	0.104043	0.895957
11	4.83614	0.0953755	0.904625
12	5.96713	0.0922044	0.907796
13	7.31295	0.091211	0.908789
14	8.97033	0.0909642	0.909036
15	10.0	0.0909269	0.909073

doplots && plot(sol1; size = (600, 200))

Mass conservation: adding the two reaction eqauations results in

$$ \partial_t ([A]+[B]) = 0,$$

therefore $[A]+[B]$ must be constant:

all(s -> isapprox(s, sum(u1_ini)), sum(sol1; dims = 1))

true

`Catalyst.@reaction_network`

Catalyst.jl provides a convenient way to define a reaction network, and the resulting reaction system. So we use this to build the same system:

rn1 = @reaction_network rn1 begin
    @combinatoric_ratelaws false
    k_p, A --> B
    k_m, B --> A
end

$$\begin{align*} \mathrm{A} &\xrightleftharpoons[k_{m}]{k_{p}} \mathrm{B} \end{align*}$$

The corresponding ODE system is:

convert(ODESystem, rn1)

$$\begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= k_{m} B\left( t \right) - k_{p} A\left( t \right) \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= - k_{m} B\left( t \right) + k_{p} A\left( t \right) \end{align}$$

Catalyst.jl adds a new method to the ODEProblem constructor which allows to pass a reaction nerwork:

prob1n = ODEProblem(rn1, u1_ini, (0, 10.0), Dict(pairs(p1)))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
 1.0
 0.0

sol1n = solve(prob1n, Rosenbrock23())

	timestamp	A(t)	B(t)
1	0.0	1.0	0.0
2	0.000113386	0.999887	0.000113379
3	0.00124725	0.998754	0.0012464
4	0.00810266	0.991933	0.00806667
5	0.0183376	0.981846	0.0181539
6	0.0399115	0.960951	0.0390486
7	0.0695373	0.933054	0.066946
8	0.117184	0.890048	0.109952
9	0.177898	0.838412	0.161588
10	0.261426	0.772769	0.227231
...

doplots && plot(sol1n; idxs = [rn1.A, rn1.B, rn1.A + rn1.B], size = (600, 200))

Unraveling `@reaction_network`

Let us try to look behind the macro voodoo - this is necessary to build networks programmatically and is behind the translation from python to Julia in CatmapInterface.jl.

It is interesting if there is a "macro - less" way to define variables, parameters and species.

@variables t

$$\begin{equation} \left[ \begin{array}{c} t \\ \end{array} \right] \end{equation}$$

@parameters k_p k_m

$$\begin{equation} \left[ \begin{array}{c} k_{p} \\ k_{m} \\ \end{array} \right] \end{equation}$$

@species A(t) B(t)

$$\begin{equation} \left[ \begin{array}{c} A\left( t \right) \\ B\left( t \right) \\ \end{array} \right] \end{equation}$$

A reaction network can be combined from several reactions:

r1p = Reaction(k_p, [A], [B], [1], [1])

k_p, A --> B

r1m = Reaction(k_m, [B], [A], [1], [1])

k_m, B --> A

rn1x = complete(ReactionSystem([r1p, r1m], t; name = :example1))

$$\begin{align*} \mathrm{A} &\xrightleftharpoons[k_{m}]{k_{p}} \mathrm{B} \end{align*}$$

Once we are here, the rest remains the same.

convert(ODESystem, rn1x)

prob1x = ODEProblem(rn1x, u1_ini, (0, 10.0), Dict(pairs(p1)))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
 1.0
 0.0

sol1x = solve(prob1x, Rosenbrock23())

	timestamp	A(t)	B(t)
1	0.0	1.0	0.0
2	0.000113386	0.999887	0.000113379
3	0.00124725	0.998754	0.0012464
4	0.00810266	0.991933	0.00806667
5	0.0183376	0.981846	0.0181539
6	0.0399115	0.960951	0.0390486
7	0.0695373	0.933054	0.066946
8	0.117184	0.890048	0.109952
9	0.177898	0.838412	0.161588
10	0.261426	0.772769	0.227231
...

doplots && plot(sol1x; size = (600, 200))

Example 2: A + 2B $\longrightleftharpoons$ AB_2

rn2 = @reaction_network rn2 begin
    @combinatoric_ratelaws false
    k_0A, ∅ --> A
    k_0B, ∅ --> B
    (k_1p, k_1m), A + 2B <--> AB_2
end

$$\begin{align*} \varnothing &\xrightarrow{k_{0A}} \mathrm{A} \\ \varnothing &\xrightarrow{k_{0B}} \mathrm{B} \\ \mathrm{A} + 2 \mathrm{B} &\xrightleftharpoons[k_{1m}]{k_{1p}} \mathrm{AB_{2}} \end{align*}$$

convert(ODESystem, rn2)

$$\begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= k_{0A} + k_{1m} \mathrm{AB}_{2}\left( t \right) - \left( B\left( t \right) \right)^{2} k_{1p} A\left( t \right) \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= k_{0B} + 2 k_{1m} \mathrm{AB}_{2}\left( t \right) - 2 \left( B\left( t \right) \right)^{2} k_{1p} A\left( t \right) \\ \frac{\mathrm{d} \mathrm{AB}_{2}\left( t \right)}{\mathrm{d}t} &= - k_{1m} \mathrm{AB}_{2}\left( t \right) + \left( B\left( t \right) \right)^{2} k_{1p} A\left( t \right) \end{align}$$

p2 = (k_0A = 0.5, k_0B = 1, k_1p = 0.1, k_1m = 1.0e-5)

(k_0A = 0.5, k_0B = 1, k_1p = 0.1, k_1m = 1.0e-5)

u2_ini = (A = 0, B = 0, AB_2 = 0)

(A = 0, B = 0, AB_2 = 0)

prob2 = ODEProblem(rn2, Dict(pairs(u2_ini)), (0, 20.0), Dict(pairs(p2)))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 20.0)
u0: 3-element Vector{Float64}:
 0.0
 0.0
 0.0

sol2 = solve(prob2, Rosenbrock23())

	timestamp	A(t)	B(t)	AB_2(t)
1	0.0	0.0	0.0	0.0
2	0.0001	5.0e-5	0.0001	6.25e-19
3	0.0011	0.00055	0.0011	1.08006e-14
4	0.0111	0.00555	0.0111	1.13501e-10
5	0.058436	0.0292179	0.0584358	9.95852e-8
6	0.0824519	0.0412254	0.0824509	5.19335e-7
7	0.141766	0.0708785	0.141757	4.69768e-6
8	0.178755	0.0893651	0.17873	1.23077e-5
9	0.241957	0.120937	0.241874	4.17015e-5
10	0.287619	0.143726	0.287451	8.40283e-5
...

doplots && plot(sol2; legend = :topleft, size = (600, 300))

Example 3: Catalysis for A + 2B $\rightleftharpoons$ AB_2

The same reaction as in example 2, but now with a catalyst C.

The reaction between A and B takes places when A and B are bound (adsorbed) to the catalyst. So we have adsorption reactions, reactions at the catalyst, and desorption. The overall of free and bound catalyst needs to be constant over time.

rn3 = @reaction_network rn3 begin
    @combinatoric_ratelaws false
    k_0A, ∅ --> A
    k_0B, ∅ --> B
    (k_1p, k_1m), A + C <--> CA
    (k_2p, k_2m), B + C <--> CB
    (k_3p, k_3m), CA + 2CB <--> CAB2 + 2C
    (k_4p, k_4m), CAB2 <--> AB2 + C
end

$$\begin{align*} \varnothing &\xrightarrow{k_{0A}} \mathrm{A} \\ \varnothing &\xrightarrow{k_{0B}} \mathrm{B} \\ \mathrm{A} + \mathrm{C} &\xrightleftharpoons[k_{1m}]{k_{1p}} \mathrm{CA} \\ \mathrm{B} + \mathrm{C} &\xrightleftharpoons[k_{2m}]{k_{2p}} \mathrm{CB} \\ \mathrm{CA} + 2 \mathrm{CB} &\xrightleftharpoons[k_{3m}]{k_{3p}} \mathrm{CAB2} + 2 \mathrm{C} \\ \mathrm{CAB2} &\xrightleftharpoons[k_{4m}]{k_{4p}} \mathrm{AB2} + \mathrm{C} \end{align*}$$

convert(ODESystem, rn3)

$$\begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= k_{0A} + k_{1m} \mathrm{CA}\left( t \right) - k_{1p} A\left( t \right) C\left( t \right) \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= k_{0B} + k_{2m} \mathrm{CB}\left( t \right) - k_{2p} B\left( t \right) C\left( t \right) \\ \frac{\mathrm{d} C\left( t \right)}{\mathrm{d}t} &= k_{1m} \mathrm{CA}\left( t \right) + k_{2m} \mathrm{CB}\left( t \right) + k_{4p} \mathrm{CAB2}\left( t \right) - k_{1p} A\left( t \right) C\left( t \right) - k_{2p} B\left( t \right) C\left( t \right) - k_{4m} C\left( t \right) \mathrm{AB2}\left( t \right) - 2 \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) + 2 \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \\ \frac{\mathrm{d} \mathrm{CA}\left( t \right)}{\mathrm{d}t} &= - k_{1m} \mathrm{CA}\left( t \right) + k_{1p} A\left( t \right) C\left( t \right) + \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) - \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \\ \frac{\mathrm{d} \mathrm{CB}\left( t \right)}{\mathrm{d}t} &= - k_{2m} \mathrm{CB}\left( t \right) + k_{2p} B\left( t \right) C\left( t \right) + 2 \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) - 2 \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \\ \frac{\mathrm{d} \mathrm{CAB2}\left( t \right)}{\mathrm{d}t} &= - k_{4p} \mathrm{CAB2}\left( t \right) + k_{4m} C\left( t \right) \mathrm{AB2}\left( t \right) - \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) + \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \\ \frac{\mathrm{d} \mathrm{AB2}\left( t \right)}{\mathrm{d}t} &= k_{4p} \mathrm{CAB2}\left( t \right) - k_{4m} C\left( t \right) \mathrm{AB2}\left( t \right) \end{align}$$

p3 = (k_0A = 0.5, k_0B = 1,
      k_1p = 10, k_1m = 0.1,
      k_2p = 10, k_2m = 0.1,
      k_3p = 10, k_3m = 0.1,
      k_4p = 10, k_4m = 0.1)

(k_0A = 0.5, k_0B = 1, k_1p = 10, k_1m = 0.1, k_2p = 10, k_2m = 0.1, k_3p = 10, k_3m = 0.1, k_4p = 10, k_4m = 0.1)

Cini=40

u3_ini = (A = 0, B = 0, CA = 0, CB = 0, CAB2 = 0, AB2 = 0, C = Cini)

(A = 0, B = 0, CA = 0, CB = 0, CAB2 = 0, AB2 = 0, C = 40)

t3end = 200

prob3 = ODEProblem(rn3, Dict(pairs(u3_ini)), (0, t3end), Dict(pairs(p3)))

ODEProblem with uType Vector{Float64} and tType Int64. In-place: true
timespan: (0, 200)
u0: 7-element Vector{Float64}:
  0.0
  0.0
 40.0
  0.0
  0.0
  0.0
  0.0

sol3 = solve(prob3, Rosenbrock23())

	timestamp	A(t)	B(t)	C(t)	CA(t)	CB(t)	CAB2(t)	AB2(t)
1	0.0	0.0	0.0	40.0	0.0	0.0	0.0	0.0
2	6.46784e-6	3.22974e-6	6.45949e-6	40.0	4.17882e-9	8.35764e-9	4.74654e-31	8.99173e-36
3	7.11463e-5	3.50726e-5	7.01452e-5	40.0	5.00556e-7	1.00111e-6	1.20711e-23	2.28831e-27
4	0.000175016	8.45196e-5	0.000169039	40.0	2.9886e-6	5.9772e-6	1.52686e-20	5.17114e-24
5	0.0003399	0.00015892	0.00031784	40.0	1.10301e-5	2.20603e-5	1.86355e-18	1.10544e-21
6	0.000559347	0.000250638	0.000501277	39.9999	2.9035e-5	5.807e-5	6.46331e-17	5.76826e-20
7	0.00085306	0.000361474	0.000722949	39.9998	6.50556e-5	0.000130111	1.1845e-15	1.55464e-18
8	0.00121011	0.000479825	0.00095965	39.9996	0.000125232	0.000250464	1.26499e-14	2.27535e-17
9	0.00165039	0.000604342	0.00120868	39.9993	0.000220853	0.000441706	9.78351e-14	2.37856e-16
10	0.00216779	0.000725254	0.00145051	39.9989	0.000358638	0.000717277	5.66062e-13	1.79861e-15
...

doplots && plot(sol3; legend = :topleft, size = (600, 300))

ctotal = rn3.C + rn3.CA + rn3.CB + rn3.CAB2

$$\begin{equation} \mathrm{CAB2}\left( t \right) + \mathrm{CA}\left( t \right) + \mathrm{CB}\left( t \right) + C\left( t \right) \end{equation}$$

sol3[ctotal]

148-element Vector{Float64}:
 40.0
 40.0
 40.0
 39.99999999999999
 39.99999999999999
 39.99999999999999
 39.99999999999999
  ⋮
 39.99999999999902
 39.99999999999902
 39.99999999999902
 39.99999999999902
 39.99999999999902
 39.99999999999926

@test sol3[ctotal]≈ fill(Cini,length(sol3))

Test Passed

Example 4: Heterogeneous catalysis

Heterogeneous catalysis assumes that the catalytic reaction takes place at surface. This means that reacting species need to be transported towards or away from the surface, and one has to model coupled transport and surface reaction.

Here we use VoronoiFVM.jl to model transport and Catalyst.jl to create the surface reaction network.

Problem specification

Assume $\Omega=(0,1)$ where a catalytic reaction takes place at $x=0$. We assume that the educts A, B, and the product AB2 are bulk species transported to the domain. At $x=1$ we set Dirichlet boundary conditions providing A,B and removing AB2.

A, B can adsorb at the catalyst at $x=0$ and react to AB2 while adsorbed. The product desorbs and is released to the bulk. So we have

Mass transport in the interior of $\Omega$:

$$\begin{aligned} \partial_t c_A + \nabla \cdot D_A \nabla c_A &=0\\ \partial_t c_B + \nabla \cdot D_B \nabla c_B &=0\\ \partial_t c_{AB2} + \nabla \cdot D_{AB2} \nabla c_{AB2} &=0 \end{aligned}$$

Coupled nonlinear robin boundary conditions at $x=0$:

$$\begin{aligned} D_A\partial_n c_A + r_1 &= 0\\ D_B\partial_n c_A + r_2 &= 0\\ D_{AB2}\partial_n c_{AB2} - r_4 &= 0\\ \end{aligned}$$

$r_1, r_2$ and $r_4$ are asorption/desorption reactions:

$$\begin{aligned} r_1&=k_{1p}c_A c_C - k_{1m}c_{CA}\\ r_2&=k_{2p}c_B c_C - k_{2m}c_{CB}\\ r_4&=k_{4p}c_{AB2} - k_{4m}c_{C_C}c_{C_{AB2}}\\ \end{aligned}$$

The free catalyst sites C and the catalyst coverages CA, CB, CAB2 behave according to:

$$\begin{equation} \frac{\mathrm{d} C\left( t \right)}{\mathrm{d}t} = k_{1m} \mathrm{CA}\left( t \right) + k_{2m} \mathrm{CB}\left( t \right) + k_{4p} \mathrm{CAB2}\left( t \right) - k_{1p} A\left( t \right) C\left( t \right) - k_{2p} B\left( t \right) C\left( t \right) - k_{4m} C\left( t \right) \mathrm{AB2}\left( t \right) - 2 \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) + 2 \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \end{equation}$$

$$\begin{equation} \frac{\mathrm{d} \mathrm{CA}\left( t \right)}{\mathrm{d}t} = - k_{1m} \mathrm{CA}\left( t \right) + k_{1p} A\left( t \right) C\left( t \right) + \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) - \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \end{equation}$$

$$\begin{equation} \frac{\mathrm{d} \mathrm{CB}\left( t \right)}{\mathrm{d}t} = - k_{2m} \mathrm{CB}\left( t \right) + k_{2p} B\left( t \right) C\left( t \right) + 2 \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) - 2 \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \end{equation}$$

$$\begin{equation} \frac{\mathrm{d} \mathrm{CAB2}\left( t \right)}{\mathrm{d}t} = - k_{4p} \mathrm{CAB2}\left( t \right) + k_{4m} C\left( t \right) \mathrm{AB2}\left( t \right) - \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) + \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \end{equation}$$

Dirichlet boundary conditions at $x=1$ :

$$\begin{aligned} c_A&=1\\ c_B&=1\\ c_{AB2}&=0 \end{aligned}$$

Finally, we set all initial concentrations to zero besides of the catalyst concenration (number of catalyst sites) $c_C|_{t=0}=C_0=1$.

Implementation

Surface reaction network

Define a reaction network under the assumption that the supply of A and B comes from the transport and does not need to be specified.

rnv = @reaction_network rnv begin
    @combinatoric_ratelaws false
    (k_1p, k_1m), A + C <--> CA
    (k_2p, k_2m), B + C <--> CB
    (k_3p, k_3m), CA + 2CB <--> CAB2 + 2C
    (k_4p, k_4m), CAB2 <--> AB2 + C
end

$$\begin{align*} \mathrm{A} + \mathrm{C} &\xrightleftharpoons[k_{1m}]{k_{1p}} \mathrm{CA} \\ \mathrm{B} + \mathrm{C} &\xrightleftharpoons[k_{2m}]{k_{2p}} \mathrm{CB} \\ \mathrm{CA} + 2 \mathrm{CB} &\xrightleftharpoons[k_{3m}]{k_{3p}} \mathrm{CAB2} + 2 \mathrm{C} \\ \mathrm{CAB2} &\xrightleftharpoons[k_{4m}]{k_{4p}} \mathrm{AB2} + \mathrm{C} \end{align*}$$

odesys=convert(ODESystem, rnv)

$$\begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= k_{1m} \mathrm{CA}\left( t \right) - k_{1p} A\left( t \right) C\left( t \right) \\ \frac{\mathrm{d} C\left( t \right)}{\mathrm{d}t} &= k_{1m} \mathrm{CA}\left( t \right) + k_{2m} \mathrm{CB}\left( t \right) + k_{4p} \mathrm{CAB2}\left( t \right) - k_{1p} A\left( t \right) C\left( t \right) - k_{2p} B\left( t \right) C\left( t \right) - k_{4m} C\left( t \right) \mathrm{AB2}\left( t \right) - 2 \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) + 2 \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \\ \frac{\mathrm{d} \mathrm{CA}\left( t \right)}{\mathrm{d}t} &= - k_{1m} \mathrm{CA}\left( t \right) + k_{1p} A\left( t \right) C\left( t \right) + \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) - \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= k_{2m} \mathrm{CB}\left( t \right) - k_{2p} B\left( t \right) C\left( t \right) \\ \frac{\mathrm{d} \mathrm{CB}\left( t \right)}{\mathrm{d}t} &= - k_{2m} \mathrm{CB}\left( t \right) + k_{2p} B\left( t \right) C\left( t \right) + 2 \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) - 2 \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \\ \frac{\mathrm{d} \mathrm{CAB2}\left( t \right)}{\mathrm{d}t} &= - k_{4p} \mathrm{CAB2}\left( t \right) + k_{4m} C\left( t \right) \mathrm{AB2}\left( t \right) - \left( C\left( t \right) \right)^{2} k_{3m} \mathrm{CAB2}\left( t \right) + \left( \mathrm{CB}\left( t \right) \right)^{2} k_{3p} \mathrm{CA}\left( t \right) \\ \frac{\mathrm{d} \mathrm{AB2}\left( t \right)}{\mathrm{d}t} &= k_{4p} \mathrm{CAB2}\left( t \right) - k_{4m} C\left( t \right) \mathrm{AB2}\left( t \right) \end{align}$$

For coupling with VoronoiFVM we need species numbers which need to correspond to the species in our network:

begin
    smap = speciesmap(rnv)
    const iA = smap[rnv.A]
    const iB = smap[rnv.B]
    const iC = smap[rnv.C]
    const iCA = smap[rnv.CA]
    const iCB = smap[rnv.CB]
    const iCAB2 = smap[rnv.CAB2]
    const iAB2 = smap[rnv.AB2]
end;

Grid:

grid = simplexgrid(0:0.01:1)

ExtendableGrids.ExtendableGrid{Float64, Int32}
      dim =       1
   nnodes =     101
   ncells =     100
  nbfaces =       2

gridplot(grid, size=(600,100))

The grid has two boundary regions: region 1 at x=0 and region 2 at x=1.

Reaction parameters:

pcat = (k_1p = 50, k_1m = 0.1,
        k_2p = 50, k_2m = 0.1,
        k_3p = 10, k_3m = 0.1,
        k_4p = 50, k_4m = 0.1)

(k_1p = 50, k_1m = 0.1, k_2p = 50, k_2m = 0.1, k_3p = 10, k_3m = 0.1, k_4p = 50, k_4m = 0.1)

Parameters for the VoronoiFVM system:

params = (D_A = 1.0,
          D_B = 1.0,
          D_AB2 = 1.0,
          pcat = pcat)

(D_A = 1.0, D_B = 1.0, D_AB2 = 1.0, pcat = (k_1p = 50, k_1m = 0.1, k_2p = 50, k_2m = 0.1, k_3p = 10, k_3m = 0.1, k_4p = 50, k_4m = 0.1))

Initial values for the reaction network (needed only for the definition of the ODE problem)

C0 = 1.0

1.0

uv_ini = (A = 0, B = 0, CA = 0, CB = 0, CAB2 = 0, AB2 = 0, C = C0)

(A = 0, B = 0, CA = 0, CB = 0, CAB2 = 0, AB2 = 0, C = 1.0)

tvend = 200.0

200.0

const probv = ODEProblem(rnv, Dict(pairs(uv_ini)), (0, tvend), Dict(pairs(pcat)))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 200.0)
u0: 7-element Vector{Float64}:
 0.0
 1.0
 0.0
 0.0
 0.0
 0.0
 0.0

Callback functions for VoronoiFVM

First, define flux and storage functions for the bulk process:

function storage(y, u, node, p)
    y[iA] = u[iA]
    y[iB] = u[iB]
    y[iAB2] = u[iAB2]
end

storage (generic function with 1 method)

function flux(y, u, edge, p)
    (; D_A, D_B, D_AB2) = p
    y[iA] = D_A * (u[iA, 1] - u[iA, 2])
    y[iB] = D_B * (u[iB, 1] - u[iB, 2])
    y[iAB2] = D_A * (u[iAB2, 1] - u[iAB2, 2])
end

flux (generic function with 1 method)

Storage term for the surface reaction:

function bstorage(y, u, bnode, p)
    y[iC] = u[iC]
    y[iCA] = u[iCA]
    y[iCB] = u[iCB]
    y[iCAB2] = u[iCAB2]
end

bstorage (generic function with 1 method)

Catalytic reaction. Here we use the right hand side function of the ODE problem generated above. In VoronoiFVM, reaction term are a the left hand side, so we need to multiply by -1.

Note that we need to pass the parameter record as generated for the ODE problem instead of pcat.

function catreaction(f, u, bnode, p)
    probv.f(f, u, probv.p, bnode.time)
    for i = 1:length(f)
        f[i] = -f[i]
    end
end

catreaction (generic function with 1 method)

Define the Dirichlet boundary condition at x=1 (region 2):

function bulkbc(f, u, bnode, p)
    v = ramp(bnode.time; du = (0.0, 1.0), dt = (0.0, 0.01))
    boundary_dirichlet!(f, u, bnode; species = iA, value = v, region = 2)
    boundary_dirichlet!(f, u, bnode; species = iB, value = v, region = 2)
    boundary_dirichlet!(f, u, bnode; species = iAB2, value = 0, region = 2)
end

bulkbc (generic function with 1 method)

Dispatch the boundary conditions

function breaction(f, u, bnode, p)
    if bnode.region == 1
        catreaction(f,u,bnode,p)
    else
 		bulkbc(f,u,bnode,p)
    end
end

breaction (generic function with 1 method)

Coupled transport-reaction system

Define a VoronoiFVM system from grid, params and the callback functions and enable the bulk and boundary species. unknown_storage = :sparse means that the solution is stored as a nspecies x nnodes sparse matrix in order to take into account that the surface species are non-existent in the bulk. unknown_storage = :dense would store a full matrix and solve dummy equations for the surface species values in the bulk.

begin
    sys = VoronoiFVM.System(grid; data = params,
                            flux, breaction, bstorage, storage,
                            unknown_storage = :sparse)
    enable_species!(sys, iA, [1])
    enable_species!(sys, iB, [1])
    enable_species!(sys, iAB2, [1])

    enable_boundary_species!(sys, iC, [1])
    enable_boundary_species!(sys, iCA, [1])
    enable_boundary_species!(sys, iCB, [1])
    enable_boundary_species!(sys, iCAB2, [1])
end;

Define an initial value for sys:

begin
    u0 = VoronoiFVM.unknowns(sys; inival = 0)
    u0[iC, 1] = C0
end;

Solution

Solve the time evolution

tsol = solve(sys; inival = u0, times = (1.0e-4, tvend));

let
    t_plot = round(10^log_t_plot; sigdigits = 3)
    vis = GridVisualizer(;size = (600, 300), flimits = (0, 1), title = "Bulk concentrations: t=$t_plot", legend = :lt)
    sol = tsol(t_plot)
    scalarplot!(vis, grid, sol[iA, :]; color = :red, label = "A")
    scalarplot!(vis, grid, sol[iB, :]; color = :green, label = "B", clear = false)
    scalarplot!(vis, grid, sol[iAB2, :]; color = :blue, label = "AB2", clear = false)
    reveal(vis)
end

Ctotalv = tsol[iC, 1, :] + tsol[iCA, 1, :] + tsol[iCB, 1, :] + tsol[iCAB2, 1, :]

248-element Vector{Float64}:
 1.0
 1.0
 0.9999999999999999
 1.0
 1.0
 0.9999999999999999
 0.9999999999999999
 ⋮
 1.0000000000000002
 1.0000000000000002
 1.0000000000000002
 1.0000000000000002
 1.0000000000000002
 1.0000000000000002

@test Ctotalv≈ ones(length(tsol))

Test Passed

let
    vis = GridVisualizer(; size = (600, 300),
                          xlabel="t",
                         flimits = (0, 1), xlimits = (1.0e-3, tvend),
                         legend = :lt, title = "Concentrations at x=0", xscale = :log10)
    t = tsol.t
    scalarplot!(vis, t, tsol[iA, 1, :]; color = :darkred, label = "A")
    scalarplot!(vis, t, tsol[iCA, 1, :]; color = :red, label = "CA")
    scalarplot!(vis, t, tsol[iB, 1, :]; color = :darkgreen, label = "B")
    scalarplot!(vis, t, tsol[iCB, 1, :]; color = :green, label = "CB")
    scalarplot!(vis, t, tsol[iAB2, 1, :]; color = :darkblue, label = "AB2")
    scalarplot!(vis, t, tsol[iCAB2, 1, :]; color = :blue, label = "CAB2")
    scalarplot!(vis, t, tsol[iC, 1, :] / C0; color = :orange, label = "C/C0")
    scalarplot!(vis, t, Ctotalv / C0; color = :darkorange, label = "Ctot/C0")

    reveal(vis)
end

Appendix

# Some tricks helping to run the notebook during VoronoiFVM.jl CI
begin
    doplots=!haskey(ENV,"VORONOIFVM_RUNTESTS")
    import Pkg as _Pkg
    haskey(ENV, "PLUTO_PROJECT") && _Pkg.activate(ENV["PLUTO_PROJECT"])
    using Revise
end

Built with Julia 1.11.1 and

Towards Heterogeneous Catalysis

Mass action kinetics

Example 1: A \(\longrightleftharpoons\) B

Solution via plain ODE problem using OrdinaryDiffEq.jl:

`Catalyst.@reaction_network`

Unraveling `@reaction_network`

Example 2: A + 2B \(\longrightleftharpoons\) AB_2

Example 3: Catalysis for A + 2B \(\rightleftharpoons\) AB_2

Example 4: Heterogeneous catalysis

Problem specification

Implementation

Surface reaction network

Callback functions for VoronoiFVM

Coupled transport-reaction system

Solution

Appendix

Towards Heterogeneous Catalysis

Mass action kinetics

Example 1: A \(\longrightleftharpoons\) B

Solution via plain ODE problem using OrdinaryDiffEq.jl:

Catalyst.@reaction_network

Unraveling @reaction_network

Example 2: A + 2B \(\longrightleftharpoons\) AB_2

Example 3: Catalysis for A + 2B \(\rightleftharpoons\) AB_2

Example 4: Heterogeneous catalysis

Problem specification

Implementation

Surface reaction network

Callback functions for VoronoiFVM

Coupled transport-reaction system

Solution

Appendix

`Catalyst.@reaction_network`

Unraveling `@reaction_network`