Update documentation

docs/src/composing_models.md

@@ -12,8 +12,8 @@ model and the Fast-JX model, with explanation to follow:
 using EarthSciMLBase, GasChem, ModelingToolkit, Dates, Unitful, DifferentialEquations
 
 
-@parameters t
-composed_ode = SuperFast(t) + FastJX(t) # Compose two models simply use the "+" operator
+@parameters t [unit = u"s", description = "Time"]
+composed_ode = couple(SuperFast(t), FastJX(t)) # Compose two models use the "couple" function
 
 start = Dates.datetime2unix(Dates.DateTime(2024, 2, 29))
 tspan = (start, start+3600*24*3)
@@ -35,9 +35,8 @@ Then, we could plot the results as:
 ```@setup 1
 using EarthSciMLBase, GasChem, ModelingToolkit, Dates, Unitful, DifferentialEquations
 
-
-@parameters t
-composed_ode = SuperFast(t) + FastJX(t) # Compose two models simply use the "+" operator
+@parameters t [unit = u"s"]
+composed_ode = couple(SuperFast(t), FastJX(t)) # Compose two models use the "couple" function
 
 start = Dates.datetime2unix(Dates.DateTime(2024, 2, 29))
 tspan = (start, start+3600*24*3)
@@ -71,12 +70,16 @@ using GasChem, EarthSciData # This will trigger the emission extension
 using Dates, ModelingToolkit, DifferentialEquations, EarthSciMLBase
 using Plots, Unitful
 
-ModelingToolkit.check_units(eqs...) = nothing 
-@parameters t
-@constants Δz = 60.0 [unit=u"m"]
-model_noemis = SuperFast(t)+FastJX(t) # A model with chemistry and photolysis, but no emissions.
-model_withemis = SuperFast(t)+FastJX(t)+ 
-    NEI2016MonthlyEmis{Float64}("mrggrid_withbeis_withrwc", t, -100.0, 30.0, 1.0, Δz) # The same model with emissions.
+@parameters t [unit = u"s", description = "Time"]
+@parameters lat = 40 
+@parameters lon = -97 
+@parameters lev = 1
+@parameters Δz = 60 [unit = u"m"]
+emis = NEI2016MonthlyEmis("mrggrid_withbeis_withrwc", t, lon, lat, lev, Δz; dtype=Float64) 
+
+
+model_noemis = couple(SuperFast(t),FastJX(t)) # A model with chemistry and photolysis, but no emissions.
+model_withemis = couple(SuperFast(t), FastJX(t), emis) # The same model with emissions.
 
 sys_noemis = structural_simplify(get_mtk(model_noemis))
 sys_withemis = structural_simplify(get_mtk(model_withemis))
@@ -86,16 +89,13 @@ tspan = (start, start+3*24*3600)
 sol_noemis = solve(ODEProblem(sys_noemis, [], tspan, []), AutoTsit5(Rosenbrock23()))
 sol_withemis = solve(ODEProblem(sys_withemis, [], tspan, []), AutoTsit5(Rosenbrock23()))
 
-plot(
-    plot(sol_noemis, title="Model without emissions"),
-    plot!(sol_withemis, title="Model with emissions"),
-)
-```
+vars = states(sys_noemis)  # Get the variables in the composed system
+var_dict = Dict(string(var) => var for var in vars)
+pols = ["O3", "OH", "NO", "NO2", "CH4", "CH3O2", "CO","CH3OOH", "CH3O", "DMS", "SO2", "ISOP"]
+var_names_p = ["superfast₊$(v)(t)" for v in pols]
 
-Here is a plot that makes it easier to see what's going on for each species:
-```@example 2
 pp = []
-for (i, v) in enumerate(states(sys_noemis))
+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", )

docs/src/geoschem/overview.md

@@ -11,7 +11,7 @@ This mechanism is the result of many journal articles which are cited in API doc
 
 ## System overview
 
-First, let's initialize the model and inspect the first few reactions in the mechanism:
+First, let's initialize the model:
 
 ```@example 1
 using GasChem, EarthSciMLBase
@@ -21,28 +21,26 @@ using Unitful, Plots
 tspan = (0.0, 360.0)
 @variables t [unit = u"s", description = "Time"]
 gc = GEOSChemGasPhase(t)
-
-equations(gc.rxn_sys)[1:5]
 ```
 
 We can also look at the first few equations after converting the reaction network to a system of ODEs:
 
 ```@example 1
-equations(gc.sys)[1:5]
+equations(gc)[1:5]
 ```
 
 You can see that each reaction has a rate constant; rate constants are specified at the end of the list of equations:
 
 ```@example 1
-equations(gc.sys)[end-3:end]
+equations(gc)[end-3:end]
 ```
 
 ## Simulation
 
 Now, let's run a simulation and plot the results:
 
 ```@example 1
-sys = structural_simplify(get_mtk(gc))
+sys = structural_simplify(gc)
 vals = ModelingToolkit.get_defaults(sys)
 for k in setdiff(states(sys),keys(vals))
     vals[k] = 0 # Set variables with no default to zero.

docs/src/geoschem/params.md

@@ -12,10 +12,9 @@ using Unitful, Plots
 
 tspan = (0.0, 60.0*60*24*4) # 4 day simulation
 @variables t #[unit = u"s", description = "Time"]
-sys = get_mtk(GEOSChemGasPhase(t))
 
 # Run a simulation with constant temperature and pressure.
-sys = structural_simplify(sys)
+sys = structural_simplify(GEOSChemGasPhase(t))
 vals = ModelingToolkit.get_defaults(sys)
 for k in setdiff(states(sys),keys(vals))
     vals[k] = 0 # Set variables with no default to zero.
@@ -24,7 +23,7 @@ prob = ODEProblem(sys, vals, tspan, vals)
 sol1 = solve(prob, AutoTsit5(Rosenbrock23()))
 
 # Now, convert parameters to variables so we can change them over time.
-sys2 = param_to_var(get_mtk(GEOSChemGasPhase(t)), :T, :num_density)
+sys2 = param_to_var(GEOSChemGasPhase(t), :T, :num_density)
 
 # Vary temperature and pressure over time.
 @unpack T, num_density = sys2
@@ -63,7 +62,7 @@ Here is a list of all of the model parameters:
 ```@example 1
 using GasChem, DataFrames, EarthSciMLBase, ModelingToolkit, Unitful
 @variables t [unit = u"s", description = "Time"]
-gc = structural_simplify(get_mtk(GEOSChemGasPhase(t)))
+gc = structural_simplify(GEOSChemGasPhase(t))
 vars = parameters(gc)
 DataFrame(
         :Name => [string(Symbolics.tosymbol(v, escape=false)) for v ∈ vars],

docs/src/geoschem/states.md

@@ -5,7 +5,7 @@ Here is a list of the chemical species in the mechanism:
 ```@example 1
 using GasChem, DataFrames, EarthSciMLBase, ModelingToolkit, Unitful
 @variables t [unit = u"s", description = "Time"]
-gc = structural_simplify(get_mtk(GEOSChemGasPhase(t)))
+gc = structural_simplify(GEOSChemGasPhase(t))
 vars = states(gc)
 DataFrame(
         :Name => [string(Symbolics.tosymbol(v, escape=false)) for v ∈ vars],

docs/src/superfast.md

@@ -16,25 +16,31 @@ using Catalyst
 @parameters t [unit = u"s", description = "Time"]
 model = SuperFast(t)
 
-model.rxn_sys
+model
 ```
 
 We can also look at them as a graph:
 
 ```@example 1
-Graph(model.rxn_sys)
+Graph(model)
 ```
 
 
 Before running any simulations with the model, we need to convert it into a system of differential equations. We can solve it using the default values for variables and parameters. However, by using the ```@unpack``` command, we can assign new values to specific variables and parameters, allowing for simulations under varied conditions.
 
 We can visualize the differential equation version of the system as follows:
 ```@example 1
+model = couple(model) # Convert the ReactionSystem in Catalyst.jl to CoupledSystem in EarthSciMLBase.jl
 sys = structural_simplify(get_mtk(model))
-@unpack O3, T = sys
+
+vars = states(sys)  # Give you the variables in the system
+var_dict = Dict(string(var) => var for var in vars)
+pars = parameters(sys) # Give you the parameters in the system
+par_dict = Dict(string(par) => par for par in pars)
+
 tspan = (0.0, 3600*24)
-u0 = [O3 => 15] # Change the initial concentration of O₃ to 15 ppb
-p0 = [T => 293] # temperature = 293K
+u0 = [var_dict["superfast₊O3(t)"] => 15] # Change the initial concentration of O₃ to 15 ppb
+p0 = [par_dict["superfast₊T"] => 293] # temperature = 293K
 prob = ODEProblem(sys, u0, tspan, p0)
 
 equations(sys)
@@ -53,17 +59,13 @@ plot(sol, ylim = (0,50), xlabel = "Time", ylabel = "Concentration (ppb)", legend
 The species included in the superfast model are: O₃, OH, HO₂, O₂, NO, NO₂, CH₄, CH₃O₂, H₂O, CH₂O, CO, CH₃OOH, CH₃O, DMS, SO₂, ISOP, O₁d, H₂O₂.
 
 The parameters in the model that are not constant are the photolysis reaction rates ```jO31D```, ```j2OH```, ```jH2O2```, ```jNO2```, ```jCH2Oa```, ```jCH3OOH``` and temperature ```T```
-```@julia
-states(sys) # Give you the variables in the system
-parameters(sys) # Give you the parameters in the system
-```
 
 Let's run some simulation with different values for parameter ```T```.
 ```@example 1
-p1 = [T => 273]
-p2 = [T => 298]
+p1 = [par_dict["superfast₊T"] => 273]
+p2 = [par_dict["superfast₊T"] => 298]
 sol1 = solve(ODEProblem(sys, [], tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0)
 sol2 = solve(ODEProblem(sys, [], tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0)
 
-plot([sol1[O3],sol2[O3]], label = ["T=273K" "T=298K"], title = "Change of O3 concentration at different temperatures", xlabel="Time (second)", ylabel="concentration (ppb)")
+plot([sol1[var_dict["superfast₊O3(t)"]],sol2[var_dict["superfast₊O3(t)"]]], label = ["T=273K" "T=298K"], title = "Change of O3 concentration at different temperatures", xlabel="Time (second)", ylabel="concentration (ppb)")
 ```

src/geoschem_fullchem.jl

@@ -51,7 +51,6 @@ Moch et al, JGR, https, * Moch2020 # //doi.org/10.1029/2020JD032706, 2020.
 * Wolfe2012:    Wolfe et al., Phys. Chem. Chem. Phys., doi: 10.1039/C2CP40388A, 2012.
 * Xie2013:      Xie et al., Atmos. Chem. Phys., doi:10.5194/acp-13-8439-2013, 2013.
 """
-
 function GEOSChemGasPhase(t)
 
     # Create reaction rate constant system constructors