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