Example using different components of EarthSciMLBase together

This example shows how to define and couple different components of EarthSciMLBase together to create a more complex model. First, we create several model components.

Our first example system is a simple reaction system:

using EarthSciMLBase
using ModelingToolkit, Catalyst, DomainSets, DifferentialEquations
using ModelingToolkit: t_nounits, D_nounits
t = t_nounits
D = D_nounits
using Plots

# Create our independent variable `t` and our partially-independent variables `x` and `y`.
@parameters x y

struct ExampleSys1Coupler
    sys
end
function ExampleSys1()
    @species c₁(t)=5.0 c₂(t)=5.0
    rs = ReactionSystem(
        [Reaction(2.0, [c₁], [c₂])],
        t; name = :Sys1, combinatoric_ratelaws = false)
    ode_model(complete(rs); metadata = Dict(CoupleType => ExampleSys1Coupler))
end

ExampleSys1()

\[ \begin{align} \frac{\mathrm{d} ~ \mathtt{c_1}\left( t \right)}{\mathrm{d}t} &= - 2 ~ \mathtt{c_1}\left( t \right) \\ \frac{\mathrm{d} ~ \mathtt{c_2}\left( t \right)}{\mathrm{d}t} &= 2 ~ \mathtt{c_1}\left( t \right) \end{align} \]

Our second example system is a simple ODE system, with the same two variables.

struct ExampleSys2Coupler
    sys
end
function ExampleSys2()
    @variables c₁(t)=5.0 c₂(t)=5.0
    @parameters p₁=1.0 p₂=0.5 x=1 y=1
    System(
        [D(c₁) ~ -p₁, D(c₂) ~ p₂],
        t, [c₁, c₂], [p₁, p₂, x, y]; name = :Sys2,
        metadata = Dict(CoupleType => ExampleSys2Coupler))
end

ExampleSys2()

\[ \begin{align} \frac{\mathrm{d} ~ \mathtt{c_1}\left( t \right)}{\mathrm{d}t} &= - \mathtt{p{_1}} \\ \frac{\mathrm{d} ~ \mathtt{c_2}\left( t \right)}{\mathrm{d}t} &= \mathtt{p{_2}} \end{align} \]

Now, we specify what should happen when we couple the two systems together. In this case, we want the the derivative of the composed system to be equal to the sum of the derivatives of the two systems. We can do that using the operator_compose function from this package.

function EarthSciMLBase.couple2(sys1::ExampleSys1Coupler, sys2::ExampleSys2Coupler)
    sys1, sys2 = sys1.sys, sys2.sys
    sys1 = convert(System, sys1)
    operator_compose(sys1, sys2)
end

Once we specify all of the above, it is simple to create our two individual systems and then couple them together.

sys1 = ExampleSys1()
sys2 = ExampleSys2()
sys = couple(sys1, sys2)

convert(System, sys)

\[ \begin{align} \frac{\mathrm{d} ~ \mathtt{Sys1.c_2}\left( t \right)}{\mathrm{d}t} &= \mathtt{Sys1.Sys2\_ddt\_c_{2}ˍt}\left( t \right) + 2 ~ \mathtt{Sys1.c_1}\left( t \right) \\ \frac{\mathrm{d} ~ \mathtt{Sys1.c_1}\left( t \right)}{\mathrm{d}t} &= \mathtt{Sys1.Sys2\_ddt\_c_{1}ˍt}\left( t \right) - 2 ~ \mathtt{Sys1.c_1}\left( t \right) \end{align} \]

At this point we have an ODE system that is composed of two other ODE systems. We can inspect its observed variables using the observed function.

simplified_sys = convert(System, sys)

\[ \begin{align} \frac{\mathrm{d} ~ \mathtt{Sys1.c_2}\left( t \right)}{\mathrm{d}t} &= \mathtt{Sys1.Sys2\_ddt\_c_{2}ˍt}\left( t \right) + 2 ~ \mathtt{Sys1.c_1}\left( t \right) \\ \frac{\mathrm{d} ~ \mathtt{Sys1.c_1}\left( t \right)}{\mathrm{d}t} &= \mathtt{Sys1.Sys2\_ddt\_c_{1}ˍt}\left( t \right) - 2 ~ \mathtt{Sys1.c_1}\left( t \right) \end{align} \]

observed(simplified_sys)

\[ \begin{align} \mathtt{Sys2.Sys2\_ddt\_c_{2}ˍt}\left( t \right) &= \mathtt{Sys2.p{_2}} \\ \mathtt{Sys2.Sys2\_ddt\_c_{1}ˍt}\left( t \right) &= - \mathtt{Sys2.p{_1}} \\ \mathtt{Sys2.c_2}\left( t \right) &= \mathtt{Sys1.c_2}\left( t \right) \\ \mathtt{Sys2.c_1}\left( t \right) &= \mathtt{Sys1.c_1}\left( t \right) \\ \mathtt{Sys1.Sys2\_ddt\_c_{2}ˍt}\left( t \right) &= \mathtt{Sys2.Sys2\_ddt\_c_{2}ˍt}\left( t \right) \\ \mathtt{Sys1.Sys2\_ddt\_c_{1}ˍt}\left( t \right) &= \mathtt{Sys2.Sys2\_ddt\_c_{1}ˍt}\left( t \right) \end{align} \]

We can also run simulations using this system:

odeprob = ODEProblem(simplified_sys, [], (0.0, 10.0))
odesol = solve(odeprob)
plot(odesol)