From 41dbb38132da69ea3f881cdae13ed29d91c7d036 Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Wed, 10 Jul 2024 22:43:18 +0800 Subject: [PATCH] Update documentation --- docs/src/composing_models.md | 38 +++++++++++++++++------------------ docs/src/geoschem/overview.md | 10 ++++----- docs/src/geoschem/params.md | 7 +++---- docs/src/geoschem/states.md | 2 +- docs/src/superfast.md | 26 +++++++++++++----------- src/geoschem_fullchem.jl | 1 - 6 files changed, 41 insertions(+), 43 deletions(-) diff --git a/docs/src/composing_models.md b/docs/src/composing_models.md index 9a91821d..c5ce1f08 100644 --- a/docs/src/composing_models.md +++ b/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", ) diff --git a/docs/src/geoschem/overview.md b/docs/src/geoschem/overview.md index 63d3dffd..fa48d559 100644 --- a/docs/src/geoschem/overview.md +++ b/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,20 +21,18 @@ 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 @@ -42,7 +40,7 @@ equations(gc.sys)[end-3:end] 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. diff --git a/docs/src/geoschem/params.md b/docs/src/geoschem/params.md index 3481cce0..23e51756 100644 --- a/docs/src/geoschem/params.md +++ b/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], diff --git a/docs/src/geoschem/states.md b/docs/src/geoschem/states.md index a526860f..40f4f1ec 100644 --- a/docs/src/geoschem/states.md +++ b/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], diff --git a/docs/src/superfast.md b/docs/src/superfast.md index 42466122..e9602455 100644 --- a/docs/src/superfast.md +++ b/docs/src/superfast.md @@ -16,13 +16,13 @@ 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) ``` @@ -30,11 +30,17 @@ Before running any simulations with the model, we need to convert it into a syst 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)") ``` \ No newline at end of file diff --git a/src/geoschem_fullchem.jl b/src/geoschem_fullchem.jl index 82a1d610..78d49b55 100644 --- a/src/geoschem_fullchem.jl +++ b/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