Coupling reaction to transport solver

[SA8416]

As a follow-up to:

https://discourse.julialang.org/t/error-when-calculating-rhs-of-reaction-ode-no-method-matching-float64-num/64420

I have a set of 22 aqueous species and 44 reactions I am modelling. The species are transported through a 2D porous medium (by diffusion, advection and electromigration), which is modelled using https://github.com/simulkade/JFVM.j

I have opted for a SNIA to handle the transport and reaction stages seperately, but am having issues using Catalyst.jl to use the transported species concentrations.

One such issue is the handling of the water species in my reactions. e.g. for one of my reactions:

(kf,kb), H₂O ⇌ H⁺ + OH⁻

which I found a workaround for by representing the equation as:

(kf*aw,kb), ∅ ⇌ H⁺ + OH⁻

instead, where the activity of water (aw) is moved from being a reactant to a multiplier of the forward rate constant (as the supply of water is effectively infinite) .

Another issue (that I suspect stems from my own inability) is inputting the transported concentrations to Catalyst.jl at each time-step.

Here is some psedocode that summarizes the main part of my problem:

m = createMeshCylindrical2D(Nr, Nz, Lr, Lz)
for t in dt:t_end
    # ∇²T = 0 , uses T_saturation as Boundary Condition
    Temperature = JFVM.solvePDE(m,LHS_T,RHS_T) 
    # ∇²P = 0 , uses the temperature gradient as a Boundary Condition
    Pressure = JFVM.solvePDE(m,LHS_P,RHS_P)
    u = f1(Pressure) # water velocity
    # ∂C/∂t = -D∇²C + u⋅∇C , uses the water velocity for the advection of solutes
    c_transport = [JFVM.solvePDE(m,LHS_c[i],RHS_c[i]) for i in all_species]
    # c_reaction to be calculated here...
    T_saturation = f2(c_reaction)
end

My code is quite bloated and illegible atm, so I am working on cleaning it up and putting together a MWE that I'll post here.

[ChrisRackauckas]

Maybe this could be a good model to use for testing the upcoming spatial modeling functionality @isaacsas?

[isaacsas]

@SA8416 what you did to keep the H20 amount constant is fine. You can also have a species within a reaction that stays constant if you do kf, H20 --> H⁺ + OH⁻ + H20. We don't support "constant" species, so you need to reduce the reaction order and move the species into the rate function (as you did), or you need to balance the reaction so the amount of the species is unchanged.

It sounds like you want to just be able to evaluate the reaction terms from the corresponding ODE model at each spatial location, is this correct? If so, suppose that u is a matrix with u[i,j] the amount of species i at location j. Then you should do something like:

rn = @reaction_network begin
            ...
        end
osys = convert(ODESystem, rn)
of = ODEFunction(osys)
f2! = of.f.f_iip

# to evaluate u in-place 
u = [current values for each species at each location]
du = copy(u)  # to store the ODE terms for each species at each location
p  = [parameter vector]
t   = 0.0 
for i = 1:num_spatial_locations
   ducol = @view du[:,i]
   ucol   = @view u[:,i]
    f2!(ducol, ucol, p, t)
end

[SA8416]

I think I've got everything formatted correctly, as per your example, but I'm getting an error:

ERROR: LoadError: UndefVarError: f2! not defined

[isaacsas]

There was a typo in my code it should be f2! = of.f.f_iip.

[yewalenikhil65]

@SA8416 what you did to keep the H20 amount constant is fine. You can also have a species within a reaction that stays constant if you do kf, H20 --> H⁺ + OH⁻ + H20. We don't support "constant" species, so you need to reduce the reaction order and move the species into the rate function (as you did), or you need to balance the reaction so the amount of the species is unchanged.

It sounds like you want to just be able to evaluate the reaction terms from the corresponding ODE model at each spatial location, is this correct? If so, suppose that u is a matrix with u[i,j] the amount of species i at location j. Then you should do something like:

rn = @reaction_network begin
            ...
        end
osys = convert(ODESystem, rn)
of = ODEFunction(osys)
f2! = of.f.f_iip

# to evaluate u in-place 
u = [current values for each species at each location]
du = copy(u)  # to store the ODE terms for each species at each location
p  = [parameter vector]
t   = 0.0 
for i = 1:num_spatial_locations
   ducol = @view du[:,i]
   ucol   = @view u[:,i]
    f2!(ducol, ucol, p, t)
end

I would like to contribute to this (Advection Reaction Diffusion PDEs) support to Catalyst. Any ideas on which packages would be good to explore(for solving PDEs) in case we are interested in Mass transport problems with catalyst?

[isaacsas]

@yewalenikhil65 this is a big project as it will involve updates to ModelingToolkit's PDE support (both to handle these types of PDEs and to handle domain geometries), design changes here to support spatial systems and specifying geometries/compartments, and more. It is actually something we are hoping to work on in the near future, but I'm not sure we're quite there yet to invest time on it at the Catalyst level.

If you are looking for more to contribute would you have any interest in filling out some network analysis-type functionality? Now that we've got these various matrix representations you added, are there some nice things one can use them to conclude about reaction networks a priori? (Calculating the deficiency index, analyzing steady-state behavior from the network, etc.)

[isaacsas]

If you really want to work towards the PDE stuff, the first step is probably learning about the ModelingToolkit PDESystems and how they work / dispatch to DiffEqOperators and such.

[yewalenikhil65]

If you are looking for more to contribute would you have any interest in filling out some network analysis-type functionality? Now that we've got these various matrix representations you added, are there some nice things one can use them to conclude about reaction networks a priori? (Calculating the deficiency index, analyzing steady-state behavior from the network, etc.)

Certainly, I would contribute to this soon. But i think there are few more functionalities to add before coming with tutorial to conclude