GitHub

Composing Models

Illustrative Example

Here is the complete example of composing, visualizing and solving the SuperFast model and the Fast-JX model, with explanation to follow:

using EarthSciMLBase, GasChem, ModelingToolkit, Dates, DynamicQuantities, DifferentialEquations
using ModelingToolkit:t

composed_ode = couple(SuperFast(), FastJX()) # Compose two models use the "couple" function

start = Dates.datetime2unix(Dates.DateTime(2024, 2, 29))
tspan = (start, start+3600*24*3)
sys = convert(ODESystem, composed_ode) # Define the coupled system
sol = solve(ODEProblem(sys, [], tspan, []), AutoTsit5(Rosenbrock23()), saveat=10.0) # Solve the coupled system
retcode: Success
Interpolation: 1st order linear
t: 25921-element Vector{Float64}:
 1.7091648e9
 1.70916481e9
 1.70916482e9
 1.70916483e9
 1.70916484e9
 1.70916485e9
 1.70916486e9
 1.70916487e9
 1.70916488e9
 1.70916489e9
 ⋮
 1.70942392e9
 1.70942393e9
 1.70942394e9
 1.70942395e9
 1.70942396e9
 1.70942397e9
 1.70942398e9
 1.70942399e9
 1.709424e9
u: 25921-element Vector{Vector{Float64}}:
 [20.0, 0.01, 0.01, 10.0, 10.0, 0.01, 0.15, 275.0, 1.6, 0.15, 2.34, 10.0]
 [19.305187270750647, 4.010708409457426e-6, 1.2044216454383465e-5, 9.236532283451524, 10.735785003994067, 4.667140999774169e-6, 0.17166113009807588, 274.9848955386636, 1.5971517961737738, 0.14628099763685165, 2.3390415080297546, 10.027682710104324]
 [18.681194115885813, 3.1012039715983065e-7, 1.103861594571725e-6, 8.612465874701202, 11.35982268892144, 8.287598327739945e-7, 0.17169312669149317, 274.9848813497871, 1.5971490457124489, 0.14627122289139596, 2.3390405740081532, 10.027711431320755]
 [18.11729674476813, 2.9285712802874957e-7, 1.1372937755482663e-6, 8.048538437222605, 11.923739164797766, 8.548866943074271e-7, 0.17171359346159062, 274.9848764721262, 1.5971480614068234, 0.1462639081490897, 2.3390402397690697, 10.027722390177278]
 [17.606177186619, 2.7180342303463286e-7, 1.152442476084169e-6, 7.5373901989043, 12.434876712123145, 8.757005328447629e-7, 0.17173330603909756, 274.984871943218, 1.5971471449846704, 0.14625685834682495, 2.339039928582638, 10.02773307829717]
 [17.139134594178415, 2.550935062751034e-7, 1.175345388227008e-6, 7.07032014262137, 12.901936358827822, 9.002422020711651e-7, 0.1717523588421338, 274.9848677031774, 1.5971462849795315, 0.14625003790350244, 2.339039636554235, 10.027743484889994]
 [16.713927054496416, 2.412131315669287e-7, 1.2038522514314455e-6, 6.64508628208911, 13.327160077165612, 9.273166518509935e-7, 0.17177078459872241, 274.9848637323308, 1.5971454774728155, 0.1462434359349939, 2.3390393623540104, 10.027753623990897]
 [16.322474318082474, 2.289320437871502e-7, 1.2340422682722444e-6, 6.253608138451942, 13.718628306949988, 9.557415283492466e-7, 0.1717886913652446, 274.984859967623, 1.5971447103980918, 0.14623701573084918, 2.339039101883418, 10.02776353465713]
 [15.963325836923904, 2.185418717604736e-7, 1.2663847162493541e-6, 5.894435126375958, 14.077791629538673, 9.859757538093207e-7, 0.17180609703082841, 274.98485639773, 1.5971439816589812, 0.1462307710472893, 2.339038854430057, 10.027773220867378]
 [15.6328039366119, 2.0897055137256102e-7, 1.299006943687936e-6, 5.563889460546513, 14.40832780393829, 1.0157624717728983e-6, 0.1718230578962959, 274.9848529969704, 1.597143285923486, 0.14622468286655338, 2.3390386181856195, 10.027782708938556]
 ⋮
 [11.368865267466012, 3.6849243850129676e-8, 3.741892420561553e-6, 0.36749983721158275, 19.07230899121884, 3.0619809040671543e-6, 0.1104472433456204, 275.1514242808599, 1.5338305173069013, 0.04638333195677815, 2.2873998647711637, 10.560086614372382]
 [11.363338240776624, 3.679254212351943e-8, 3.79462058158918e-6, 0.36196850832346045, 19.07783814403022, 3.107056804662112e-6, 0.11045069885179601, 275.1514237209915, 1.5338303992831277, 0.04638209791713522, 2.28739982397639, 10.560088785950086]
 [11.357619217199636, 3.671588885038599e-8, 3.84827533841261e-6, 0.3562451916756907, 19.08355928919722, 3.153310783712789e-6, 0.11045415173926479, 275.15142316125224, 1.5338302819405025, 0.04638086463679477, 2.287399783427608, 10.560090952930448]
 [11.351708196735046, 3.661928403072934e-8, 3.902856691031845e-6, 0.35032988726827347, 19.08947242671983, 3.200742841219185e-6, 0.11045760200802672, 275.15142260164225, 1.5338301652790258, 0.0463796321157568, 2.2873997431248174, 10.560093115313464]
 [11.345605179382858, 3.650272766454948e-8, 3.958364639446884e-6, 0.34422259510120873, 19.095577556598055, 3.2493529771813004e-6, 0.11046104965808182, 275.1514220421614, 1.533830049298698, 0.046378400354021306, 2.287399703068018, 10.560095273099137]
 [11.339310165143068, 3.6366219751846425e-8, 4.0147991836577255e-6, 0.3379233151744966, 19.101874678831898, 3.2991411915991347e-6, 0.11046449468943008, 275.1514214828098, 1.5338299339995187, 0.0463771693515883, 2.2873996632572107, 10.560097426287465]
 [11.33282315401568, 3.620976029262016e-8, 4.072160323664372e-6, 0.3314320474881369, 19.108363793421354, 3.3501074844726884e-6, 0.1104679371020715, 275.1514209235873, 1.5338298193814879, 0.046375939108457775, 2.287399623692395, 10.560099574878448]
 [11.326144146000692, 3.603334928687069e-8, 4.130448059466822e-6, 0.32474879204212986, 19.115044900366424, 3.402251855801961e-6, 0.1104713768960061, 275.1514203644941, 1.5338297054446057, 0.04637470962462972, 2.2873995843735706, 10.560101718872088]
 [11.319576432745334, 3.5426494042188535e-8, 4.190075640343612e-6, 0.3181768501115421, 19.121614701422036, 3.4596003893355346e-6, 0.11047480896653658, 275.15141980739264, 1.5338295920622287, 0.04637348096546815, 2.2873995452480416, 10.56010385523839]

In the composed system, the variable name for O₃ is not O3 but superfast₊O3(t). So we need some preparation of the result before visualizing. 

vars = unknowns(sys)  # Get the variables in the composed system
var_dict = Dict(string(var) => var for var in vars)
pols = ["O3", "OH", "NO", "NO2", "CH3O2", "CO","CH3OOH", "CH2O", "ISOP"]
var_names_p = ["SuperFast₊$(v)(t)" for v in pols]

x_t = unix2datetime.(sol[t]) # Convert from unixtime to date time for visualizing
25921-element Vector{Dates.DateTime}:
 2024-02-29T00:00:00
 2024-02-29T00:00:10
 2024-02-29T00:00:20
 2024-02-29T00:00:30
 2024-02-29T00:00:40
 2024-02-29T00:00:50
 2024-02-29T00:01:00
 2024-02-29T00:01:10
 2024-02-29T00:01:20
 2024-02-29T00:01:30
 ⋮
 2024-03-02T23:58:40
 2024-03-02T23:58:50
 2024-03-02T23:59:00
 2024-03-02T23:59:10
 2024-03-02T23:59:20
 2024-03-02T23:59:30
 2024-03-02T23:59:40
 2024-03-02T23:59:50
 2024-03-03T00:00:00

Then, we could plot the results as:

using Plots
pp = []
for (i, v) in enumerate(var_names_p)
    name = pols[i]
    push!(pp, Plots.plot(x_t,sol[var_dict[v]],label = "$name", size = (1000, 600), xrotation=45))
end
Plots.plot(pp..., layout=(3, 3))
Example block output

Adding Emission Data

GasChem.jl incorporates an emissions model that utilizes data from the US National Emissions Inventory for the year 2016. This model is activated as an extension when the EarthSciData package is used. Here's a simple example:

using GasChem, EarthSciData # This will trigger the emission extension
using Dates, ModelingToolkit, DifferentialEquations, EarthSciMLBase
using Plots, DynamicQuantities
using ModelingToolkit:t

domain = DomainInfo(
    Dates.DateTime(2016, 5, 1), Dates.DateTime(2016, 5, 4);
    lonrange = deg2rad(-129):deg2rad(1):deg2rad(-61),
    latrange = deg2rad(11):deg2rad(1):deg2rad(59),
    levrange = 1:1:3,
    dtype = Float64)

emis = NEI2016MonthlyEmis("mrggrid_withbeis_withrwc", domain)

model_noemis = couple(SuperFast(),FastJX()) # A model with chemistry and photolysis, but no emissions.
model_withemis = couple(SuperFast(), FastJX(), emis) # The same model with emissions.

sys_noemis = convert(ODESystem, model_noemis)
sys_withemis = convert(ODESystem, model_withemis)

tspan = EarthSciMLBase.get_tspan(domain)
sol_noemis = solve(ODEProblem(sys_noemis, [], tspan, []), AutoTsit5(Rosenbrock23()))
sol_withemis = solve(ODEProblem(sys_withemis, [], tspan, []), AutoTsit5(Rosenbrock23()))

vars = unknowns(sys_noemis)  # Get the variables in the composed system
var_dict = Dict(string(var) => var for var in vars)
pols = ["O3", "OH", "NO", "NO2", "CH3O2", "CO","CH3OOH", "CH2O", "ISOP"]
var_names_p = ["SuperFast₊$(v)(t)" for v in pols]

pp = []
for (i, v) in enumerate(var_names_p)
    p = plot(unix2datetime.(sol_noemis[t]),sol_noemis[var_dict[v]], label="No Emissions", title = v,
        size = (1000, 600), xrotation=45)
    plot!(unix2datetime.(sol_withemis[t]),sol_withemis[var_dict[v]], label="With Emissions", )
    push!(pp, p)
end
Plots.plot(pp..., layout=(3, 3))
Example block output