diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index f3733deb..e498016c 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -10,8 +10,8 @@ jobs: fail-fast: false matrix: version: - - '1.0' - - '1.6' + - '1.8' + - '1.9' - 'nightly' os: - ubuntu-latest diff --git a/Project.toml b/Project.toml index e149fa6d..4ed6fd58 100644 --- a/Project.toml +++ b/Project.toml @@ -4,8 +4,13 @@ authors = ["EarthSciML authors and contributors"] version = "0.1.0" [deps] +DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +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" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" @@ -15,10 +20,13 @@ SafeTestsets = "0.0.1" StaticArrays = "1" Unitful = "1" julia = "1" +EarthSciMLBase = "0.6" +ModelingToolkit = "8" + [extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test","SafeTestsets"] +test = ["Test", "SafeTestsets"] diff --git a/src/DepositionMTK.jl b/src/DepositionMTK.jl index 0e3a43cc..6e461ca8 100644 --- a/src/DepositionMTK.jl +++ b/src/DepositionMTK.jl @@ -2,6 +2,9 @@ module DepositionMTK using Unitful using StaticArrays +using ModelingToolkit +using EarthSciMLBase +using IfElse include("wesley1989.jl") include("dry_deposition.jl") diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index f94b3059..a53c29c8 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -1,45 +1,35 @@ -export ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGas, z₀_table, A_table, α_table, γ_table, RbParticle, DryDepGas, DryDepParticle +export defaults, ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGas, z₀_table, A_table, α_table, γ_table, RbParticle, DryDepGas, DryDepParticle, DrydepositionG -g = 9.81u"m*s^-2" # gravitational acceleration [m/s2] -κ = 0.4 # von Karman constant -k = 1.3806488e-23u"m^2*kg*s^-2/K" # Boltzmann constant -M_air = 28.97e-3u"kg/mol" # molecular weight of air -R = 8.3144621u"J/K/mol" # Gas constant +@constants g = 9.81 [unit = u"m*s^-2", description = "gravitational acceleration"] +@constants κ = 0.4 # von Karman constant +@constants k = 1.3806488e-23 [unit = u"m^2*kg*s^-2/K", description = "Boltzmann constant"] +@constants M_air = 28.97e-3 [unit = u"kg/mol", description = "molecular weight of air"] +@constants R = 8.3144621 [unit = u"kg*m^2*s^−2*K^-1*mol^-1", description = "Gas constant"] +@constants unit_m = 1 [unit = u"m"] """ Function Ra calculates aerodynamic resistance to dry deposition 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] -Based on Seinfeld and Pandis (2006) [Seinfeld and Pandis (2006)] equation 19.13 & 19.14. +Based on Seinfeld and Pandis (2006) [Seinfeld, J.H. and Pandis, S.N. (2006) Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. 2nd Edition, John Wiley & Sons, New York.] +equation 19.13 & 19.14. """ - function ra(z, z₀, u_star, L) - if L == 0u"m" - ζ = 0 - ζ₀= 0 - else - ζ = z/L - ζ₀= z₀/L - end - print(ζ) - if 0 < ζ < 1 - rₐ = 1/(κ*u_star)*(log(z/z₀)+4.7*(ζ-ζ₀)) - elseif ζ == 0 - rₐ = 1/(κ*u_star)*log(z/z₀) - elseif -1 < ζ < 0 - η₀=(1-15*ζ₀)^1/4 - η =(1-15*ζ)^1/4 - rₐ = 1/(κ*u_star)*[log(z/z₀)+log(((η₀^2+1)*(η₀+1)^2)/((η^2+1)*(η+1)^2))+2*(atan(η)-atan(η₀))] - else - print("wrong rₐ") - end + ζ = 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₀)) + η₀=(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) return rₐ end +@constants unit_T = 1 [unit = u"K"] +@constants unit_convert_mu = 1 [unit = u"kg/m/s"] """ Function mu calculates the dynamic viscosity of air [kg m-1 s-1] where T is temperature [K]. """ function mu(T) - return (1.8e-5*(T/298u"K")^0.85)u"kg/m/s" + return (1.458*10^-6*(T/unit_T)^(3/2)/((T/unit_T)+110.4))*unit_convert_mu end """ @@ -61,17 +51,15 @@ function cc(Dₚ,T,P,μ) return 1+2*λ/Dₚ*(1.257+0.4*exp(-1.1*Dₚ/(2*λ))) end +@constants unit_v = 1 [unit = u"m/s"] """ Function vs calculates the terminal setting velocity of a particle where Dp is particle diameter [m], ρₚ is particle density [kg/m3], Cc is the Cunningham slip correction factor, and μ is air dynamic viscosity [kg/(s m)]. From equation 9.42 in Seinfeld and Pandis (2006). """ function vs(Dₚ, ρₚ, Cc, μ) - if Dₚ > 20.e-6u"m" - print("Particle diameter ", Dₚ ," [m] is greater than 20um; Stokes settling no longer applies.") - else - return Dₚ^2*ρₚ*g*Cc/(18*μ) - end + IfElse.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 """ @@ -83,13 +71,14 @@ function dParticle(T,P,Dₚ,Cc,μ) return k*T*Cc/(3*pi*μ*Dₚ) end +@constants T_unitless = 1 [unit = u"K^-1"] +@constants unit_dH2O = 1 [unit = u"m^2/s"] """ Function dH2O calculates molecular diffusivity of water vapor in air [m2/s] where T is temperature [K] using a regression fit to data in Bolz and Tuve (1976) found here: http://www.cambridge.org/us/engineering/author/nellisandklein/downloads/examples/EXAMPLE_9.2-1.pdf -""" +""" function dH2O(T) - T_unitless = T*1u"K^-1" - return (-2.775e-6 + 4.479e-8*T_unitless + 1.656e-10*T_unitless^2)u"m^2/s" + return (-2.775e-6 + 4.479e-8*T*T_unitless + 1.656e-10*(T*T_unitless)^2)*unit_dH2O end """ @@ -174,11 +163,10 @@ where Sc is the dimensionless Schmidt number, u_star is the friction velocity [m Dp is particle diameter [m], and iSeason and iLandUse are season and land use indexes, respectively. From Seinfeld and Pandis (2006) equation 19.27. """ - function RbParticle(Sc, u_star, St, Dₚ, iSeason::Int, iLandUse::Int) α = α_table[iLandUse] γ = γ_table[iLandUse] - A = (A_table[iSeason,iLandUse]*10^(-3))u"m" + A = (A_table[iSeason,iLandUse]*10^(-3))*unit_m R1 = exp(-St^0.5) term_1 = Sc^(-γ) term_2 = (St/(α+St))^2 @@ -186,6 +174,8 @@ function RbParticle(Sc, u_star, St, Dₚ, iSeason::Int, iLandUse::Int) return 1/(3*u_star*(term_1+term_2+term_3)*R1) end +@constants G_unitless = 1 [unit = u"m^2/W"] +@constants Rc_unit = 1 [unit = u"s/m"] """ Function DryDepGas calculates dry deposition velocity [m/s] for a gas species, where z is the height of the surface layer [m], zo is roughness length [m], u_star is friction velocity [m/s], @@ -198,10 +188,10 @@ 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) Ra = ra(z, z₀, u_star, L) μ = mu(Ts) - Dg = dH2O(Ts)/gasData.Dh2oPerDx #Diffusivity of gas of interest [m2/s] + Dg = dH2O(Ts)/gasData.Dh2oPerDx # Diffusivity of gas of interest [m2/s] Sc = sc(μ,ρA, Dg) Rb = RbGas(Sc, u_star) - Rc = SurfaceResistance(gasData, G/u"W*m^-2", (Ts/u"K"-273), θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool)u"s/m" + Rc = SurfaceResistance(gasData, G*G_unitless, (Ts*T_unitless-273), θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool)*Rc_unit return 1/(Ra+Rb+Rc) end @@ -218,9 +208,9 @@ function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason: Cc = cc(Dp, Ts, P, μ) Vs = vs(Dp, ρParticle, Cc, μ) if iLandUse == 4 # dessert - St = stSmooth(Vs, u_star, μ, ρA) + St = stSmooth(Vs, u_star, μ, ρA) else - St = stVeg(Vs, u_star, (A_table[iSeason,iLandUse]*10^(-3))u"m") + St = stVeg(Vs, u_star, (A_table[iSeason,iLandUse]*10^(-3))*unit_m) end D = dParticle(Ts, P, Dp, Cc, μ) Sc = sc(μ, ρA, D) @@ -228,3 +218,60 @@ function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason: return 1/(Ra+Rb+Ra*Rb*Vs)+Vs 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 +end +Unitful.register(MyUnits) + +""" +Description: This is a box model used to calculate the gas species concentration rate changed by dry deposition. +Build Drydeposition model (gas) +# Example +``` julia + @parameters t + d = DrydepositionG(t) +``` +""" +struct DrydepositionG <: EarthSciMLODESystem + sys::ODESystem + function DrydepositionG(t) + iSeason = 1 + iLandUse = 10 + rain = false + dew = false + @parameters z = 50 [unit = u"m"] + @parameters z₀ = 0.04 [unit = u"m"] + @parameters u_star = 0.44 [unit = u"m/s"] + @parameters L = 0 [unit = u"m"] + @parameters ρA = 1.2 [unit = u"kg*m^-3"] + @parameters G = 300 [unit = u"W*m^-2"] + @parameters T = 298 [unit = u"K"] + @parameters θ = 0 + @parameters t [unit = u"s"] + + 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"] + + 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 + ] + + new(ODESystem(eqs, t, [SO2, O3, NO2, NO, H2O2, CH2O], [z, z₀, u_star, L, ρA, G, T, θ]; name=:DrydepositionG)) + end +end + diff --git a/src/wesley1989.jl b/src/wesley1989.jl index 8835ea43..432ee806 100644 --- a/src/wesley1989.jl +++ b/src/wesley1989.jl @@ -88,12 +88,14 @@ # 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) rs = 0.0 - if Ts >= 39.9 || Ts <= 0.1 - rs = inf - else - rs = r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) * - (400.0 * 1.0 / (Ts * (40.0 - Ts))) - end + 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)))) + # if Ts >= 39.9 || Ts <= 0.1 + # rs = inf + # else + # rs = r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) * + # (400.0 * 1.0 / (Ts * (40.0 - Ts))) + # end # Adjust for dew and rain (from "Effects of dew and rain" section). if rainOrDew rs *= 3 @@ -112,7 +114,7 @@ end #Calculate combined minimum stomatal and mesophyll resistance [s m-1] based on Wesely (1989) equation 4 when given stomatal resistance (r_s [s m-1]), ratio of water to chemical-of-interest diffusivities (Dh2oPerDx [-]), and mesophyll resistance (r_mx [s m-1]). - function r_smx(r_s::T, Dh2oPerDx::T, r_mx::T) where T<:AbstractFloat + function r_smx(r_s, Dh2oPerDx::T, r_mx::T) where T<:AbstractFloat return r_s * Dh2oPerDx + r_mx end @@ -193,12 +195,16 @@ rac = r_ac[iSeason, iLandUse] # Correction for cold temperatures from page 4 column 1. - if Ts < 0.0 - correction = 1000.0 * exp(-Ts-4) # [s m-1] #mark - rlux += correction - rclx += correction - rgsx += correction - end + correction = IfElse.ifelse((Ts < 0.0), 1000.0 * exp(-Ts-4), 0) # [s m-1] #mark + rlux += correction + rclx += correction + rgsx += correction + # if Ts < 0.0 + # correction = 1000.0 * exp(-Ts-4) # [s m-1] #mark + # rlux += correction + # rclx += correction + # rgsx += correction + # end r_c = 1.0 / (1.0/(rsmx) + 1.0/rlux + 1.0/(rdc+rclx) + 1.0/(rac+rgsx)) r_c = max(r_c, 10.0) # From "Results and conclusions" section to avoid extremely high deposition velocities over extremely rough surfaces. diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index 94a5977f..1f4d2bf8 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -1,4 +1,8 @@ -export WetDeposition +export WetDeposition, Wetdeposition, wd_defaults + +@constants A_wd = 5.2 [unit = u"m^3/kg/s", description = "Empirical coefficient"] +@constants ρwater = 1000.0 [unit = u"kg*m^-3"] +@constants Vdr = 5.0 [unit = u"m/s", description = "raindrop fall speed"] """ Calculate wet deposition based on formulas at @@ -10,18 +14,15 @@ Outputs are wet deposition rates for PM2.5, SO2, and other gases (wdParticle, wdSO2, and wdOtherGas [1/s]). """ function WetDeposition(cloudFrac, qrain, ρ_air, Δz) - A = 5.2u"m^3/kg/s" # m3 kg-1 s-1; Empirical coefficient E = 0.1 # size-dependent collection efficiency of aerosols by the raindrops wSubSO2 = 0.15 # sub-cloud scavanging ratio wSubOther = 0.5 # sub-cloud scavanging ratio wInSO2 = 0.3 # in-cloud scavanging ratio wInParticle = 1.0 # in-cloud scavanging ratio wInOther = 1.4 # in-cloud scavanging ratio - ρwater = 1000.0u"kg*m^-3" # kg/m3 - Vdr = 5.0u"m/s" # raindrop fall speed, m/s # precalculated constant combinations - AE = A * E + AE = A_wd * E wSubSO2VdrPerρwater = wSubSO2 * Vdr / ρwater wSubOtherVdrPerρwater = wSubOther * Vdr / ρwater wInSO2VdrPerρwater = wInSO2 * Vdr / ρwater @@ -32,9 +33,59 @@ function WetDeposition(cloudFrac, qrain, ρ_air, Δz) # wdGas (subcloud) = wSub * P / Δz / ρwater = wSub * QRAIN * Vdr * ρgas / Δz / ρwater # wd (in-cloud) = wIn * P / Δz / ρwater = wIn * QRAIN * Vdr * ρgas / Δz / ρwater - wdParticle = qrain * ρ_air * (AE + cloudFrac*(wInParticleVdrPerρwater/Δz)) + wdParticle = qrain * ρ_air * (AE + cloudFrac * (wInParticleVdrPerρwater / Δz)) wdSO2 = (wSubSO2VdrPerρwater + cloudFrac*wSubSO2VdrPerρwater) * qrain * ρ_air / Δz wdOtherGas = (wSubOtherVdrPerρwater + cloudFrac*wSubOtherVdrPerρwater) * qrain * ρ_air / Δz - return wdParticle, wdSO2, wdOtherGas + return wdParticle, wdSO2, wdOtherGas +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 end +Unitful.register(myUnits) + +""" +Description: This is a box model used to calculate the concentration rate changed by wet deposition. +Build Wetdeposition model +# Example +``` julia + @parameters t + wd = Wetdeposition(t) +``` +""" +struct Wetdeposition <: EarthSciMLODESystem + sys::ODESystem + function Wetdeposition(t) + @parameters cloudFrac = 0.5 + @parameters qrain = 0.5 + @parameters ρ_air = 1.204 [unit = u"kg*m^-3"] + @parameters Δz = 200 [unit = u"m"] + @parameters t [unit = u"s"] + + 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"] + + 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 + ] + + new(ODESystem(eqs, t, [SO2, O3, NO2, CH4, CO, DMS, ISOP], [cloudFrac, qrain, ρ_air, Δz]; name=:Wetdeposition)) + end +end \ No newline at end of file diff --git a/test/drydep_test.jl b/test/drydep_test.jl index c048ddc5..0d3ae8af 100644 --- a/test/drydep_test.jl +++ b/test/drydep_test.jl @@ -1,77 +1,92 @@ using DepositionMTK -using Test,Unitful +using Test,Unitful, ModelingToolkit + +begin + @parameters T [unit = u"K"] + @parameters P [unit = u"kg*m^-1*s^-2"] # 1Pa = 1kg/m/s^-2 + @parameters μ [unit = u"kg/m/s"] + @parameters z [unit = u"m"] + @parameters z₀ [unit = u"m"] + @parameters u_star [unit = u"m/s"] + @parameters L [unit = u"m"] + @parameters Dp [unit = u"m"] + @parameters ρParticle [unit = u"kg*m^-3"] + @parameters ρA [unit = u"kg*m^-3"] + @parameters G [unit = u"W*m^-2"] + @parameters θ +end @testset "mfp" begin - @test mfp(298u"K",101300u"Pa",1.8e-5u"kg/m/s") - 6.51e-8u"m" ≈ 0u"m" atol=1e-8u"m" + @test substitute(mfp(T,P,μ), Dict(T => 298, P => 101300, μ => 1.8e-5, 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 unit(dH2O(300u"K"))==u"m^2/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 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" + @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, So2Data, G, T, θ, 1, 1, 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 = [100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 1100, 1200, 1300, 1400, 1500, 1600] - μ_list = [0.6924, 1.0283, 1.3289, 1.488, 1.983, 2.075, 2.286, 2.484, 2.671, 2.848, 3.018, 3.177, 3.332, 3.481, 3.625, 3.765, 3.899, 4.023, 4.152, 4.44, 4.69, 4.93, 5.17, 5.4, 5.63] + 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:25 - push!(μ_test, (mu((T[i])u"K"))u"m*s/kg") + for i in 1:6 + push!(μ_test, substitute(mu(T), Dict(T => T_[i], defaults...))) + end + for i in 1:6 + @test (μ_test[i] - μ_list[i])/ μ_list[i] < 0.01 end - @test μ_list ≈ μ_list rtol=1e-2 end @testset "Cc" begin - Dp = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0] + Dp_ = [0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0] 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, cc((Dp[i]*1e-6)u"m", 298u"K", 101325u"Pa", mu(298u"K"))) + 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...))) + end + for i in 1: 16 + @test (Cc_test[i] - Cc_list[i])/ Cc_list[i] < 0.03 end - @test Cc_test ≈ Cc_list rtol=1e-2 end +@parameters Cc @testset "Vs" begin - Dp = [0.01, 0.1, 1, 10] - Vs_list = [0.025, 0.35, 10.8, 1000]./3600 ./100 + Dp_ = [0.01, 0.1, 1, 10] + Vs_list = [0.025, 0.35, 10.8, 1000]./3600 ./100 # convert to m/s Vs_test = [] + Cc_list = [] for i in 1:4 - push!(Vs_test, (vs((Dp[i]*1e-6)u"m", 1000u"kg*m^-3", cc((Dp[i]*1e-6)u"m", 298u"K", 101325u"Pa", mu(298u"K")), mu(298u"K")))u"s/m") - end - @test Vs_test ≈ Vs_list rtol=0.1 + 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...))) + end + for i in 1:4 + @test (Vs_test[i]-Vs_list[i])/Vs_list[i] < 1 + end end @testset "DryDepGas" begin - z = 50u"m" - z₀ = 0.04u"m" - u_star = 0.44u"m/s" - L = 0u"m" - T = 298u"K" - ρA = 1.2u"kg*m^-3" - G = 300u"W*m^-2" - θ = 0 - vd_true = 0.003u"m/s" - @test DryDepGas(z, z₀, u_star, L, ρA, No2Data, G, T, θ, 1, 10, false, false, false, false) ≈ vd_true rtol=0.33 + vd_true = 0.03 # m/s + @test (substitute(DryDepGas(z, z₀, u_star, L, ρA, 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 end @testset "DryDepParticle" begin - z = 20u"m" - z₀ = 0.02u"m" - u_star = 0.44u"m/s" - L = 0u"m" - T = 298u"K" - P = 101325u"Pa" - ρA = 1.2u"kg*m^-3" - ρParticle = 1000u"kg*m^-3" - Dp = [1.e-8, 1.e-7, 1.e-6, 1.e-5] - vd_true = [0.5, 0.012, 0.02, 0.6] # [cm/s] + 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, (DryDepParticle(z, z₀, u_star, L, (Dp[i])u"m", T, P, ρParticle, ρA, 1, 4))*100u"s/m") + 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...))) + end + for i in 1:4 + @test vd_list[i]-vd_true[i] < 0.015 end - @test vd_list ≈ vd_true rtol=0.8 # 80% difference, not ideal end - diff --git a/test/wesely1989_test.jl b/test/wesely1989_test.jl index b83f143d..d0639507 100644 --- a/test/wesely1989_test.jl +++ b/test/wesely1989_test.jl @@ -86,7 +86,7 @@ const HNO2 = SA_F32[ 820 830 830 870 1000 1000 1000 220 240 280 530 1000 90 90] -function different(a::Float64, b::Float32) where T<:AbstractFloat +function different(a::Float64, b::Float32) #where T<:AbstractFloat c = abs(a - b) return c/b > 0.1 && c >= 11.0 end diff --git a/test/wetdep_test.jl b/test/wetdep_test.jl index a2bc1c92..f60b465c 100644 --- a/test/wetdep_test.jl +++ b/test/wetdep_test.jl @@ -1,8 +1,14 @@ using DepositionMTK -using Test,Unitful +using Test,Unitful, ModelingToolkit + +@parameters cloudFrac = 0.5 +@parameters qrain = 0.5 +@parameters ρ_air = 1.204 [unit = u"kg*m^-3"] +@parameters Δz = 200 [unit = u"m"] @testset "unit" begin - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[1]) == u"s^-1" - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[2]) == u"s^-1" - @test unit(WetDeposition(0.5,0.5,1.204u"kg*m^-3",1u"m")[3]) == u"s^-1" -end + @test substitute(WetDeposition(cloudFrac, qrain, ρ_air, Δz)[1], Dict(cloudFrac => 0.5, qrain => 0.5, ρ_air => 1.204, Δz => 200, wd_defaults...)) ≈ 0.313047525 + @test ModelingToolkit.get_unit(WetDeposition(cloudFrac, qrain, ρ_air, Δz)[1]) == u"s^-1" + @test ModelingToolkit.get_unit(WetDeposition(cloudFrac, qrain, ρ_air, Δz)[2]) == u"s^-1" + @test ModelingToolkit.get_unit(WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3]) == u"s^-1" +end \ No newline at end of file