From 3a25667e0427d9799a3629aab2d8d47dd404440e Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Tue, 20 Aug 2024 22:24:31 -0700 Subject: [PATCH 01/15] Update with EarthSciMLBasev0.15 --- Project.toml | 8 ++++-- ext/EarthSciDataExt.jl | 51 ++++++++++++++++++++++++++++++++++ ext/GasChemExt.jl | 63 +++++++++++++++--------------------------- src/dry_deposition.jl | 44 +++++++++++++++-------------- src/wet_deposition.jl | 38 +++++++++++++------------ test/connector_test.jl | 31 +++++++++++++++++---- 6 files changed, 149 insertions(+), 86 deletions(-) create mode 100644 ext/EarthSciDataExt.jl diff --git a/Project.toml b/Project.toml index 2e2d7eec..f70640de 100644 --- a/Project.toml +++ b/Project.toml @@ -11,18 +11,21 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [weakdeps] +EarthSciData = "a293c155-435f-439d-9c11-a083b6b47337" GasChem = "58070593-4751-4c87-a5d1-63807d11d76c" [extensions] GasChemExt = "GasChem" +EarthSciDataExt = "EarthSciData" [compat] DifferentialEquations = "7" -EarthSciMLBase = "0.10" +EarthSciMLBase = "0.15" GasChem = "0.6" IfElse = "0.1" ModelingToolkit = "8" @@ -33,9 +36,10 @@ julia = "1" [extras] Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +EarthSciData = "a293c155-435f-439d-9c11-a083b6b47337" GasChem = "58070593-4751-4c87-a5d1-63807d11d76c" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "SafeTestsets", "GasChem", "Dates"] +test = ["Test", "SafeTestsets", "GasChem", "Dates", "EarthSciData"] diff --git a/ext/EarthSciDataExt.jl b/ext/EarthSciDataExt.jl new file mode 100644 index 00000000..08babae6 --- /dev/null +++ b/ext/EarthSciDataExt.jl @@ -0,0 +1,51 @@ +module EarthSciDataExt + +using AtmosphericDeposition, EarthSciData, EarthSciMLBase, Unitful, ModelingToolkit + +function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler, g::EarthSciData.GEOSFPCoupler) + d, g = d.sys, g.sys + + @constants( + PaPerhPa = 100, [unit = u"Pa/hPa", description="Conversion factor from hPa to Pa"], + MW_air = 29, [unit = u"g/mol", description="dry air molar mass"], + kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"], + R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"], + ) + + # ρA in DrydepositionG(t) are in units of "kg/m3". + # ρ = P*M/(R*T)= Pa*(g/mol)/(m3*Pa/mol/K*K) = g/m3 + # Overall, ρA = P*M/(R*T)*kgperg + + d = param_to_var(d, :T, :z, :z₀, :u_star, :G, :ρA) + ConnectorSystem([ + d.T ~ g.I3₊T, + d.z ~ g.A1₊PBLH, + d.z₀ ~ g.A1₊Z0M, + d.u_star ~ g.A1₊USTAR, + d.G ~ g. A1₊SWGDN, + d.ρA ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air, + ], d, g) +end + +function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g::EarthSciData.GEOSFPCoupler) + d, g = d.sys, g.sys + + @constants( + PaPerhPa = 100, [unit = u"Pa/hPa", description="Conversion factor from hPa to Pa"], + MW_air = 29, [unit = u"g/mol", description="dry air molar mass"], + kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"], + R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"], + ) + + # ρ_air in Wetdeposition(t) are in units of "kg/m3". + # ρ = P*M/(R*T)= Pa*(g/mol)/(m3*Pa/mol/K*K) = g/m3 + # Overall, ρ_air = P*M/(R*T)*kgperg + + d = param_to_var(d, :cloudFrac, :ρ_air) + ConnectorSystem([ + d.cloudFrac ~ g.A3cld₊CLOUD, + d.ρ_air ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air, + ], d, g) +end + +end \ No newline at end of file diff --git a/ext/GasChemExt.jl b/ext/GasChemExt.jl index 00c3b792..f7420fca 100644 --- a/ext/GasChemExt.jl +++ b/ext/GasChemExt.jl @@ -2,48 +2,31 @@ module GasChemExt using AtmosphericDeposition, GasChem, EarthSciMLBase, Unitful, ModelingToolkit -export register_couplings_ext - -# Use a global flag to track initialization and ensure that the coupling between SuperFast and NEI Emission is only initialized once -const couplings_registered_ext = Ref(false) - -function register_couplings_ext() - if couplings_registered_ext[] - println("Couplings have already been registered.") - return - end - - println("Registering couplings in ext") - @parameters t [unit = u"s"] - - register_coupling(SuperFast(t), Wetdeposition(t)) do c, e - operator_compose(convert(ODESystem, c), e, Dict( - c.SO2 => e.SO2, - c.NO2 => e.NO2, - c.O3 => e.O3, - c.CH4 => e.CH4, - c.CO => e.CO, - c.DMS => e.DMS, - c.ISOP => e.ISOP)) - end - - register_coupling(SuperFast(t), DrydepositionG(t)) do c, e - operator_compose(convert(ODESystem, c), e, Dict( - c.SO2 => e.SO2, - c.NO2 => e.NO2, - c.O3 => e.O3, - c.NO => e.NO, - c.H2O2 => e.H2O2, - c.CH2O => e.CH2O)) - end - - couplings_registered_ext[] = true - println("Coupling registry after registration: ", EarthSciMLBase.coupling_registry) +function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDeposition.DrydepositionGCoupler) + c, d = c.sys, d.sys + + operator_compose(convert(ODESystem, c), d, Dict( + c.SO2 => d.SO2, + c.NO2 => d.NO2, + c.O3 => d.O3, + c.NO => d.NO, + c.H2O2 => d.H2O2, + c.CH2O => d.CH2O, + )) end -function __init__() - println("Initializing EmissionExt module") - register_couplings_ext() +function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDeposition.WetdepositionCoupler) + c, d = c.sys, d.sys + + operator_compose(convert(ODESystem, c), d, Dict( + c.SO2 => d.SO2, + c.NO2 => d.NO2, + c.O3 => d.O3, + c.CH4 => d.CH4, + c.CO => d.CO, + c.DMS => d.DMS, + c.ISOP => d.ISOP, + )) end end \ No newline at end of file diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 9698f332..b54789ae 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -220,12 +220,9 @@ end defaults = [g => 9.81, κ => 0.4, k => 1.3806488e-23, M_air => 28.97e-3, R => 8.3144621, unit_T => 1, unit_convert_mu => 1, T_unitless => 1, unit_dH2O => 1, unit_m => 1, G_unitless => 1, Rc_unit => 1, unit_v => 1] -# Add unit "ppb" to Unitful -module MyUnits -using Unitful -@unit ppb "ppb" Number 1 / 1000000000 false +struct DrydepositionGCoupler + sys end -Unitful.register(MyUnits) """ Description: This is a box model used to calculate the gas species concentration rate changed by dry deposition. @@ -236,28 +233,32 @@ Build Drydeposition model (gas) d = DrydepositionG(t) ``` """ -function DrydepositionG(t) +function DrydepositionG(t; name=:DrydepositionG) iSeason = 1 iLandUse = 10 rain = false dew = false - @parameters z = 50 [unit = u"m", description = "top of the surface layer"] - @parameters z₀ = 0.04 [unit = u"m", description = "roughness lenght"] - @parameters u_star = 0.44 [unit = u"m/s", description = "friction velocity"] - @parameters L = 0 [unit = u"m", description = "Monin-Obukhov length"] - @parameters ρA = 1.2 [unit = u"kg*m^-3", description = "air density"] - @parameters G = 300 [unit = u"W*m^-2", description = "solar irradiation"] - @parameters T = 298 [unit = u"K", description = "temperature"] - @parameters θ = 0 [description = "slope of the local terrain, in unit radians"] + params = @parameters( + z = 50, [unit = u"m", description = "top of the surface layer"], + z₀ = 0.04, [unit = u"m", description = "roughness lenght"], + u_star = 0.44, [unit = u"m/s", description = "friction velocity"], + L = 0, [unit = u"m", description = "Monin-Obukhov length"], + ρA = 1.2, [unit = u"kg*m^-3", description = "air density"], + G = 300, [unit = u"W*m^-2", description = "solar irradiation"], + T = 298, [unit = u"K", description = "temperature"], + θ = 0, [description = "slope of the local terrain, in unit radians"], + ) D = Differential(t) - @variables SO2(t) [unit = u"ppb"] - @variables O3(t) [unit = u"ppb"] - @variables NO2(t) [unit = u"ppb"] - @variables NO(t) [unit = u"ppb"] - @variables H2O2(t) [unit = u"ppb"] - @variables CH2O(t) [unit = u"ppb"] + vars = @variables( + SO2(t), [unit = u"nmol/mol"], + O3(t),[unit = u"nmol/mol"], + NO2(t), [unit = u"nmol/mol"], + NO(t), [unit = u"nmol/mol"], + H2O2(t), [unit = u"nmol/mol"], + CH2O(t), [unit = u"nmol/mol"], + ) eqs = [ D(SO2) ~ -DryDepGas(z, z₀, u_star, L, ρA, So2Data, G, T, θ, iSeason, iLandUse, rain, dew, true, false) / z * SO2 @@ -268,6 +269,7 @@ function DrydepositionG(t) D(CH2O) ~ -DryDepGas(z, z₀, u_star, L, ρA, HchoData, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * CH2O ] - ODESystem(eqs, t, [SO2, O3, NO2, NO, H2O2, CH2O], [z, z₀, u_star, L, ρA, G, T, θ]; name=:DrydepositionG) + ODESystem(eqs, t, vars, params; name=name, + metadata=Dict(:coupletype => DrydepositionGCoupler)) end diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index 2efd507e..d143195c 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -41,12 +41,9 @@ function WetDeposition(cloudFrac, qrain, ρ_air, Δz) end wd_defaults = [A_wd => 5.2, ρwater => 1000.0, Vdr => 5.0] -# Add unit "ppb" to Unitful -module myUnits -using Unitful -@unit ppb "ppb" Number 1 / 1000000000 false +struct WetdepositionCoupler + sys end -Unitful.register(myUnits) """ Description: This is a box model used to calculate wet deposition based on formulas at EMEP model. @@ -57,21 +54,25 @@ Build Wetdeposition model wd = Wetdeposition(t) ``` """ -function Wetdeposition(t) - @parameters cloudFrac = 0.5 [description = "fraction of grid cell covered by clouds"] - @parameters qrain = 0.5 [description = "rain mixing ratio"] - @parameters ρ_air = 1.204 [unit = u"kg*m^-3", description = "air density"] - @parameters Δz = 200 [unit = u"m", description = "fall distance"] +function Wetdeposition(t; name=:Wetdeposition) + params = @parameters( + cloudFrac = 0.5, [description = "fraction of grid cell covered by clouds"], + qrain = 0.5, [description = "rain mixing ratio"], + ρ_air = 1.204, [unit = u"kg*m^-3", description = "air density"], + Δz = 200, [unit = u"m", description = "fall distance"], + ) D = Differential(t) - @variables SO2(t) [unit = u"ppb"] - @variables O3(t) [unit = u"ppb"] - @variables NO2(t) [unit = u"ppb"] - @variables CH4(t) [unit = u"ppb"] - @variables CO(t) [unit = u"ppb"] - @variables DMS(t) [unit = u"ppb"] - @variables ISOP(t) [unit = u"ppb"] + vars = @variables( + SO2(t), [unit = u"nmol/mol"], + O3(t), [unit = u"nmol/mol"], + NO2(t), [unit = u"nmol/mol"], + CH4(t), [unit = u"nmol/mol"], + CO(t), [unit = u"nmol/mol"], + DMS(t), [unit = u"nmol/mol"], + ISOP(t), [unit = u"nmol/mol"], + ) eqs = [ D(SO2) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[2] * SO2 @@ -83,5 +84,6 @@ function Wetdeposition(t) D(ISOP) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * ISOP ] - ODESystem(eqs, t, [SO2, O3, NO2, CH4, CO, DMS, ISOP], [cloudFrac, qrain, ρ_air, Δz]; name=:WetDeposition) + ODESystem(eqs, t, vars, params; name=name, + metadata=Dict(:coupletype => WetdepositionCoupler)) end diff --git a/test/connector_test.jl b/test/connector_test.jl index 5c813e74..7ee7f918 100644 --- a/test/connector_test.jl +++ b/test/connector_test.jl @@ -1,16 +1,37 @@ using AtmosphericDeposition -using Test, Unitful, ModelingToolkit, GasChem, Dates, EarthSciMLBase +using Test, Unitful, ModelingToolkit, GasChem, Dates, EarthSciMLBase, EarthSciData -@testset "Connector" begin - ModelingToolkit.check_units(eqs...) = nothing +@testset "GasChemExt" begin start = Dates.datetime2unix(Dates.DateTime(2016, 5, 1)) - @parameters t + @parameters t [unit = u"s"] composed_ode = couple(SuperFast(t), FastJX(t), DrydepositionG(t), Wetdeposition(t)) tspan = (start, start+3600*24*3) sys = structural_simplify(get_mtk(composed_ode)) @test length(states(sys)) ≈ 18 eqs = string(equations(sys)) - wanteqs = ["Differential(t)(superfast₊O3(t)) ~ superfast₊DrydepositionG_ddt_O3ˍt(t) + superfast₊WetDeposition_ddt_O3ˍt(t)"] + wanteqs = ["Differential(t)(SuperFast₊O3(t)) ~ SuperFast₊DrydepositionG_ddt_O3ˍt(t) + SuperFast₊Wetdeposition_ddt_O3ˍt(t)"] @test contains(string(eqs), wanteqs[1]) end + +@testset "EarthSciDataExt" begin + @parameters t [unit = u"s"] + + @parameters lat = 40 + @parameters lon = -97 + @parameters lev = 1 + geosfp = GEOSFP("4x5", t) + + model = couple(SuperFast(t), FastJX(t), geosfp, Wetdeposition(t), DrydepositionG(t)) + + sys = structural_simplify(get_mtk(model)) + @test length(states(sys)) ≈ 18 + + eqs = string(equations(get_mtk(model))) + wanteq = "DrydepositionG₊G(t) ~ GEOSFP₊A1₊SWGDN(t)" + @test contains(eqs, wanteq) + wanteq = "Wetdeposition₊cloudFrac(t) ~ GEOSFP₊A3cld₊CLOUD(t)" + @test contains(eqs, wanteq) + wanted = "Wetdeposition₊ρA(t) ~ GEOSFP₊P * PaPerhPa/(GEOSFP₊I3₊T*R)*kgperg*MW_air" + @test contains(eqs, wanteq) +end \ No newline at end of file From 97d1540ed68dcea918527b6a56d87b845fab8d9b Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Mon, 26 Aug 2024 01:17:55 +0000 Subject: [PATCH 02/15] CompatHelper: add new compat entry for EarthSciData in [weakdeps] at version 0.9, (keep existing compat) --- Project.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f70640de..5cf6eaa2 100644 --- a/Project.toml +++ b/Project.toml @@ -20,11 +20,12 @@ EarthSciData = "a293c155-435f-439d-9c11-a083b6b47337" GasChem = "58070593-4751-4c87-a5d1-63807d11d76c" [extensions] -GasChemExt = "GasChem" EarthSciDataExt = "EarthSciData" +GasChemExt = "GasChem" [compat] DifferentialEquations = "7" +EarthSciData = "0.9" EarthSciMLBase = "0.15" GasChem = "0.6" IfElse = "0.1" From 46379366642c9cc22c12fee02c1d31ff4ae7e8b1 Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Fri, 6 Sep 2024 18:05:45 -0500 Subject: [PATCH 03/15] Chang iSeason & iLandUse as parameters and add connections with geosfp to get qrain and L --- Project.toml | 3 ++- docs/src/wesley1989.md | 18 ++++++++----- src/AtmosphericDeposition.jl | 1 + src/dry_deposition.jl | 12 ++++----- src/wesley1989.jl | 51 ++++++++++++++++++++---------------- test/drydep_test.jl | 6 ++--- 6 files changed, 53 insertions(+), 38 deletions(-) diff --git a/Project.toml b/Project.toml index f70640de..977787e1 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["EarthSciML authors and contributors"] version = "0.2.1" [deps] +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0" IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" @@ -20,8 +21,8 @@ EarthSciData = "a293c155-435f-439d-9c11-a083b6b47337" GasChem = "58070593-4751-4c87-a5d1-63807d11d76c" [extensions] -GasChemExt = "GasChem" EarthSciDataExt = "EarthSciData" +GasChemExt = "GasChem" [compat] DifferentialEquations = "7" diff --git a/docs/src/wesley1989.md b/docs/src/wesley1989.md index 95ffc6bf..c5404633 100644 --- a/docs/src/wesley1989.md +++ b/docs/src/wesley1989.md @@ -39,13 +39,19 @@ a. For season index ```iSeason```, there're five seasonal categories 4. Winter, snow on ground 5. Transitional -b. For land use index ```iLandUse```, there're five land use categories +b. For land use index ```iLandUse```, there're eleven land use categories - 1. Evergreen–needleleaf trees - 2. Deciduous broadleaf trees - 3. Grass - 4. Desert - 5. Shrubs and interrupted woodlands + 1. Urban land + 2. agricultural land + 3. range land + 4. deciduous forest + 5. coniferous forest + 6. mixed forest including wetland + 7. water, both salt and fresh + 8. barren land, mostly desert + 9. nonforested wetland + 10. mixed agricultural and range land + 11. rocky open areas with low-growing shrubs c. For gasData ```gasData```, we include the following species: diff --git a/src/AtmosphericDeposition.jl b/src/AtmosphericDeposition.jl index 8414bf38..a7ccfe4f 100644 --- a/src/AtmosphericDeposition.jl +++ b/src/AtmosphericDeposition.jl @@ -5,6 +5,7 @@ using StaticArrays using ModelingToolkit using EarthSciMLBase using IfElse +using DataInterpolations include("wesley1989.jl") include("dry_deposition.jl") diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index b54789ae..3d1234ff 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -185,13 +185,13 @@ irradiation [W m-2], Θ is the slope of the local terrain [radians], iSeason and dew and rain indicate whether there is dew or rain on the ground, and isSO2 and isO3 indicate whether the gas species of interest is either SO2 or O3, respectively. Based on Seinfeld and Pandis (2006) equation 19.2. """ -function DryDepGas(z, z₀, u_star, L, ρA, gasData::GasData, G, Ts, θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) +function DryDepGas(z, z₀, u_star, L, ρA, gasData::GasData, G, Ts, θ, iSeason, iLandUse, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) Ra = ra(z, z₀, u_star, L) μ = mu(Ts) Dg = dH2O(Ts) / gasData.Dh2oPerDx # Diffusivity of gas of interest [m2/s] Sc = sc(μ, ρA, Dg) Rb = RbGas(Sc, u_star) - Rc = WesleySurfaceResistance(gasData, G * G_unitless, (Ts * T_unitless - 273), θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) * Rc_unit + Rc = WesleySurfaceResistance(gasData, G * G_unitless, (Ts * T_unitless - 273), θ, iSeason, iLandUse, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) * Rc_unit return 1 / (Ra + Rb + Rc) end @@ -202,7 +202,7 @@ Dp is particle diameter [m], Ts is surface air temperature [K], P is pressure [P and iSeason and iLandUse are indexes for the season and land use. Based on Seinfeld and Pandis (2006) equation 19.7. """ -function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason::Int, iLandUse::Int) +function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason, iLandUse) Ra = ra(z, z₀, u_star, L) μ = mu(Ts) Cc = cc(Dp, Ts, P, μ) @@ -210,7 +210,7 @@ function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason: if iLandUse == 4 # dessert St = stSmooth(Vs, u_star, μ, ρA) else - St = stVeg(Vs, u_star, (A_table[iSeason, iLandUse] * 10^(-3)) * unit_m) + St = stVeg(Vs, u_star, (obtain_value(iSeason, iLandUse,A_table) * 10^(-3)) * unit_m) end D = dParticle(Ts, P, Dp, Cc, μ) Sc = sc(μ, ρA, D) @@ -234,11 +234,11 @@ Build Drydeposition model (gas) ``` """ function DrydepositionG(t; name=:DrydepositionG) - iSeason = 1 - iLandUse = 10 rain = false dew = false params = @parameters( + iSeason = 1, [description = "Index for season"], + iLandUse = 10, [description = "Index for land-use"], z = 50, [unit = u"m", description = "top of the surface layer"], z₀ = 0.04, [unit = u"m", description = "roughness lenght"], u_star = 0.44, [unit = u"m/s", description = "friction velocity"], diff --git a/src/wesley1989.jl b/src/wesley1989.jl index 62f64e9e..11fc7468 100644 --- a/src/wesley1989.jl +++ b/src/wesley1989.jl @@ -84,10 +84,17 @@ const Nh3Data = GasData(0.97, 2.e4, 0) # Changed according to Walmsley (1996) const PanData = GasData(2.6, 3.6, 0.1) # Peroxyacetyl nitrate const Hno2Data = GasData(1.6, 1.e5, 0.1) # Nitrous acid +# Obtain values from matrix using symbolic parameter iSeason and iLandUse +function obtain_value(iSeason, iLandUse, matrix) + index=(iLandUse-1)*5+iSeason + interpolate_r_i = DataInterpolations.LinearInterpolation(vec(matrix),1:55) + interpolate_r_i(index) +end + # Calculate bulk canopy stomatal resistance [s m-1] based on Wesely (1989) equation 3 when given the solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the season index (iSeason), the land use index (iLandUse), and whether there is currently rain or dew. -function r_s(G, Ts, iSeason::Int, iLandUse::Int, rainOrDew::Bool) +function r_s(G, Ts, iSeason, iLandUse, rainOrDew::Bool) rs = 0.0 - rs = IfElse.ifelse((Ts >= 39.9) & (Ts <= 0.1), inf, r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) * (400.0 * 1.0 / (Ts * (40.0 - Ts)))) + rs = IfElse.ifelse((Ts >= 39.9) & (Ts <= 0.1), inf, obtain_value(iSeason, iLandUse,r_i) * (1 + (200.0 * 1.0 / (G + 1.0))^2) * (400.0 * 1.0 / (Ts * (40.0 - Ts)))) # Adjust for dew and rain (from "Effects of dew and rain" section). if rainOrDew rs *= 3 @@ -111,7 +118,7 @@ function r_smx(r_s, Dh2oPerDx::T, r_mx::T) where {T<:AbstractFloat} end # Calculate the resistance of the outer surfaces in the upper canopy (leaf cuticular resistance in healthy vegetation) based on Wesely (1989) equations 7 and 10-14 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), the land use index (iLandUse), whether there is currently rain or dew, and whether the chemical of interest is either SO2 or O3. -function r_lux(Hstar::T, fo::T, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) where {T<:AbstractFloat} +function r_lux(Hstar::T, fo::T, iSeason, iLandUse, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) where {T<:AbstractFloat} rlux = 0.0 if dew && (iSeason != 4) # Dew doesn't have any effect in the winter if isSO2 @@ -122,10 +129,10 @@ function r_lux(Hstar::T, fo::T, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bo end elseif isO3 # equation 11 - rlux = 1.0 / (1.0 / 3000.0 + 1.0 / (3 * r_lu[iSeason, iLandUse])) + rlux = 1.0 / (1.0 / 3000.0 + 1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu))) else - rluO = 1.0 / (1.0 / 3000.0 + 1.0 / (3 * r_lu[iSeason, iLandUse])) #equation 11 - rlux = 1.0 / (1.0 / (3 * r_lu[iSeason, iLandUse] / (1.0e-5 * Hstar + fo)) + rluO = 1.0 / (1.0 / 3000.0 + 1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu))) #equation 11 + rlux = 1.0 / (1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu) / (1.0e-5 * Hstar + fo)) + 1.0e-7 * Hstar + fo / rluO) # equation 14, modified to match Walmsley eq. 5g end elseif rain && (iSeason != 4) @@ -134,31 +141,31 @@ function r_lux(Hstar::T, fo::T, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bo rlux = 50 #equation 13 and a half else # equation 12 - rlux = 1.0 / (1.0 / 5000.0 + 1.0 / (3 * r_lu[iSeason, iLandUse])) + rlux = 1.0 / (1.0 / 5000.0 + 1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu))) end elseif isO3 # equation 13 - rlux = 1.0 / (1.0 / 1000.0 + 1.0 / (3 * r_lu[iSeason, iLandUse])) + rlux = 1.0 / (1.0 / 1000.0 + 1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu))) else - rluO = 1.0 / (1.0 / 1000.0 + 1.0 / (3 * r_lu[iSeason, iLandUse])) # equation 13 - rlux = 1.0 / (1.0 / (3 * r_lu[iSeason, iLandUse] / (1.e-5 * Hstar + fo)) + 1.0e-7 * Hstar + fo / rluO) # equation 14, modified to match Walmsley eq. 5g + rluO = 1.0 / (1.0 / 1000.0 + 1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu))) # equation 13 + rlux = 1.0 / (1.0 / (3 * obtain_value(iSeason, iLandUse,r_lu) / (1.e-5 * Hstar + fo)) + 1.0e-7 * Hstar + fo / rluO) # equation 14, modified to match Walmsley eq. 5g end else - rlux = r_lu[iSeason, iLandUse] / (1.0e-5 * Hstar + fo) + rlux = obtain_value(iSeason, iLandUse,r_lu) / (1.0e-5 * Hstar + fo) end return rlux end # Calculate the resistance of the exposed surfaces in the lower portions of structures (canopies, buildings) above the ground based on Wesely (1989) equation 8 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), and the land use index (iLandUse). -function r_clx(Hstar::T, fo::T, iSeason::Int, iLandUse::Int) where {T<:AbstractFloat} - return 1.0 / (Hstar / (1.0e5 * r_clS[iSeason, iLandUse]) + - fo / r_clO[iSeason, iLandUse]) +function r_clx(Hstar::T, fo::T, iSeason, iLandUse) where {T<:AbstractFloat} + return 1.0 / (Hstar / (1.0e5 * obtain_value(iSeason, iLandUse,r_clS)) + + fo / obtain_value(iSeason, iLandUse,r_clO)) end # Calculate the resistance to uptake at the 'ground' surface based on Wesely (1989) equation 9 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), and the land use index (iLandUse). -function r_gsx(Hstar::T, fo::T, iSeason::Int, iLandUse::Int) where {T<:AbstractFloat} - return 1.0 / (Hstar / (1.0e5 * r_gsS[iSeason, iLandUse]) + - fo / r_gsO[iSeason, iLandUse]) +function r_gsx(Hstar::T, fo::T, iSeason, iLandUse) where {T<:AbstractFloat} + return 1.0 / (Hstar / (1.0e5 * obtain_value(iSeason, iLandUse,r_gsS)) + + fo / obtain_value(iSeason, iLandUse,r_gsO)) end # Calculates surface resistance to dry depostion [s m-1] based on Wesely (1989) equation 2 when given information on the chemical of interest (gasData), solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the slope of the local terrain (Θ [radians]), the season index (iSeason), the land use index (iLandUse), whether there is currently rain or dew, and whether the chemical of interest is either SO2 (isSO2) or O3 (isO3). @@ -174,17 +181,17 @@ function WesleySurfaceResistance(gasData::GasData, G, Ts, θ, rclx = 0.0 rgsx = 0.0 if isSO2 - rclx = r_clS[iSeason, iLandUse] - rgsx = r_gsS[iSeason, iLandUse] + rclx = obtain_value(iSeason, iLandUse,r_clS) + rgsx = obtain_value(iSeason, iLandUse,r_gsS) elseif isO3 - rclx = r_clO[iSeason, iLandUse] - rgsx = r_gsO[iSeason, iLandUse] + rclx = obtain_value(iSeason, iLandUse,r_clO) + rgsx = obtain_value(iSeason, iLandUse,r_gsO) else rclx = r_clx(gasData.Hstar, gasData.Fo, iSeason, iLandUse) rgsx = r_gsx(gasData.Hstar, gasData.Fo, iSeason, iLandUse) end - rac = r_ac[iSeason, iLandUse] + rac = obtain_value(iSeason, iLandUse,r_ac) # Correction for cold temperatures from page 4 column 1. correction = IfElse.ifelse((Ts < 0.0), 1000.0 * exp(-Ts - 4), 0) # [s m-1] #mark diff --git a/test/drydep_test.jl b/test/drydep_test.jl index 517e3df4..935d319e 100644 --- a/test/drydep_test.jl +++ b/test/drydep_test.jl @@ -13,7 +13,7 @@ begin @parameters ρParticle [unit = u"kg*m^-3"] @parameters ρA [unit = u"kg*m^-3"] @parameters G [unit = u"W*m^-2"] - @parameters θ + @parameters θ, iLandUse, iSeason end @testset "mfp" begin @@ -28,7 +28,7 @@ end #@test unit(dH2O(300u"K"))==u"m^2/s" @test ModelingToolkit.get_unit(DryDepParticle(z, z₀, u_star, L, Dp, T, P, ρParticle, ρA, 1, 1)) == u"m/s" #@test unit(DryDepParticle(0.4u"m",0.3u"m",1u"m/s", 1u"m", 1e-6u"m", 300u"K", 10300u"Pa", 1u"kg*m^-3",0.001u"kg*m^-3",1,1)) == u"m/s" - @test ModelingToolkit.get_unit(DryDepGas(z, z₀, u_star, L, ρA, AtmosphericDeposition.So2Data, G, T, θ, 1, 1, false, false, true, false)) == u"m/s" + @test ModelingToolkit.get_unit(DryDepGas(z, z₀, u_star, L, ρA, AtmosphericDeposition.So2Data, G, T, θ, iSeason, iLandUse, false, false, true, false)) == u"m/s" #@test unit(DryDepGas(0.4u"m",0.3u"m",1u"m/s", 1u"m", 0.001u"kg*m^-3", So2Data, 800u"W*m^-2", 300u"K", 0, 1, 1, false, false, true, false)) == u"m/s" end @@ -73,7 +73,7 @@ end @testset "DryDepGas" begin vd_true = 0.03 # m/s - @test (substitute(DryDepGas(z, z₀, u_star, L, ρA, AtmosphericDeposition.No2Data, G, T, 0, 1, 10, false, false, false, false), Dict(z => 50, z₀ => 0.04, u_star => 0.44, L => 0, T => 298, ρA => 1.2, G => 300, defaults...)) - vd_true) / vd_true < 0.33 + @test (substitute(DryDepGas(z, z₀, u_star, L, ρA, AtmosphericDeposition.No2Data, G, T, 0, iSeason, iLandUse, false, false, false, false), Dict(z => 50, z₀ => 0.04, u_star => 0.44, L => 0, T => 298, ρA => 1.2, G => 300, iSeason => 1, iLandUse => 10, defaults...)) - vd_true) / vd_true < 0.33 end @testset "DryDepParticle" begin From 69aa7a91947d6a8e42b074489a36091a67674dff Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Fri, 6 Sep 2024 18:06:46 -0500 Subject: [PATCH 04/15] Add connections with geosfp to calculate L and qrain --- ext/EarthSciDataExt.jl | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/ext/EarthSciDataExt.jl b/ext/EarthSciDataExt.jl index 08babae6..40632acc 100644 --- a/ext/EarthSciDataExt.jl +++ b/ext/EarthSciDataExt.jl @@ -10,13 +10,18 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler, MW_air = 29, [unit = u"g/mol", description="dry air molar mass"], kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"], R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"], + vK = 0.4, [description = "von Karman's constant"], + Cp = 1000, [unit = u"W*s/kg/K", description="specific heat at constant pressure"], + gg = 9.81, [unit = u"m*s^-2", description="gravitational acceleration"], ) # ρA in DrydepositionG(t) are in units of "kg/m3". # ρ = P*M/(R*T)= Pa*(g/mol)/(m3*Pa/mol/K*K) = g/m3 # Overall, ρA = P*M/(R*T)*kgperg - d = param_to_var(d, :T, :z, :z₀, :u_star, :G, :ρA) + # Monin-Obhukov length = -Air density * Cp * T(surface air) * Ustar^3/（vK * g * Sensible Heat flux） + + d = param_to_var(d, :T, :z, :z₀, :u_star, :G, :ρA, :L) ConnectorSystem([ d.T ~ g.I3₊T, d.z ~ g.A1₊PBLH, @@ -24,6 +29,7 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler, d.u_star ~ g.A1₊USTAR, d.G ~ g. A1₊SWGDN, d.ρA ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air, + d.L ~ -g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air * Cp * g.A1₊TS * (g.A1₊USTAR)^3/(vK*gg*g.A1₊HFLUX), ], d, g) end @@ -35,16 +41,21 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g MW_air = 29, [unit = u"g/mol", description="dry air molar mass"], kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"], R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"], + Vdr = 5.0, [unit = u"m/s", description="droplet velocity"], ) # ρ_air in Wetdeposition(t) are in units of "kg/m3". # ρ = P*M/(R*T)= Pa*(g/mol)/(m3*Pa/mol/K*K) = g/m3 # Overall, ρ_air = P*M/(R*T)*kgperg - d = param_to_var(d, :cloudFrac, :ρ_air) + # From EMEP algorithm: P = QRAIN * Vdr * ρgas => QRAIN = P / Vdr / ρgas + # kg*m-2*s-1/(m/s)/(kg/m3) + + d = param_to_var(d, :cloudFrac, :ρ_air, :qrain) ConnectorSystem([ d.cloudFrac ~ g.A3cld₊CLOUD, d.ρ_air ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air, + d.qrain ~ (g.A3mstE₊PFLCU + g.A3mstE₊PFLLSAN) / Vdr / (g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air), ], d, g) end From e8d87bfcb2ea5812ae05b6a4d37a5d03ce06fb3c Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Thu, 26 Sep 2024 15:54:06 -0500 Subject: [PATCH 05/15] Update Deposition to be compatible with MTKv9 --- Project.toml | 13 +++++-------- ext/EarthSciDataExt.jl | 20 ++++++++------------ ext/GasChemExt.jl | 2 +- src/AtmosphericDeposition.jl | 6 ++++-- src/dry_deposition.jl | 26 ++++++++++++++------------ src/wesley1989.jl | 5 +++-- src/wet_deposition.jl | 16 ++++++++-------- test/connector_test.jl | 32 +++++++++++++++----------------- test/drydep_test.jl | 22 +++++++++------------- test/wetdep_test.jl | 2 +- 10 files changed, 68 insertions(+), 76 deletions(-) diff --git a/Project.toml b/Project.toml index 977787e1..bcf6e927 100644 --- a/Project.toml +++ b/Project.toml @@ -6,15 +6,13 @@ version = "0.2.1" [deps] DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821" EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0" -IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [weakdeps] EarthSciData = "a293c155-435f-439d-9c11-a083b6b47337" @@ -26,13 +24,12 @@ GasChemExt = "GasChem" [compat] DifferentialEquations = "7" -EarthSciMLBase = "0.15" -GasChem = "0.6" -IfElse = "0.1" -ModelingToolkit = "8" +EarthSciData = "0.9" +EarthSciMLBase = "0.16" +GasChem = "0.7" +ModelingToolkit = "9" SafeTestsets = "0.1" StaticArrays = "1" -Unitful = "1" julia = "1" [extras] diff --git a/ext/EarthSciDataExt.jl b/ext/EarthSciDataExt.jl index 40632acc..e7bc628c 100644 --- a/ext/EarthSciDataExt.jl +++ b/ext/EarthSciDataExt.jl @@ -1,14 +1,12 @@ module EarthSciDataExt -using AtmosphericDeposition, EarthSciData, EarthSciMLBase, Unitful, ModelingToolkit +using AtmosphericDeposition, EarthSciData, EarthSciMLBase, DynamicQuantities, ModelingToolkit function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler, g::EarthSciData.GEOSFPCoupler) d, g = d.sys, g.sys @constants( - PaPerhPa = 100, [unit = u"Pa/hPa", description="Conversion factor from hPa to Pa"], - MW_air = 29, [unit = u"g/mol", description="dry air molar mass"], - kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"], + MW_air = 0.029, [unit = u"kg/mol", description="dry air molar mass"], R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"], vK = 0.4, [description = "von Karman's constant"], Cp = 1000, [unit = u"W*s/kg/K", description="specific heat at constant pressure"], @@ -27,9 +25,9 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler, d.z ~ g.A1₊PBLH, d.z₀ ~ g.A1₊Z0M, d.u_star ~ g.A1₊USTAR, - d.G ~ g. A1₊SWGDN, - d.ρA ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air, - d.L ~ -g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air * Cp * g.A1₊TS * (g.A1₊USTAR)^3/(vK*gg*g.A1₊HFLUX), + d.G ~ g.A1₊SWGDN, + d.ρA ~ g.P/(g.I3₊T*R)*MW_air, + d.L ~ -g.P/(g.I3₊T*R)*MW_air * Cp * g.A1₊TS * (g.A1₊USTAR)^3/(vK*gg*g.A1₊HFLUX), ], d, g) end @@ -37,9 +35,7 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g d, g = d.sys, g.sys @constants( - PaPerhPa = 100, [unit = u"Pa/hPa", description="Conversion factor from hPa to Pa"], - MW_air = 29, [unit = u"g/mol", description="dry air molar mass"], - kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"], + MW_air = 0.029, [unit = u"kg/mol", description="dry air molar mass"], R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"], Vdr = 5.0, [unit = u"m/s", description="droplet velocity"], ) @@ -54,8 +50,8 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g d = param_to_var(d, :cloudFrac, :ρ_air, :qrain) ConnectorSystem([ d.cloudFrac ~ g.A3cld₊CLOUD, - d.ρ_air ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air, - d.qrain ~ (g.A3mstE₊PFLCU + g.A3mstE₊PFLLSAN) / Vdr / (g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air), + d.ρ_air ~ g.P/(g.I3₊T*R)*MW_air, + d.qrain ~ (g.A3mstE₊PFLCU + g.A3mstE₊PFLLSAN) / Vdr / (g.P/(g.I3₊T*R)*MW_air), ], d, g) end diff --git a/ext/GasChemExt.jl b/ext/GasChemExt.jl index f7420fca..be976bf6 100644 --- a/ext/GasChemExt.jl +++ b/ext/GasChemExt.jl @@ -1,6 +1,6 @@ module GasChemExt -using AtmosphericDeposition, GasChem, EarthSciMLBase, Unitful, ModelingToolkit +using AtmosphericDeposition, GasChem, EarthSciMLBase, DynamicQuantities, ModelingToolkit function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDeposition.DrydepositionGCoupler) c, d = c.sys, d.sys diff --git a/src/AtmosphericDeposition.jl b/src/AtmosphericDeposition.jl index a7ccfe4f..98100af7 100644 --- a/src/AtmosphericDeposition.jl +++ b/src/AtmosphericDeposition.jl @@ -1,11 +1,13 @@ module AtmosphericDeposition -using Unitful +using DynamicQuantities using StaticArrays using ModelingToolkit using EarthSciMLBase -using IfElse using DataInterpolations +using ModelingToolkit:t + +@register_unit ppb 1 include("wesley1989.jl") include("dry_deposition.jl") diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 3d1234ff..1f18c4ad 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -14,12 +14,14 @@ Based on Seinfeld and Pandis (2006) [Seinfeld, J.H. and Pandis, S.N. (2006) Atmo equation 19.13 & 19.14. """ function ra(z, z₀, u_star, L) - ζ = IfElse.ifelse((L / unit_m == 0), 0, z / L) - ζ₀ = IfElse.ifelse((L / unit_m == 0), 0, z₀ / L) - rₐ_1 = IfElse.ifelse((0 < ζ) & (ζ < 1), 1 / (κ * u_star) * (log(z / z₀) + 4.7 * (ζ - ζ₀)), 1 / (κ * u_star) * log(z / z₀)) + ζ = ifelse((L / unit_m == 0), 0, z / L) + ζ₀ = ifelse((L / unit_m == 0), 0, z₀ / L) + rₐ_1 = ifelse((0 < ζ), 1 / (κ * u_star) * (log(z / z₀) + 4.7 * (ζ - ζ₀)), 1 / (κ * u_star) * log(z / z₀)) + rₐ_1 = ifelse((ζ < 1), 1 / (κ * u_star) * (log(z / z₀) + 4.7 * (ζ - ζ₀)), 1 / (κ * u_star) * log(z / z₀)) η₀ = (1 - 15 * ζ₀)^1 / 4 η = (1 - 15 * ζ)^1 / 4 - rₐ = IfElse.ifelse((-1 < ζ) & (ζ < 0), 1 / (κ * u_star) * [log(z / z₀) + log(((η₀^2 + 1) * (η₀ + 1)^2) / ((η^2 + 1) * (η + 1)^2)) + 2 * (atan(η) - atan(η₀))][1], rₐ_1) + rₐ = ifelse((-1 < ζ), 1 / (κ * u_star) * [log(z / z₀) + log(((η₀^2 + 1) * (η₀ + 1)^2) / ((η^2 + 1) * (η + 1)^2)) + 2 * (atan(η) - atan(η₀))][1], rₐ_1) + rₐ = ifelse((ζ < 0), 1 / (κ * u_star) * [log(z / z₀) + log(((η₀^2 + 1) * (η₀ + 1)^2) / ((η^2 + 1) * (η + 1)^2)) + 2 * (atan(η) - atan(η₀))][1], rₐ_1) return rₐ end @@ -58,7 +60,7 @@ particle where Dp is particle diameter [m], ρₚ is particle density [kg/m3], C From equation 9.42 in Seinfeld and Pandis (2006). """ function vs(Dₚ, ρₚ, Cc, μ) - IfElse.ifelse((Dₚ > 20.e-6 * unit_m), 99999999 * unit_v, Dₚ^2 * ρₚ * g * Cc / (18 * μ)) + ifelse((Dₚ > 20.e-6 * unit_m), 99999999 * unit_v, Dₚ^2 * ρₚ * g * Cc / (18 * μ)) # Particle diameter Dₚ greater than 20um; Stokes settling no longer applies. end @@ -233,7 +235,7 @@ Build Drydeposition model (gas) d = DrydepositionG(t) ``` """ -function DrydepositionG(t; name=:DrydepositionG) +function DrydepositionG(; name=:DrydepositionG) rain = false dew = false params = @parameters( @@ -252,12 +254,12 @@ function DrydepositionG(t; name=:DrydepositionG) D = Differential(t) vars = @variables( - SO2(t), [unit = u"nmol/mol"], - O3(t),[unit = u"nmol/mol"], - NO2(t), [unit = u"nmol/mol"], - NO(t), [unit = u"nmol/mol"], - H2O2(t), [unit = u"nmol/mol"], - CH2O(t), [unit = u"nmol/mol"], + SO2(t), [unit = u"ppb"], + O3(t),[unit = u"ppb"], + NO2(t), [unit = u"ppb"], + NO(t), [unit = u"ppb"], + H2O2(t), [unit = u"ppb"], + CH2O(t), [unit = u"ppb"], ) eqs = [ diff --git a/src/wesley1989.jl b/src/wesley1989.jl index 11fc7468..e2661a36 100644 --- a/src/wesley1989.jl +++ b/src/wesley1989.jl @@ -94,7 +94,8 @@ end # Calculate bulk canopy stomatal resistance [s m-1] based on Wesely (1989) equation 3 when given the solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the season index (iSeason), the land use index (iLandUse), and whether there is currently rain or dew. function r_s(G, Ts, iSeason, iLandUse, rainOrDew::Bool) rs = 0.0 - rs = IfElse.ifelse((Ts >= 39.9) & (Ts <= 0.1), inf, obtain_value(iSeason, iLandUse,r_i) * (1 + (200.0 * 1.0 / (G + 1.0))^2) * (400.0 * 1.0 / (Ts * (40.0 - Ts)))) + rs = ifelse((Ts >= 39.9), inf, obtain_value(iSeason, iLandUse,r_i) * (1 + (200.0 * 1.0 / (G + 1.0))^2) * (400.0 * 1.0 / (Ts * (40.0 - Ts)))) + rs = ifelse((Ts <= 0.1), inf, obtain_value(iSeason, iLandUse,r_i) * (1 + (200.0 * 1.0 / (G + 1.0))^2) * (400.0 * 1.0 / (Ts * (40.0 - Ts)))) # Adjust for dew and rain (from "Effects of dew and rain" section). if rainOrDew rs *= 3 @@ -194,7 +195,7 @@ function WesleySurfaceResistance(gasData::GasData, G, Ts, θ, rac = obtain_value(iSeason, iLandUse,r_ac) # Correction for cold temperatures from page 4 column 1. - correction = IfElse.ifelse((Ts < 0.0), 1000.0 * exp(-Ts - 4), 0) # [s m-1] #mark + correction = ifelse((Ts < 0.0), 1000.0 * exp(-Ts - 4), 0) # [s m-1] #mark rlux += correction rclx += correction rgsx += correction diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index d143195c..a19a8859 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -54,7 +54,7 @@ Build Wetdeposition model wd = Wetdeposition(t) ``` """ -function Wetdeposition(t; name=:Wetdeposition) +function Wetdeposition(; name=:Wetdeposition) params = @parameters( cloudFrac = 0.5, [description = "fraction of grid cell covered by clouds"], qrain = 0.5, [description = "rain mixing ratio"], @@ -65,13 +65,13 @@ function Wetdeposition(t; name=:Wetdeposition) D = Differential(t) vars = @variables( - SO2(t), [unit = u"nmol/mol"], - O3(t), [unit = u"nmol/mol"], - NO2(t), [unit = u"nmol/mol"], - CH4(t), [unit = u"nmol/mol"], - CO(t), [unit = u"nmol/mol"], - DMS(t), [unit = u"nmol/mol"], - ISOP(t), [unit = u"nmol/mol"], + SO2(t), [unit = u"ppb"], + O3(t), [unit = u"ppb"], + NO2(t), [unit = u"ppb"], + CH4(t), [unit = u"ppb"], + CO(t), [unit = u"ppb"], + DMS(t), [unit = u"ppb"], + ISOP(t), [unit = u"ppb"], ) eqs = [ diff --git a/test/connector_test.jl b/test/connector_test.jl index 7ee7f918..264ce166 100644 --- a/test/connector_test.jl +++ b/test/connector_test.jl @@ -1,13 +1,13 @@ using AtmosphericDeposition -using Test, Unitful, ModelingToolkit, GasChem, Dates, EarthSciMLBase, EarthSciData +using Test, DynamicQuantities, ModelingToolkit, GasChem, Dates, EarthSciMLBase, EarthSciData +using ModelingToolkit:t @testset "GasChemExt" begin start = Dates.datetime2unix(Dates.DateTime(2016, 5, 1)) - @parameters t [unit = u"s"] - composed_ode = couple(SuperFast(t), FastJX(t), DrydepositionG(t), Wetdeposition(t)) - tspan = (start, start+3600*24*3) - sys = structural_simplify(get_mtk(composed_ode)) - @test length(states(sys)) ≈ 18 + composed_ode = couple(SuperFast(), FastJX(), DrydepositionG(), Wetdeposition()) + combined_mtk = convert(ODESystem, composed_ode) + sys = structural_simplify(combined_mtk) + @test length(unknowns(sys)) ≈ 18 eqs = string(equations(sys)) wanteqs = ["Differential(t)(SuperFast₊O3(t)) ~ SuperFast₊DrydepositionG_ddt_O3ˍt(t) + SuperFast₊Wetdeposition_ddt_O3ˍt(t)"] @@ -15,23 +15,21 @@ using Test, Unitful, ModelingToolkit, GasChem, Dates, EarthSciMLBase, EarthSciDa end @testset "EarthSciDataExt" begin - @parameters t [unit = u"s"] - - @parameters lat = 40 - @parameters lon = -97 + @parameters lat = deg2rad(40.0f0) [unit=u"rad"] + @parameters lon = deg2rad(-97.0f0) [unit=u"rad"] @parameters lev = 1 - geosfp = GEOSFP("4x5", t) - - model = couple(SuperFast(t), FastJX(t), geosfp, Wetdeposition(t), DrydepositionG(t)) + geosfp = GEOSFP("4x5") + + model = couple(SuperFast(), FastJX(), geosfp, Wetdeposition(), DrydepositionG()) - sys = structural_simplify(get_mtk(model)) - @test length(states(sys)) ≈ 18 + sys = structural_simplify(convert(ODESystem, model)) + @test length(unknowns(sys)) ≈ 18 - eqs = string(equations(get_mtk(model))) + eqs = string(equations(convert(ODESystem, model))) wanteq = "DrydepositionG₊G(t) ~ GEOSFP₊A1₊SWGDN(t)" @test contains(eqs, wanteq) wanteq = "Wetdeposition₊cloudFrac(t) ~ GEOSFP₊A3cld₊CLOUD(t)" @test contains(eqs, wanteq) - wanted = "Wetdeposition₊ρA(t) ~ GEOSFP₊P * PaPerhPa/(GEOSFP₊I3₊T*R)*kgperg*MW_air" + wanted = "Wetdeposition₊ρA(t) ~ GEOSFP₊P/(GEOSFP₊I3₊T*R)*kgperg*MW_air" @test contains(eqs, wanteq) end \ No newline at end of file diff --git a/test/drydep_test.jl b/test/drydep_test.jl index 935d319e..2594be5c 100644 --- a/test/drydep_test.jl +++ b/test/drydep_test.jl @@ -1,5 +1,5 @@ using AtmosphericDeposition -using Test, Unitful, ModelingToolkit +using Test, DynamicQuantities, ModelingToolkit begin @parameters T [unit = u"K"] @@ -17,19 +17,15 @@ begin end @testset "mfp" begin - @test substitute(mfp(T, P, μ), Dict(T => 298, P => 101300, μ => 1.8e-5, defaults...)) ≈ 6.512893276888993e-8 + @test substitute(mfp(T, P, μ), Dict(T => 298, P => 101300, μ => 1.8e-5, AtmosphericDeposition.defaults...)) ≈ 6.512893276888993e-8 @test ModelingToolkit.get_unit(mfp(T, P, μ)) == u"m" - #@test unit(ModelingToolkit.get_unit(mfp(T,P,μ))*1 - 1u"m") == u"m" - #@test mfp(298u"K",101300u"Pa",1.8e-5u"kg/m/s") - 6.51e-8u"m" ≈ 0u"m" atol=1e-8u"m" end + @testset "unit" begin @test ModelingToolkit.get_unit(dH2O(T)) == u"m^2/s" - #@test unit(dH2O(300u"K"))==u"m^2/s" @test ModelingToolkit.get_unit(DryDepParticle(z, z₀, u_star, L, Dp, T, P, ρParticle, ρA, 1, 1)) == u"m/s" - #@test unit(DryDepParticle(0.4u"m",0.3u"m",1u"m/s", 1u"m", 1e-6u"m", 300u"K", 10300u"Pa", 1u"kg*m^-3",0.001u"kg*m^-3",1,1)) == u"m/s" @test ModelingToolkit.get_unit(DryDepGas(z, z₀, u_star, L, ρA, AtmosphericDeposition.So2Data, G, T, θ, iSeason, iLandUse, false, false, true, false)) == u"m/s" - #@test unit(DryDepGas(0.4u"m",0.3u"m",1u"m/s", 1u"m", 0.001u"kg*m^-3", So2Data, 800u"W*m^-2", 300u"K", 0, 1, 1, false, false, true, false)) == u"m/s" end @testset "viscosity" begin @@ -37,7 +33,7 @@ end μ_list = [1.725, 1.846, 1.962, 2.075, 2.181, 2.286] .* 10^-5 μ_test = [] for i in 1:6 - push!(μ_test, substitute(mu(T), Dict(T => T_[i], defaults...))) + push!(μ_test, substitute(mu(T), Dict(T => T_[i], AtmosphericDeposition.defaults...))) end for i in 1:6 @test (μ_test[i] - μ_list[i]) / μ_list[i] < 0.01 @@ -49,7 +45,7 @@ end Cc_list = [216, 108, 43.6, 22.2, 11.4, 4.95, 2.85, 1.865, 1.326, 1.164, 1.082, 1.032, 1.016, 1.008, 1.003, 1.0016] Cc_test = [] for i in 1:16 - push!(Cc_test, substitute(cc(Dp, T, P, μ), Dict(Dp => Dp_[i] * 1e-6, T => 298, P => 101325, μ => substitute(mu(T), Dict(T => 298, defaults...)), defaults...))) + push!(Cc_test, substitute(cc(Dp, T, P, μ), Dict(Dp => Dp_[i] * 1e-6, T => 298, P => 101325, μ => substitute(mu(T), Dict(T => 298, AtmosphericDeposition.defaults...)), AtmosphericDeposition.defaults...))) end for i in 1:16 @test (Cc_test[i] - Cc_list[i]) / Cc_list[i] < 0.03 @@ -63,8 +59,8 @@ end Vs_test = [] Cc_list = [] for i in 1:4 - push!(Cc_list, substitute(cc(Dp, T, P, μ), Dict(Dp => Dp_[i] * 1e-6, T => 298, P => 101325, μ => substitute(mu(T), Dict(T => 298, defaults...)), defaults...))) - push!(Vs_test, substitute(vs(Dp, ρParticle, Cc, μ), Dict(Dp => Dp_[i] * 1e-6, ρParticle => 1000, Cc => Cc_list[i], μ => 1.836522217711828e-5, defaults...))) + push!(Cc_list, substitute(cc(Dp, T, P, μ), Dict(Dp => Dp_[i] * 1e-6, T => 298, P => 101325, μ => substitute(mu(T), Dict(T => 298, AtmosphericDeposition.defaults...)), AtmosphericDeposition.defaults...))) + push!(Vs_test, substitute(vs(Dp, ρParticle, Cc, μ), Dict(Dp => Dp_[i] * 1e-6, ρParticle => 1000, Cc => Cc_list[i], μ => 1.836522217711828e-5, AtmosphericDeposition.defaults...))) end for i in 1:4 @test (Vs_test[i] - Vs_list[i]) / Vs_list[i] < 1 @@ -73,7 +69,7 @@ end @testset "DryDepGas" begin vd_true = 0.03 # m/s - @test (substitute(DryDepGas(z, z₀, u_star, L, ρA, AtmosphericDeposition.No2Data, G, T, 0, iSeason, iLandUse, false, false, false, false), Dict(z => 50, z₀ => 0.04, u_star => 0.44, L => 0, T => 298, ρA => 1.2, G => 300, iSeason => 1, iLandUse => 10, defaults...)) - vd_true) / vd_true < 0.33 + @test (substitute(DryDepGas(z, z₀, u_star, L, ρA, AtmosphericDeposition.No2Data, G, T, 0, iSeason, iLandUse, false, false, false, false), Dict(z => 50, z₀ => 0.04, u_star => 0.44, L => 0, T => 298, ρA => 1.2, G => 300, iSeason => 1, iLandUse => 10, AtmosphericDeposition.defaults...)) - vd_true) / vd_true < 0.33 end @testset "DryDepParticle" begin @@ -81,7 +77,7 @@ end vd_true = [0.5, 0.012, 0.02, 0.6] ./ 100 # [m/s] vd_list = [] for i in 1:4 - push!(vd_list, substitute(DryDepParticle(z, z₀, u_star, L, Dp, T, P, ρParticle, ρA, 1, 4), Dict(z => 20, z₀ => 0.02, u_star => 0.44, L => 0, T => 298, P => 101325, ρA => 1.2, ρParticle => 1000, Dp => Dp_[i], defaults...))) + push!(vd_list, substitute(DryDepParticle(z, z₀, u_star, L, Dp, T, P, ρParticle, ρA, 1, 4), Dict(z => 20, z₀ => 0.02, u_star => 0.44, L => 0, T => 298, P => 101325, ρA => 1.2, ρParticle => 1000, Dp => Dp_[i], AtmosphericDeposition.defaults...))) end for i in 1:4 @test vd_list[i] - vd_true[i] < 0.015 diff --git a/test/wetdep_test.jl b/test/wetdep_test.jl index 8157aeb6..aece7d7c 100644 --- a/test/wetdep_test.jl +++ b/test/wetdep_test.jl @@ -1,5 +1,5 @@ using AtmosphericDeposition -using Test, Unitful, ModelingToolkit +using Test, DynamicQuantities, ModelingToolkit @parameters cloudFrac = 0.5 @parameters qrain = 0.5 From 98193f8cff557f9485fa0c237f87145b9a1d2ccb Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Mon, 30 Sep 2024 16:47:23 -0500 Subject: [PATCH 06/15] Update documentation andinclude species HNO3, H2O2, CH2O --- docs/Project.toml | 2 +- docs/src/dry_deposition.md | 26 +++++++++++++------------- docs/src/emep.md | 16 ++++++++-------- ext/EarthSciDataExt.jl | 4 ++-- ext/GasChemExt.jl | 2 ++ src/dry_deposition.jl | 2 ++ src/wet_deposition.jl | 6 ++++++ test/connector_test.jl | 7 +++++-- 8 files changed, 39 insertions(+), 26 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 9bfbf718..c75f724b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,5 +4,5 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" -Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" \ No newline at end of file diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md index fc2bcfec..15fbc70a 100644 --- a/docs/src/dry_deposition.md +++ b/docs/src/dry_deposition.md @@ -3,7 +3,7 @@ This is an implementation of a box model used to calculate changes in gas species concentration due to dry deposition. ## Running the model -Here's an example of how concentration of different species, such as SO₂, O₃, NO₂, NO, H₂O₂ and CH₂O change due to dry deposition. +Here's an example of how concentration of different species, such as SO₂, O₃, NO₂, NO, H₂O₂, CH₂O and HNO₃change due to dry deposition. We can create an instance of the model in the following manner: ```julia @@ -11,16 +11,16 @@ using AtmosphericDeposition using ModelingToolkit using DifferentialEquations using EarthSciMLBase -using Unitful +using DynamicQuantities +using ModelingToolkit:t -@parameters t [unit = u"s", description="Time"] -model = DrydepositionG(t) +model = DrydepositionG() ``` Before running any simulations with the model we need to convert it into a system of differential equations. ```julia sys = structural_simplify(model) tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,5,2.34,0.15] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O +u0 = [2.0,10.0,5,5,2.34,0.15,10] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O, HNO₃ sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` @@ -29,14 +29,14 @@ using AtmosphericDeposition using ModelingToolkit using DifferentialEquations using EarthSciMLBase -using Unitful +using DynamicQuantities +using ModelingToolkit:t -@parameters t [unit = u"s", description="Time"] -model = DrydepositionG(t) +model = DrydepositionG() sys = structural_simplify(model) tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,5,2.34,0.15] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O +u0 = [2.0,10.0,5,5,2.34,0.15,10] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O, HNO₃ sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` @@ -49,16 +49,16 @@ plot(sol, xlabel="Time (second)", ylabel="concentration (ppb)", legend=:outerrig ## Parameters The parameters in the model are: ```julia -parameters(sys) # [z, z₀, u_star, L, ρA, G, T, θ] +parameters(sys) # [iLandUse, z, z₀, u_star, L, ρA, G, iSeason, T, θ] ``` -where ```z``` is the top of the surface layer [m], ```z₀``` is the roughness length [m], ```u_star``` is friction velocity [m/s], and ```L``` is Monin-Obukhov length [m], ```ρA``` is air density [kg/m3], ```T``` is surface air temperature [K], ```G``` is solar irradiation [W m-2], ```Θ``` is the slope of the local terrain [radians]. +where ```iSeason``` and ```iLandUse``` are the index for the season and landuse,```z``` is the top of the surface layer [m], ```z₀``` is the roughness length [m], ```u_star``` is friction velocity [m/s], and ```L``` is Monin-Obukhov length [m], ```ρA``` is air density [kg/m3], ```T``` is surface air temperature [K], ```G``` is solar irradiation [W m-2], ```Θ``` is the slope of the local terrain [radians]. Let's run some simulation with different value for parameter ```z```. ```@example 1 @unpack O3 = sys -p1 = [50,0.04,0.44,0,1.2,300,298,0] -p2 = [10,0.04,0.44,0,1.2,300,298,0] +p1 = [10, 50,0.04,0.44,0,1.2,300,1,298,0] +p2 = [10, 10,0.04,0.44,0,1.2,300,1,298,0] sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0) sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0) diff --git a/docs/src/emep.md b/docs/src/emep.md index 4a9205c2..463b3027 100644 --- a/docs/src/emep.md +++ b/docs/src/emep.md @@ -8,17 +8,17 @@ using AtmosphericDeposition using ModelingToolkit using DifferentialEquations using EarthSciMLBase -using Unitful +using DynamicQuantities +using ModelingToolkit:t -@parameters t [unit = u"s", description="Time"] -model = Wetdeposition(t) +model = Wetdeposition() ``` Before running any simulations with the model we need to convert it into a system of differential equations. ```julia sys = structural_simplify(model) tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,1400,275,50,0.15] # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP +u0 = [2.0,10.0,5,1400,275,50,0.15,2.34,10,0.15] # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP, H₂O₂, HNO₃, CH₂O prob = ODEProblem(sys, u0, tspan, []) sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` @@ -28,14 +28,14 @@ using AtmosphericDeposition using ModelingToolkit using DifferentialEquations using EarthSciMLBase -using Unitful +using DynamicQuantities +using ModelingToolkit:t -@parameters t [unit = u"s", description="Time"] -model = Wetdeposition(t) +model = Wetdeposition() sys = structural_simplify(model) tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,1400,275,50,0.15] # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP +u0 = [2.0,10.0,5,1400,275,50,0.15,2.34,10,0.15] # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP, H₂O₂, HNO₃, CH₂O prob = ODEProblem(sys, u0, tspan, []) sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` diff --git a/ext/EarthSciDataExt.jl b/ext/EarthSciDataExt.jl index e7bc628c..e5d19ef3 100644 --- a/ext/EarthSciDataExt.jl +++ b/ext/EarthSciDataExt.jl @@ -13,7 +13,7 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler, gg = 9.81, [unit = u"m*s^-2", description="gravitational acceleration"], ) - # ρA in DrydepositionG(t) are in units of "kg/m3". + # ρA in DrydepositionG() are in units of "kg/m3". # ρ = P*M/(R*T)= Pa*(g/mol)/(m3*Pa/mol/K*K) = g/m3 # Overall, ρA = P*M/(R*T)*kgperg @@ -40,7 +40,7 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g Vdr = 5.0, [unit = u"m/s", description="droplet velocity"], ) - # ρ_air in Wetdeposition(t) are in units of "kg/m3". + # ρ_air in Wetdeposition() are in units of "kg/m3". # ρ = P*M/(R*T)= Pa*(g/mol)/(m3*Pa/mol/K*K) = g/m3 # Overall, ρ_air = P*M/(R*T)*kgperg diff --git a/ext/GasChemExt.jl b/ext/GasChemExt.jl index be976bf6..6add2acf 100644 --- a/ext/GasChemExt.jl +++ b/ext/GasChemExt.jl @@ -26,6 +26,8 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos c.CO => d.CO, c.DMS => d.DMS, c.ISOP => d.ISOP, + c.H2O2 => d.H2O2, + c.CH2O => d.CH2O, )) end diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 1f18c4ad..7adff712 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -260,6 +260,7 @@ function DrydepositionG(; name=:DrydepositionG) NO(t), [unit = u"ppb"], H2O2(t), [unit = u"ppb"], CH2O(t), [unit = u"ppb"], + HNO3(t), [unit = u"ppb"], ) eqs = [ @@ -269,6 +270,7 @@ function DrydepositionG(; name=:DrydepositionG) D(NO) ~ -DryDepGas(z, z₀, u_star, L, ρA, NoData, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * NO D(H2O2) ~ -DryDepGas(z, z₀, u_star, L, ρA, H2o2Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * H2O2 D(CH2O) ~ -DryDepGas(z, z₀, u_star, L, ρA, HchoData, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * CH2O + D(HNO3) ~ -DryDepGas(z, z₀, u_star, L, ρA, Hno3Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * HNO3 ] ODESystem(eqs, t, vars, params; name=name, diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index a19a8859..83af024e 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -72,6 +72,9 @@ function Wetdeposition(; name=:Wetdeposition) CO(t), [unit = u"ppb"], DMS(t), [unit = u"ppb"], ISOP(t), [unit = u"ppb"], + H2O2(t), [unit = u"ppb"], + HNO3(t), [unit = u"ppb"], + CH2O(t), [unit = u"ppb"], ) eqs = [ @@ -82,6 +85,9 @@ function Wetdeposition(; name=:Wetdeposition) D(CO) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * CO D(DMS) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * DMS D(ISOP) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * ISOP + D(H2O2) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * H2O2 + D(HNO3) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * HNO3 + D(CH2O) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * CH2O ] ODESystem(eqs, t, vars, params; name=name, diff --git a/test/connector_test.jl b/test/connector_test.jl index 264ce166..20badd7b 100644 --- a/test/connector_test.jl +++ b/test/connector_test.jl @@ -7,7 +7,9 @@ using ModelingToolkit:t composed_ode = couple(SuperFast(), FastJX(), DrydepositionG(), Wetdeposition()) combined_mtk = convert(ODESystem, composed_ode) sys = structural_simplify(combined_mtk) - @test length(unknowns(sys)) ≈ 18 + print(unknowns(sys)) + @test length(unknowns(sys)) ≈ 20 + #TODO Change 20 to 18 after the latest GasChem package include species HNO3 eqs = string(equations(sys)) wanteqs = ["Differential(t)(SuperFast₊O3(t)) ~ SuperFast₊DrydepositionG_ddt_O3ˍt(t) + SuperFast₊Wetdeposition_ddt_O3ˍt(t)"] @@ -23,7 +25,8 @@ end model = couple(SuperFast(), FastJX(), geosfp, Wetdeposition(), DrydepositionG()) sys = structural_simplify(convert(ODESystem, model)) - @test length(unknowns(sys)) ≈ 18 + @test length(unknowns(sys)) ≈ 20 + #TODO Change 20 to 18 after the latest GasChem package include species HNO3 eqs = string(equations(convert(ODESystem, model))) wanteq = "DrydepositionG₊G(t) ~ GEOSFP₊A1₊SWGDN(t)" From aa9566d23c713d37ebbc45e69ea7f5fe8f954805 Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Mon, 30 Sep 2024 16:59:04 -0500 Subject: [PATCH 07/15] Update Documentation --- docs/src/dry_deposition.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md index 15fbc70a..26390908 100644 --- a/docs/src/dry_deposition.md +++ b/docs/src/dry_deposition.md @@ -51,14 +51,14 @@ The parameters in the model are: ```julia parameters(sys) # [iLandUse, z, z₀, u_star, L, ρA, G, iSeason, T, θ] ``` -where ```iSeason``` and ```iLandUse``` are the index for the season and landuse,```z``` is the top of the surface layer [m], ```z₀``` is the roughness length [m], ```u_star``` is friction velocity [m/s], and ```L``` is Monin-Obukhov length [m], ```ρA``` is air density [kg/m3], ```T``` is surface air temperature [K], ```G``` is solar irradiation [W m-2], ```Θ``` is the slope of the local terrain [radians]. +where ```iLandUse``` is the index for the landuse, ```iSeason``` is the index for the season, ```z``` is the top of the surface layer [m], ```z₀``` is the roughness length [m], ```u_star``` is friction velocity [m/s], and ```L``` is Monin-Obukhov length [m], ```ρA``` is air density [kg/m3], ```T``` is surface air temperature [K], ```G``` is solar irradiation [W m-2], ```Θ``` is the slope of the local terrain [radians]. Let's run some simulation with different value for parameter ```z```. ```@example 1 @unpack O3 = sys -p1 = [10, 50,0.04,0.44,0,1.2,300,1,298,0] -p2 = [10, 10,0.04,0.44,0,1.2,300,1,298,0] +p1 = [10,50,0.04,0.44,0,1.2,300,1,298,0] +p2 = [10,10,0.04,0.44,0,1.2,300,1,298,0] sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0) sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0) From 382c6b164c86cc777858fbee7f2b05c0f6c3bfde Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Mon, 30 Sep 2024 17:27:54 -0500 Subject: [PATCH 08/15] Correct the examples in documentation --- docs/src/dry_deposition.md | 6 +++--- docs/src/emep.md | 14 +++++++------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md index 26390908..8b449319 100644 --- a/docs/src/dry_deposition.md +++ b/docs/src/dry_deposition.md @@ -55,10 +55,10 @@ where ```iLandUse``` is the index for the landuse, ```iSeason``` is the index fo Let's run some simulation with different value for parameter ```z```. ```@example 1 -@unpack O3 = sys +@unpack O3,z = sys -p1 = [10,50,0.04,0.44,0,1.2,300,1,298,0] -p2 = [10,10,0.04,0.44,0,1.2,300,1,298,0] +p1 = [z=>50] +p2 = [z=>10] sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0) sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0) diff --git a/docs/src/emep.md b/docs/src/emep.md index 463b3027..865688e4 100644 --- a/docs/src/emep.md +++ b/docs/src/emep.md @@ -49,16 +49,16 @@ plot(sol, xlabel="Time (second)", ylabel="concentration (ppb)", legend=:outerrig ## Parameters The parameters in the model are: ```julia @example 1 -parameters(sys) # [cloudFrac, qrain, ρ_air, Δz] +parameters(sys) #[ρ_air, qrain, Δz, cloudFrac] ``` where ```cloudFrac``` is fraction of grid cell covered by clouds, ```qrain``` is rain mixing ratio, ```ρ_air``` is air density [kg/m3], and ```Δz``` is fall distance [m]. Let's run some simulation with different value for parameter ```cloudFrac```. ```@example 1 -@unpack O3 = sys +@unpack O3,cloudFrac,qrain = sys -p1 = [0.3,0.5,1.204,200] -p2 = [0.6,0.5,1.204,200] +p1 = [cloudFrac=>0.3] +p2 = [cloudFrac=>0.6] sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0) sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0) @@ -68,12 +68,12 @@ From the plot we could see that with larger cloud fraction, the wet deposition r Let's run some simulation with different value for parameter ```qrain``` ```@example 1 -p3 = [0.5,0.3,1.204,200] -p4 = [0.5,0.6,1.204,200] +p3 = [qrain=>0.3] +p4 = [qrain=>0.6] sol3 = solve(ODEProblem(sys, u0, tspan, p3),AutoTsit5(Rosenbrock23()), saveat=10.0) sol4 = solve(ODEProblem(sys, u0, tspan, p4),AutoTsit5(Rosenbrock23()), saveat=10.0) -plot([sol3[O3],sol4[O3]], label = ["cloudFrac=0.3" "cloudFrac=0.6"], title = "Change of O3 concentration due to wet deposition", xlabel="Time (second)", ylabel="concentration (ppb)") +plot([sol3[O3],sol4[O3]], label = ["qrain=0.3" "qrain=0.6"], title = "Change of O3 concentration due to wet deposition", xlabel="Time (second)", ylabel="concentration (ppb)") ``` The graph indicates that an increase in the rain mixing ratio leads to a corresponding rise in the rate of wet deposition. From d058416012925fe6954452837db5a0a461121a71 Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Tue, 1 Oct 2024 13:34:49 -0500 Subject: [PATCH 09/15] Change deposition species --- docs/src/dry_deposition.md | 10 +++++----- docs/src/emep.md | 3 +-- ext/GasChemExt.jl | 5 ----- src/dry_deposition.jl | 14 ++++++-------- src/wet_deposition.jl | 20 ++++++-------------- 5 files changed, 18 insertions(+), 34 deletions(-) diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md index 8b449319..38959c98 100644 --- a/docs/src/dry_deposition.md +++ b/docs/src/dry_deposition.md @@ -19,9 +19,9 @@ model = DrydepositionG() Before running any simulations with the model we need to convert it into a system of differential equations. ```julia sys = structural_simplify(model) -tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,5,2.34,0.15,10] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O, HNO₃ -sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters +tspan = (0.0, 3600*24) +prob = ODEProblem(sys, [], tspan, []) # default initial concentration of SO₂, O₃, NO₂, H₂O₂, HNO₃, CH₂O +sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` ```@setup 1 @@ -36,8 +36,8 @@ model = DrydepositionG() sys = structural_simplify(model) tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,5,2.34,0.15,10] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O, HNO₃ -sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters +prob = ODEProblem(sys, [], tspan, []) # default initial concentration of SO₂, O₃, NO₂, H₂O₂, HNO₃, CH₂O +sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` which we can plot as diff --git a/docs/src/emep.md b/docs/src/emep.md index 865688e4..d1711e90 100644 --- a/docs/src/emep.md +++ b/docs/src/emep.md @@ -35,8 +35,7 @@ model = Wetdeposition() sys = structural_simplify(model) tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,1400,275,50,0.15,2.34,10,0.15] # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP, H₂O₂, HNO₃, CH₂O -prob = ODEProblem(sys, u0, tspan, []) +prob = ODEProblem(sys, [], tspan, []) # default initial concentration of SO₂, O₃, NO₂, H₂O₂, HNO₃, CH₂O sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` diff --git a/ext/GasChemExt.jl b/ext/GasChemExt.jl index 6add2acf..0e9f96a3 100644 --- a/ext/GasChemExt.jl +++ b/ext/GasChemExt.jl @@ -9,7 +9,6 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos c.SO2 => d.SO2, c.NO2 => d.NO2, c.O3 => d.O3, - c.NO => d.NO, c.H2O2 => d.H2O2, c.CH2O => d.CH2O, )) @@ -22,10 +21,6 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos c.SO2 => d.SO2, c.NO2 => d.NO2, c.O3 => d.O3, - c.CH4 => d.CH4, - c.CO => d.CO, - c.DMS => d.DMS, - c.ISOP => d.ISOP, c.H2O2 => d.H2O2, c.CH2O => d.CH2O, )) diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 7adff712..cc4bdda2 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -254,20 +254,18 @@ function DrydepositionG(; name=:DrydepositionG) D = Differential(t) vars = @variables( - SO2(t), [unit = u"ppb"], - O3(t),[unit = u"ppb"], - NO2(t), [unit = u"ppb"], - NO(t), [unit = u"ppb"], - H2O2(t), [unit = u"ppb"], - CH2O(t), [unit = u"ppb"], - HNO3(t), [unit = u"ppb"], + SO2(t) = 2, [unit = u"ppb"], + O3(t) = 10,[unit = u"ppb"], + NO2(t) = 10, [unit = u"ppb"], + H2O2(t) = 2.34, [unit = u"ppb"], + CH2O(t) = 0.15, [unit = u"ppb"], + HNO3(t) = 10, [unit = u"ppb"], ) eqs = [ D(SO2) ~ -DryDepGas(z, z₀, u_star, L, ρA, So2Data, G, T, θ, iSeason, iLandUse, rain, dew, true, false) / z * SO2 D(O3) ~ -DryDepGas(z, z₀, u_star, L, ρA, O3Data, G, T, θ, iSeason, iLandUse, rain, dew, false, true) / z * O3 D(NO2) ~ -DryDepGas(z, z₀, u_star, L, ρA, No2Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * NO2 - D(NO) ~ -DryDepGas(z, z₀, u_star, L, ρA, NoData, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * NO D(H2O2) ~ -DryDepGas(z, z₀, u_star, L, ρA, H2o2Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * H2O2 D(CH2O) ~ -DryDepGas(z, z₀, u_star, L, ρA, HchoData, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * CH2O D(HNO3) ~ -DryDepGas(z, z₀, u_star, L, ρA, Hno3Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * HNO3 diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index 83af024e..6f69695a 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -65,26 +65,18 @@ function Wetdeposition(; name=:Wetdeposition) D = Differential(t) vars = @variables( - SO2(t), [unit = u"ppb"], - O3(t), [unit = u"ppb"], - NO2(t), [unit = u"ppb"], - CH4(t), [unit = u"ppb"], - CO(t), [unit = u"ppb"], - DMS(t), [unit = u"ppb"], - ISOP(t), [unit = u"ppb"], - H2O2(t), [unit = u"ppb"], - HNO3(t), [unit = u"ppb"], - CH2O(t), [unit = u"ppb"], + SO2(t) = 2, [unit = u"ppb"], + O3(t) = 10, [unit = u"ppb"], + NO2(t) = 10, [unit = u"ppb"], + H2O2(t) = 2.34, [unit = u"ppb"], + HNO3(t) = 10, [unit = u"ppb"], + CH2O(t) = 0.15, [unit = u"ppb"], ) eqs = [ D(SO2) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[2] * SO2 D(O3) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * O3 D(NO2) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * NO2 - D(CH4) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * CH4 - D(CO) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * CO - D(DMS) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * DMS - D(ISOP) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * ISOP D(H2O2) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * H2O2 D(HNO3) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * HNO3 D(CH2O) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * CH2O From d45d5b7e94ed6db6e6bd32b27541fa8157d17dd2 Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Tue, 1 Oct 2024 13:52:56 -0500 Subject: [PATCH 10/15] Fix typo --- docs/src/dry_deposition.md | 4 ++-- docs/src/emep.md | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md index 38959c98..31e58ad6 100644 --- a/docs/src/dry_deposition.md +++ b/docs/src/dry_deposition.md @@ -59,8 +59,8 @@ Let's run some simulation with different value for parameter ```z```. p1 = [z=>50] p2 = [z=>10] -sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0) -sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0) +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 = ["z=50m" "z=10m"], title = "Change of O3 concentration due to dry deposition", xlabel="Time (second)", ylabel="concentration (ppb)") ``` diff --git a/docs/src/emep.md b/docs/src/emep.md index d1711e90..81f7af50 100644 --- a/docs/src/emep.md +++ b/docs/src/emep.md @@ -58,8 +58,8 @@ Let's run some simulation with different value for parameter ```cloudFrac```. p1 = [cloudFrac=>0.3] p2 = [cloudFrac=>0.6] -sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0) -sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0) +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 = ["cloudFrac=0.3" "cloudFrac=0.6"], title = "Change of O3 concentration due to wet deposition", xlabel="Time (second)", ylabel="concentration (ppb)") ``` @@ -69,8 +69,8 @@ Let's run some simulation with different value for parameter ```qrain``` ```@example 1 p3 = [qrain=>0.3] p4 = [qrain=>0.6] -sol3 = solve(ODEProblem(sys, u0, tspan, p3),AutoTsit5(Rosenbrock23()), saveat=10.0) -sol4 = solve(ODEProblem(sys, u0, tspan, p4),AutoTsit5(Rosenbrock23()), saveat=10.0) +sol3 = solve(ODEProblem(sys, [], tspan, p3),AutoTsit5(Rosenbrock23()), saveat=10.0) +sol4 = solve(ODEProblem(sys, [], tspan, p4),AutoTsit5(Rosenbrock23()), saveat=10.0) plot([sol3[O3],sol4[O3]], label = ["qrain=0.3" "qrain=0.6"], title = "Change of O3 concentration due to wet deposition", xlabel="Time (second)", ylabel="concentration (ppb)") ``` From ebccc8e8acb247aa7799aefcc5acdb652b87f7de Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Fri, 1 Nov 2024 16:57:54 -0500 Subject: [PATCH 11/15] Update packages version --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index bcf6e927..4970a66a 100644 --- a/Project.toml +++ b/Project.toml @@ -24,8 +24,8 @@ GasChemExt = "GasChem" [compat] DifferentialEquations = "7" -EarthSciData = "0.9" -EarthSciMLBase = "0.16" +EarthSciData = "0.10" +EarthSciMLBase = "0.19" GasChem = "0.7" ModelingToolkit = "9" SafeTestsets = "0.1" From a2af6b52ae22e44b62daf461da9b6da7f1b57b4d Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Fri, 1 Nov 2024 17:36:55 -0500 Subject: [PATCH 12/15] Change geos-fp --- test/connector_test.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/test/connector_test.jl b/test/connector_test.jl index 20badd7b..90752cf4 100644 --- a/test/connector_test.jl +++ b/test/connector_test.jl @@ -2,6 +2,11 @@ using AtmosphericDeposition using Test, DynamicQuantities, ModelingToolkit, GasChem, Dates, EarthSciMLBase, EarthSciData using ModelingToolkit:t +domain = DomainInfo(DateTime(2022, 1, 1), DateTime(2022, 1, 3); + latrange=deg2rad(-85.0f0):deg2rad(2):deg2rad(85.0f0), + lonrange=deg2rad(-180.0f0):deg2rad(2.5):deg2rad(175.0f0), + levrange=1:10, dtype=Float64) + @testset "GasChemExt" begin start = Dates.datetime2unix(Dates.DateTime(2016, 5, 1)) composed_ode = couple(SuperFast(), FastJX(), DrydepositionG(), Wetdeposition()) @@ -20,7 +25,8 @@ end @parameters lat = deg2rad(40.0f0) [unit=u"rad"] @parameters lon = deg2rad(-97.0f0) [unit=u"rad"] @parameters lev = 1 - geosfp = GEOSFP("4x5") + + geosfp = GEOSFP("4x5", domain) model = couple(SuperFast(), FastJX(), geosfp, Wetdeposition(), DrydepositionG()) From 1aa17d883921ab54868aa206b13a8176f0b08cca Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Wed, 20 Nov 2024 01:29:05 +0000 Subject: [PATCH 13/15] CompatHelper: add new compat entry for DataInterpolations at version 6, (keep existing compat) --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 4970a66a..639c3e1f 100644 --- a/Project.toml +++ b/Project.toml @@ -23,6 +23,7 @@ EarthSciDataExt = "EarthSciData" GasChemExt = "GasChem" [compat] +DataInterpolations = "6" DifferentialEquations = "7" EarthSciData = "0.10" EarthSciMLBase = "0.19" From 727865437732f42c8bcf125ef990e2229bbfcb69 Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Wed, 20 Nov 2024 01:29:10 +0000 Subject: [PATCH 14/15] CompatHelper: add new compat entry for DynamicQuantities at version 1, (keep existing compat) --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 4970a66a..5a32e904 100644 --- a/Project.toml +++ b/Project.toml @@ -24,6 +24,7 @@ GasChemExt = "GasChem" [compat] DifferentialEquations = "7" +DynamicQuantities = "1" EarthSciData = "0.10" EarthSciMLBase = "0.19" GasChem = "0.7" From a829c703633301f59e23bfc9291ed5612a17c1d1 Mon Sep 17 00:00:00 2001 From: Christopher Tessum Date: Fri, 6 Dec 2024 17:59:00 -0600 Subject: [PATCH 15/15] Update dependencies and connections --- .github/workflows/CI.yml | 2 +- Project.toml | 8 ++++---- ext/GasChemExt.jl | 6 +++--- test/connector_test.jl | 17 +++++++---------- 4 files changed, 15 insertions(+), 18 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 65810317..da0cff09 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -10,8 +10,8 @@ jobs: fail-fast: false matrix: version: - - '1.9' - '1.10' + - '1.11' os: - ubuntu-latest arch: diff --git a/Project.toml b/Project.toml index 4970a66a..bb00117c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "AtmosphericDeposition" uuid = "1a52f20c-0d16-41d8-a00a-b2996d86a462" authors = ["EarthSciML authors and contributors"] -version = "0.2.1" +version = "0.2.2" [deps] DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" @@ -24,9 +24,9 @@ GasChemExt = "GasChem" [compat] DifferentialEquations = "7" -EarthSciData = "0.10" -EarthSciMLBase = "0.19" -GasChem = "0.7" +EarthSciData = "0.12" +EarthSciMLBase = "0.20" +GasChem = "0.9" ModelingToolkit = "9" SafeTestsets = "0.1" StaticArrays = "1" diff --git a/ext/GasChemExt.jl b/ext/GasChemExt.jl index 0e9f96a3..e81e6dc0 100644 --- a/ext/GasChemExt.jl +++ b/ext/GasChemExt.jl @@ -6,7 +6,7 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos c, d = c.sys, d.sys operator_compose(convert(ODESystem, c), d, Dict( - c.SO2 => d.SO2, + #c.SO2 => d.SO2, # SuperFast does not currently have SO2 c.NO2 => d.NO2, c.O3 => d.O3, c.H2O2 => d.H2O2, @@ -18,7 +18,7 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos c, d = c.sys, d.sys operator_compose(convert(ODESystem, c), d, Dict( - c.SO2 => d.SO2, + #c.SO2 => d.SO2, # SuperFast does not currently have SO2 c.NO2 => d.NO2, c.O3 => d.O3, c.H2O2 => d.H2O2, @@ -26,4 +26,4 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos )) end -end \ No newline at end of file +end diff --git a/test/connector_test.jl b/test/connector_test.jl index 90752cf4..d4d08e79 100644 --- a/test/connector_test.jl +++ b/test/connector_test.jl @@ -10,11 +10,9 @@ domain = DomainInfo(DateTime(2022, 1, 1), DateTime(2022, 1, 3); @testset "GasChemExt" begin start = Dates.datetime2unix(Dates.DateTime(2016, 5, 1)) composed_ode = couple(SuperFast(), FastJX(), DrydepositionG(), Wetdeposition()) - combined_mtk = convert(ODESystem, composed_ode) - sys = structural_simplify(combined_mtk) + sys = convert(ODESystem, composed_ode) print(unknowns(sys)) - @test length(unknowns(sys)) ≈ 20 - #TODO Change 20 to 18 after the latest GasChem package include species HNO3 + @test length(unknowns(sys)) ≈ 15 eqs = string(equations(sys)) wanteqs = ["Differential(t)(SuperFast₊O3(t)) ~ SuperFast₊DrydepositionG_ddt_O3ˍt(t) + SuperFast₊Wetdeposition_ddt_O3ˍt(t)"] @@ -27,18 +25,17 @@ end @parameters lev = 1 geosfp = GEOSFP("4x5", domain) - + model = couple(SuperFast(), FastJX(), geosfp, Wetdeposition(), DrydepositionG()) - sys = structural_simplify(convert(ODESystem, model)) - @test length(unknowns(sys)) ≈ 20 - #TODO Change 20 to 18 after the latest GasChem package include species HNO3 + sys = convert(ODESystem, model) + @test length(unknowns(sys)) ≈ 15 - eqs = string(equations(convert(ODESystem, model))) + eqs = string(observed(sys)) wanteq = "DrydepositionG₊G(t) ~ GEOSFP₊A1₊SWGDN(t)" @test contains(eqs, wanteq) wanteq = "Wetdeposition₊cloudFrac(t) ~ GEOSFP₊A3cld₊CLOUD(t)" @test contains(eqs, wanteq) wanted = "Wetdeposition₊ρA(t) ~ GEOSFP₊P/(GEOSFP₊I3₊T*R)*kgperg*MW_air" @test contains(eqs, wanteq) -end \ No newline at end of file +end