Update to MTKv9

docs/Project.toml

@@ -4,23 +4,23 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 EarthSciData = "a293c155-435f-439d-9c11-a083b6b47337"
 EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0"
 GasChem = "58070593-4751-4c87-a5d1-63807d11d76c"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 
-ModelingToolkit = "8"
+ModelingToolkit = "9"
 Documenter = "1"
-GasChem = "0.6"
+DynamicQuantities = "0.13"
+GasChem = "0.7"
 DifferentialEquations = "7"
-EarthSciMLBase = "0.15"
+EarthSciMLBase = "0.16"
 DataFrames = "1"
-Catalyst = "13"
-Unitful = "1"
+Catalyst = "14"
 Plots = "1"
 
 

docs/src/composing_models.md

@@ -9,22 +9,21 @@ Here is the complete example of composing, visualizing and solving the SuperFast
 model and the Fast-JX model, with explanation to follow:
 
 ```julia 
-using EarthSciMLBase, GasChem, ModelingToolkit, Dates, Unitful, DifferentialEquations
+using EarthSciMLBase, GasChem, ModelingToolkit, Dates, DynamicQuantities, DifferentialEquations
+using ModelingToolkit:t
 
-
-@parameters t [unit = u"s", description = "Time"]
-composed_ode = couple(SuperFast(t), FastJX(t)) # Compose two models use the "couple" function
+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 = structural_simplify(get_mtk(composed_ode)) # Define the coupled system  
+sys = structural_simplify(convert(ODESystem, composed_ode)) # Define the coupled system  
 sol = solve(ODEProblem(sys, [], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # Solve the coupled system
 ```
 
 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. 
 
 ```julia
-vars = states(sys)  # Get the variables in the composed system
+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", "CH4", "CH3O2", "CO","CH3OOH", "CH3O", "DMS", "SO2", "ISOP"]
 var_names_p = ["superfast₊$(v)(t)" for v in pols]
@@ -33,17 +32,17 @@ x_t = unix2datetime.(sol[t]) # Convert from unixtime to date time for visualizin
 ```
 Then, we could plot the results as:
 ```@setup 1
-using EarthSciMLBase, GasChem, ModelingToolkit, Dates, Unitful, DifferentialEquations
+using EarthSciMLBase, GasChem, ModelingToolkit, Dates, DynamicQuantities, DifferentialEquations
+using ModelingToolkit:t
 
-@parameters t [unit = u"s"]
-composed_ode = couple(SuperFast(t), FastJX(t)) # Compose two models use the "couple" function
+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 = structural_simplify(get_mtk(composed_ode)) # Define the coupled system  
+sys = structural_simplify(convert(ODESystem, composed_ode)) # Define the coupled system  
 sol = solve(ODEProblem(sys, [], tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # Solve the coupled system
 
-vars = states(sys)  # Get the variables in the composed system
+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", "CH4", "CH3O2", "CO","CH3OOH", "CH3O", "DMS", "SO2", "ISOP"]
 var_names_p = ["SuperFast₊$(v)(t)" for v in pols]
@@ -68,27 +67,27 @@ Here's a simple example:
 ```@example 2 
 using GasChem, EarthSciData # This will trigger the emission extension
 using Dates, ModelingToolkit, DifferentialEquations, EarthSciMLBase
-using Plots, Unitful
+using Plots, DynamicQuantities
+using ModelingToolkit:t
 
-@parameters t [unit = u"s", description = "Time"]
 @parameters lat = 40 
 @parameters lon = -97 
 @parameters lev = 1
-emis = NEI2016MonthlyEmis("mrggrid_withbeis_withrwc", t, lon, lat, lev; dtype=Float64) 
+emis = NEI2016MonthlyEmis("mrggrid_withbeis_withrwc", lon, lat, lev; 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.
+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 = structural_simplify(get_mtk(model_noemis))
-sys_withemis = structural_simplify(get_mtk(model_withemis))
+sys_noemis = structural_simplify(convert(ODESystem, model_noemis))
+sys_withemis = structural_simplify(convert(ODESystem, model_withemis))
 
 start = Dates.datetime2unix(Dates.DateTime(2016, 5, 1))
 tspan = (start, start+3*24*3600)
 sol_noemis = solve(ODEProblem(sys_noemis, [], tspan, []), AutoTsit5(Rosenbrock23()))
 sol_withemis = solve(ODEProblem(sys_withemis, [], tspan, []), AutoTsit5(Rosenbrock23()))
 
-vars = states(sys_noemis)  # Get the variables in the composed system
+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", "CH4", "CH3O2", "CO","CH3OOH", "CH3O", "DMS", "SO2", "ISOP"]
 var_names_p = ["SuperFast₊$(v)(t)" for v in pols]

docs/src/geoschem/overview.md

@@ -16,11 +16,11 @@ First, let's initialize the model and we can also look at the first few ODE equa
 ```@example 1
 using GasChem, EarthSciMLBase
 using DifferentialEquations, ModelingToolkit
-using Unitful, Plots
+using DynamicQuantities, Plots
+using ModelingToolkit:t
 
 tspan = (0.0, 360.0)
-@variables t [unit = u"s", description = "Time"]
-gc = GEOSChemGasPhase(t)
+gc = GEOSChemGasPhase()
 equations(gc)[1:5]
 ```
 
@@ -37,11 +37,11 @@ Now, let's run a simulation and plot the results:
 ```@example 1
 sys = structural_simplify(gc)
 vals = ModelingToolkit.get_defaults(sys)
-for k in setdiff(states(sys),keys(vals))
+for k in setdiff(unknowns(sys),keys(vals))
     vals[k] = 0 # Set variables with no default to zero.
 end
 prob = ODEProblem(sys, vals, tspan, vals)
 sol = solve(prob, AutoTsit5(Rosenbrock23()))
 plot(sol, legend = :outertopright, xlabel = "Time (s)", 
-        ylabel = "Concentration (nmol/mol)")
+        ylabel = "Concentration (ppb)")
 ```

docs/src/geoschem/params.md

@@ -8,27 +8,27 @@ We can explore what happens when we change them:
 ```@example
 using GasChem, EarthSciMLBase
 using DifferentialEquations, ModelingToolkit
-using Unitful, Plots
+using DynamicQuantities, Plots
+using ModelingToolkit:t
 
 tspan = (0.0, 60.0*60*24*4) # 4 day simulation
-@variables t #[unit = u"s", description = "Time"]
 
 # Run a simulation with constant temperature and pressure.
-sys = structural_simplify(GEOSChemGasPhase(t))
+sys = structural_simplify(GEOSChemGasPhase())
 vals = ModelingToolkit.get_defaults(sys)
-for k in setdiff(states(sys),keys(vals))
+for k in setdiff(unknowns(sys),keys(vals))
     vals[k] = 0 # Set variables with no default to zero.
 end
 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(GEOSChemGasPhase(t), :T, :num_density)
+sys2 = param_to_var(GEOSChemGasPhase(), :T, :num_density)
 
 # Vary temperature and pressure over time.
 @unpack T, num_density = sys2
 @constants T_0 = 300 [unit=u"K"]
-@constants t_0 = 1 #[unit=u"s"]
+@constants t_0 = 1 [unit=u"s"]
 eqs = [
     T ~ T_0 + T_0 / 1.5 * sin(2π*t/t_0/(60*60*24)),
     num_density ~ 2.7e19 - 2.5e19*t/t_0/(60*60*24*4),
@@ -38,15 +38,15 @@ sys2 = extend(sys2,ODESystem(eqs, t; name=:var_T))
 # Run the simulation again.
 sys2 = structural_simplify(sys2)
 vals = ModelingToolkit.get_defaults(sys2)
-for k in setdiff(states(sys2),keys(vals))
+for k in setdiff(unknowns(sys2),keys(vals))
     vals[k] = 0 # Set variables with no default to zero.
 end
 prob = ODEProblem(sys2, vals, tspan, vals)
 sol2 = solve(prob, AutoTsit5(Rosenbrock23()))
 
 # Plot the results
 p1 = plot(sol1.t, sol1[sys2.O3], xticks=:none, label="Constant T and P",
-        ylabel = "Concentration (nmol/mol)")
+        ylabel = "Concentration (ppb)")
 plot!(p1, sol2.t, sol2[sys2.O3], label="Varying T and P")
 
 p2 = plot(sol2.t, sol2[sys2.T], label = "T (K)", xticks=:none)
@@ -60,9 +60,9 @@ plot(p1, p2, p3, layout=grid(3, 1, heights=[0.7, 0.15, 0.15]))
 Here is a list of all of the model parameters:
 
 ```@example 1
-using GasChem, DataFrames, EarthSciMLBase, ModelingToolkit, Unitful
+using GasChem, DataFrames, EarthSciMLBase, ModelingToolkit, DynamicQuantities
 @variables t [unit = u"s", description = "Time"]
-gc = structural_simplify(GEOSChemGasPhase(t))
+gc = structural_simplify(GEOSChemGasPhase())
 vars = parameters(gc)
 DataFrame(
         :Name => [string(Symbolics.tosymbol(v, escape=false)) for v ∈ vars],

docs/src/geoschem/states.md

@@ -3,10 +3,10 @@
 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(GEOSChemGasPhase(t))
-vars = states(gc)
+using GasChem, DataFrames, EarthSciMLBase, ModelingToolkit, DynamicQuantities
+using ModelingToolkit:t
+gc = structural_simplify(GEOSChemGasPhase())
+vars = unknowns(gc)
 DataFrame(
         :Name => [string(Symbolics.tosymbol(v, escape=false)) for v ∈ vars],
         :Units => [ModelingToolkit.get_unit(v) for v ∈ vars],

docs/src/superfast.md

@@ -11,25 +11,25 @@ First, we can look at the reaction equations:
 
 ```@example 1
 using GasChem, EarthSciMLBase, ModelingToolkit
-using Unitful, DifferentialEquations
+using DynamicQuantities, DifferentialEquations
 using Catalyst
 using Plots
+using ModelingToolkit:t
 
-@parameters t [unit = u"s", description = "Time"]
-model = SuperFast(t)
+model = SuperFast()
 ```
 
 We can also look at them as a graph:
 
 ```@example 1
-Graph(SuperFast(t, rxn_sys=true))
+Graph(SuperFast(;name=:SuperFast, rxn_sys=true))
 ```
 
 ## Variables and parameters
 The chemical species included in the superfast model are:
 
 ```@example 1
-vars = states(model)
+vars = unknowns(model)
 using DataFrames
 DataFrame(
         :Name => [string(Symbolics.tosymbol(v, escape=false)) for v ∈ vars],

ext/EarthSciDataExt.jl

@@ -1,31 +1,31 @@
 module EarthSciDataExt
-using GasChem, EarthSciData, ModelingToolkit, EarthSciMLBase, Unitful
+using GasChem, EarthSciData, ModelingToolkit, EarthSciMLBase, DynamicQuantities
+@register_unit ppb 1
 
 function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, e::EarthSciData.NEI2016MonthlyEmisCoupler)
     c, e = c.sys, e.sys
 
     @constants(
-        MW_NO2 = 46.0055, [unit = u"g/mol", description="NO2 molar mass"],
-        MW_NO = 30.01, [unit = u"g/mol", description="NO molar mass"],
-        MW_FORM = 30.0260, [unit = u"g/mol", description="Formaldehyde molar mass"],
-        MW_CH4 = 16.0425, [unit = u"g/mol", description="Methane molar mass"],
-        MW_CO = 28.0101, [unit = u"g/mol", description="Carbon monoxide molar mass"],
-        MW_SO2 = 64.0638, [unit = u"g/mol", description="Sulfur dioxide molar mass"],
-        MW_ISOP = 68.12, [unit = u"g/mol", description="Isoprene molar mass"],
-
-        gperkg = 1e3, [unit = u"g/kg", description="Conversion factor from kg to g"],
-        nmolpermol = 1e9, [unit = u"nmol/mol", description="Conversion factor from mol to nmol"],
+        MW_NO2 = 46.0055*10^-3, [unit = u"kg/mol", description="NO2 molar mass"],
+        MW_NO = 30.01*10^-3, [unit = u"kg/mol", description="NO molar mass"],
+        MW_FORM = 30.0260*10^-3, [unit = u"kg/mol", description="Formaldehyde molar mass"],
+        MW_CH4 = 16.0425*10^-3, [unit = u"kg/mol", description="Methane molar mass"],
+        MW_CO = 28.0101*10^-3, [unit = u"kg/mol", description="Carbon monoxide molar mass"],
+        MW_SO2 = 64.0638*10^-3, [unit = u"kg/mol", description="Sulfur dioxide molar mass"],
+        MW_ISOP = 68.12*10^-3, [unit = u"kg/mol", description="Isoprene molar mass"],
+
+        nmolpermol = 1e9, [unit = u"ppb", description="noml/mol, Conversion factor from mol to nmol"],
         R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"],
     )
 
     # Emissions are in units of "kg/m3/s" and need to be converted to "ppb/s" or "nmol/mol/s".
     # To do this we need to convert kg of emissions to nmol of emissions,
     # and we need to convert m3 of air to mol of air.
-    # nmol_emissions = kg_emissions * gperkg / MW_emission * nmolpermol = kg * g/kg / g/mol * nmol/mol = nmol
+    # nmol_emissions = kg_emissions * gperkg / MW_emission * nmolpermol = kg / kg/mol * nmol/mol = nmol
     # mol_air = m3_air / R / T * P = m3 / (m3*Pa/mol/K) / K * Pa = mol
     # So, the overall conversion is:
-    # nmol_emissions / mol_air = (kg_emissions * gperkg / MW_emission * nmolpermol) / (m3_air / R / T * P)
-    uconv = gperkg * nmolpermol * R * c.T / c.P # Conversion factor with MW factored out.
+    # nmol_emissions / mol_air = (kg_emissions / MW_emission * nmolpermol) / (m3_air / R / T * P)
+    uconv =  nmolpermol * R * c.T / c.P # Conversion factor with MW factored out.
     operator_compose(convert(ODESystem, c), e, Dict(
         c.NO2 => e.NO2 => uconv / MW_NO2,
         c.NO => e.NO => uconv / MW_NO,
@@ -40,11 +40,10 @@ end
 function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, g::EarthSciData.GEOSFPCoupler)
     c, g = c.sys, g.sys
 
-    @constants PaPerhPa = 100 [unit = u"Pa/hPa", description="Conversion factor from hPa to Pa"]
     c = param_to_var(c, :T, :P)
     ConnectorSystem([
             c.T ~ g.I3₊T,
-            c.P ~ g.P * PaPerhPa,
+            c.P ~ g.P,
         ], c, g)
 end
 

src/Fast-JX.jl

@@ -127,13 +127,13 @@ function j_mean(σ_interp, ϕ::Float32, time::T2, lat::T, long::T, Temperature::
 end
 
 # Dummy functions for unit validation. Basically ModelingToolkit 
-# will call the function with a Unitful.Quantity or an integer to 
+# will call the function with a DynamicQuantities.Quantity or an integer to 
 # get information about the type and units of the output.
-j_mean_H2O2(t, lat, long, T::Unitful.Quantity) = 0.0u"s^-1"
-j_mean_CH2Oa(t, lat, long, T::Unitful.Quantity) = 0.0u"s^-1"
-j_mean_CH3OOH(t, lat, long, T::Unitful.Quantity) = 0.0u"s^-1"
-j_mean_NO2(t, lat, long, T::Unitful.Quantity) = 0.0u"s^-1"
-j_mean_o31D(t, lat, long, T::Unitful.Quantity) = 0.0u"s^-1"
+j_mean_H2O2(t, lat, long, T::DynamicQuantities.Quantity) = 0.0u"s^-1"
+j_mean_CH2Oa(t, lat, long, T::DynamicQuantities.Quantity) = 0.0u"s^-1"
+j_mean_CH3OOH(t, lat, long, T::DynamicQuantities.Quantity) = 0.0u"s^-1"
+j_mean_NO2(t, lat, long, T::DynamicQuantities.Quantity) = 0.0u"s^-1"
+j_mean_o31D(t, lat, long, T::DynamicQuantities.Quantity) = 0.0u"s^-1"
 
 j_mean_H2O2(t, lat, long, T) = j_mean(σ_H2O2_interp, ϕ_H2O2_jx, t, lat, long, T)
 @register_symbolic j_mean_H2O2(t, lat, long, T)
@@ -157,11 +157,10 @@ Description: This is a box model used to calculate the photolysis reaction rate
 Build Fast-JX model
 # Example
 ``` julia
-    @parameters t [unit = u"s"]
-    fj = FastJX(t)
+    fj = FastJX()
 ```
 """
-function FastJX(t; name=:FastJX)
+function FastJX(;name=:FastJX)
     @parameters T = 298.0 [unit = u"K", description = "Temperature"]
     @parameters lat = 40.0 [description = "Latitude (Degrees)"]
     @parameters long = -97.0 [description = "Longitude (Degrees)"]

src/GasChem.jl

@@ -4,9 +4,12 @@ using EarthSciMLBase
 using ModelingToolkit
 using Catalyst
 using Dates
-using Unitful
+using DynamicQuantities
 using StaticArrays
 using Interpolations
+using ModelingToolkit:t
+
+@register_unit ppb 1
 
 include("SuperFast.jl")
 include("geoschem_ratelaws.jl")