Skip to content

Commit

Permalink
Update Deposition to be compatible with MTKv9
Browse files Browse the repository at this point in the history
  • Loading branch information
@jialinl6
jialinl6 committed Sep 26, 2024
1 parent 69aa7a9 commit e8d87bf
Show file tree
Hide file tree
Showing 10 changed files with 68 additions and 76 deletions.
13 changes: 5 additions & 8 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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]
Expand Down
20 changes: 8 additions & 12 deletions ext/EarthSciDataExt.jl
View file
Original file line number Diff line number Diff line change
@@ -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"],
Expand All @@ -27,19 +25,17 @@ 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

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"],
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"],
)
Expand All @@ -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

Expand Down
2 changes: 1 addition & 1 deletion ext/GasChemExt.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down
6 changes: 4 additions & 2 deletions src/AtmosphericDeposition.jl
View file
Original file line number Diff line number Diff line change
@@ -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")
Expand Down
26 changes: 14 additions & 12 deletions src/dry_deposition.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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(
Expand All @@ -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 = [
Expand Down
5 changes: 3 additions & 2 deletions src/wesley1989.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
16 changes: 8 additions & 8 deletions src/wet_deposition.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"],
Expand All @@ -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 = [
Expand Down
32 changes: 15 additions & 17 deletions test/connector_test.jl
View file
Original file line number Diff line number Diff line change
@@ -1,37 +1,35 @@
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)"]
@test contains(string(eqs), wanteqs[1])
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
22 changes: 9 additions & 13 deletions test/drydep_test.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
using AtmosphericDeposition
using Test, Unitful, ModelingToolkit
using Test, DynamicQuantities, ModelingToolkit

begin
@parameters T [unit = u"K"]
Expand All @@ -17,27 +17,23 @@ 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
T_ = [275, 300, 325, 350, 375, 400]
μ_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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -73,15 +69,15 @@ 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
Dp_ = [1.e-8, 1.e-7, 1.e-6, 1.e-5]
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
Expand Down
2 changes: 1 addition & 1 deletion test/wetdep_test.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
using AtmosphericDeposition
using Test, Unitful, ModelingToolkit
using Test, DynamicQuantities, ModelingToolkit

@parameters cloudFrac = 0.5
@parameters qrain = 0.5
Expand Down

0 comments on commit e8d87bf

Please sign in to comment.