Merge pull request #27 from jialinl6/Fix23

Fix23

Project.toml

@@ -35,6 +35,8 @@ julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 
 [targets]
-test = ["Test"]
+test = ["Test","NonlinearSolve"]
+

docs/make.jl

@@ -16,6 +16,7 @@ makedocs(;
     repo = Remotes.GitHub("EarthSciML", "Aerosol.jl"),
     pages=[
         "Home" => "index.md",
+        "VBS" => "VBS.md",
         "Thermodynamics" => [
             "Isorropia" => [
                 "Overview" => "isorropia/overview.md",

docs/src/VBS.md

@@ -0,0 +1,108 @@
+# VBS: Volatility Basis Set (VBS) model that deals with the partitioning of Secondary Organic Aerosol (SOA)
+
+## Model
+
+This is an implementation of the model to simulate the behavior of semivolatile organic compounds in the atmosphere, as described in Donahue et al. (2006):
+
+> [Donahue, N. M., Robinson, A. L., Stanier, C. O., & Pandis, S. N. (2006). Coupled partitioning, dilution, and chemical aging of semivolatile organics. *Environmental science & technology*, 40(8), 2635-2643.](https://pubs.acs.org/doi/10.1021/es052297c)
+
+We can create an instance of the model in the following manner. ```Ci``` is the input of the function, which is the organic compound concentration in all phases (air and aerosol phase), values were read from the original Figure 1.(a):
+
+```@example 2
+using Aerosol
+using NonlinearSolve, ModelingToolkit
+
+Ci = [2.5, 1.8, 4.0, 4.0, 5.8, 4.8, 6.3, 8.0] 
+ns = VBS(Ci)
+nothing # hide
+```
+We could solve and visualize the results with the following codes, Let's try reproducing Figure 1.(a): 
+```@example 2
+simplens = structural_simplify(ns)
+@unpack C_OA,ξ = simplens
+guess = [29] # initial guess of the C_OA, the total organic compound in aerosol phase
+params = []
+prob = NonlinearProblem(simplens, guess, params)
+sol = solve(prob, TrustRegion())
+nothing # hide
+```
+Figure 1.(a) shows a partitioning in ideal situations that compounds are distributed without any potential effects. The green part is the concentration of organic compounds in aerosol phase and white part is the concentration in gas phase.
+```@example 2
+using Plots
+p1 = begin 
+    x = log10.([0.01,0.1,1,10,100,1000,10000,100000])
+    p = bar(x,Ci,label="total",color=:white)
+    bar!(x,sol[ξ].*Ci,label="condensed-phase",color=:green)
+    ylabel!("Organic Mass, μg/m³")
+    xlabel!("log10(C*)")
+    title!("Typical Ambient Partitioning")
+end
+```
+Let's try reproducing Figure 1.(b): Semi-volatile emissions as they might appear near the output of a primarysource, before substantial dilution into the background atmosphere (only enough dilution to cool the emissions to ambient temperatureis assumed). The high loading leads to partitioning well into the highC* end of the distribution, as shown in brown. Note the scale oftheyaxis (mg/m³).
+```@example 2
+Ci_2 = [0.4, 0.8, 1.25, 1.7, 2.1, 2.5, 3, 3.4].*1000 # The unit of Ci is μg/m³, and the values read from the original Figure 1(b) are in units of mg/m³
+ns2 = VBS(Ci_2)
+@unpack C_OA,ξ = structural_simplify(ns2)
+prob2 = NonlinearProblem(structural_simplify(ns2), [10000], [])
+sol2 = solve(prob2, TrustRegion())
+p2 = begin
+    x = log10.([0.01,0.1,1,10,100,1000,10000,100000])
+    p = bar(x,Ci_2,label="total",color=:white)
+    bar!(x,sol2[ξ].*Ci_2,label="condensed-phase",color=:red)
+    ylabel!("Organic Mass, mg/m³")
+    xlabel!("log10(C*)")
+    title!("Cooled Fresh Emissions")
+end
+```
+Let's try reproducing Figure 1.(c): The effect of dilution by pure air on the emissionsdepicted above. The dilution factor of 1000 is indicated with a horizontal black arrow. Dilution by a factor of 1000 reduces the aerosolmass by a factor of 4000 because of repartitioning into the vapor phase
+```@example 2
+Ci_3 = [0.4, 0.8, 1.25, 1.7, 2.1, 2.5, 3, 3.4]
+@unpack C_OA,ξ = structural_simplify(VBS(Ci_3))
+sol3 = solve(NonlinearProblem(structural_simplify(VBS(Ci_3)), [10], []), TrustRegion())
+p3 = begin
+    x = log10.([0.01,0.1,1,10,100,1000,10000,100000])
+    p = bar(x,Ci_3,label="total",color=:white)
+    bar!(x,sol3[ξ].*Ci_3,label="condensed-phase",color=:red)
+    ylabel!("Organic Mass, μg/m³")
+    xlabel!("log10(C*)")
+    title!("Diluted Emissions (dilution factor = 1000)")
+end
+```
+Let's try reproducing Figure 1.(d): The effect of dilution, as depicted in panel b above but nowinto background air represented in panel a above. The partitioning of the background organic material and the fresh emissions are keptseparate,  in  green  and  brown,  only  for  illustrative  purposes.  The  vapor  portions  of  the  background  and  the  fresh  emissions  are  alsoseparated, though each is shown with a white bar. 
+```@example 2
+p4 = begin
+    x = log10.([0.4, 0.8, 1.25, 1.7, 2.1, 2.5, 3, 3.4])
+    Ci_sum = Ci .+ Ci_3
+    p = bar(x,Ci_sum,label="total",color=:white)
+    both = sol[ξ].*Ci .+ sol3[ξ].*Ci_3
+    bar!(x,both,label="diluted emission",color=:red)
+    bar!(x,sol[ξ].*Ci,label="background",color=:green)
+    ylabel!("Organic Mass, μg/m³")
+    xlabel!("log10(C*)")
+    title!("Emission & Background Mixed")
+end
+```
+## Temperature dependence
+Temperature dependece was taken into consideration. A  basis  set  of  saturationconcentrations (ranging from 1 ng to 100 mg at 300 K) changes withtemperature  according  to  the  Clausius  Clapeyron  equation. Let's try reproducing Figure 2:
+```@example 2
+Ci_t = [0.4,0.55,0.7,1.1,1.25,1.45,1.7,1.8,2] # values were read from the original Figure
+ns_t = VBS(Ci_t,[0.001,0.01,0.1,1,10,100,1000,10000,100000]) # The second input of the VBS function is different basic set at standard temperatures, ranging from 1 ng to 100 mg at 300 K
+@unpack C_OA,ξ,T = structural_simplify(ns_t)
+params1 = [T=>310] # Temperature is an adjustable parameter in the model.
+params2 = [T=>285]
+prob_310 = NonlinearProblem(structural_simplify(ns_t), [29], params1)
+prob_285 = NonlinearProblem(structural_simplify(ns_t), [29], params2)
+sol_310 = solve(prob_310, TrustRegion())
+sol_285 = solve(prob_285, TrustRegion())
+
+# Helpful visualize function
+function pl_t(sol,string::String)
+    x = log10.([0.001,0.01,0.1,1,10,100,1000,10000,100000])
+    p = bar(x,Ci_t,label="total",color=:white)
+    bar!(x,sol[ξ].*Ci_t,label="condensed-phase",color=:green)
+    ylabel!("Organic Mass, μg/m³")
+    xlabel!("log10(C*)")
+    title!(string)
+end
+plot(pl_t(sol_310,"310K"), pl_t(sol_285,"285K"), layout=(1, 2))
+```

src/Aerosol.jl

@@ -2,6 +2,7 @@ module Aerosol
 
 using Reexport
 
+include("VBS.jl")
 include("isorropia/isorropia.jl")
 @reexport using .ISORROPIA
 @reexport using .ISORROPIA.MyUnits

src/VBS.jl

@@ -0,0 +1,49 @@
+using ModelingToolkit, Unitful
+export VBS
+
+# The function to calculate Ci_star at different temperatures
+function calc_Ci_star(T, Ci_star_standard)
+    T1 = 300 # K 
+    ΔH = 100000 # J/mol
+    R = 8.314 # J/K/mol
+    return Ci_star_standard / exp(ΔH/R*(1/(T)-1/T1))
+end
+@register_symbolic calc_Ci_star(T, Ci_star_standard)
+
+# The function create a nonlinear system with input Ci (the organic compound concentration in all phases)
+function VBS(Ci)
+    @parameters T = 300 [unit = u"K", description = "temperature"]
+    @parameters T_unit = 1 [unit = u"K"]
+    @variables Ci_star[1:8] [description = "saturation concentrations"]
+
+    Ci_star_standard = [0.01,0.1,1,10,100,1000,10000,100000] # standard saturation concentrations
+    @parameters C_OA_unit = 1 [unit = u"m^3/μg", description = "make C_OA unitless"]
+    @variables C_OA [unit = u"μg/(m^3)", description = "organic compound in aerosol phase"]
+    @variables ξ[1:8] [description = "partitioning coefficient at corresponding Ci_star"]
+
+    eqs = [
+        [ξ[i] ~ 1/(1+calc_Ci_star(T/T_unit,Ci_star_standard[i])/(C_OA*C_OA_unit)) for i ∈ 1:8]
+        C_OA*C_OA_unit ~ sum(ξ.*Ci)
+    ]
+
+    NonlinearSystem(eqs, [C_OA, collect(ξ)...], [T,C_OA_unit,T_unit]; name=:VBS)
+end
+
+# The function create a nonlinear system with input Ci (the organic compound concentration in all phases) and Ci_star_standard(basic set at standard temperatures)
+function VBS(Ci,Ci_star_standard)
+    @parameters T = 300 [unit = u"K", description = "temperature"]
+    @parameters T_unit = 1 [unit = u"K"]
+    n = length(Ci_star_standard)
+    @variables Ci_star[1:n] [description = "saturation concentrations"]
+
+    @parameters C_OA_unit = 1 [unit = u"m^3/μg", description = "make C_OA unitless"]
+    @variables C_OA [unit = u"μg/(m^3)", description = "organic compound in aerosol phase"]
+    @variables ξ[1:n] [description = "partitioning coefficient at corresponding Ci_star"]
+
+    eqs = [
+        [ξ[i] ~ 1/(1+calc_Ci_star(T/T_unit,Ci_star_standard[i])/(C_OA*C_OA_unit)) for i ∈ 1:n]
+        C_OA*C_OA_unit ~ sum(ξ.*Ci)
+    ]
+
+    NonlinearSystem(eqs, [C_OA, collect(ξ)...], [T,C_OA_unit,T_unit]; name=:VBS)
+end

test/VBS_test.jl

@@ -0,0 +1,26 @@
+using Aerosol
+using NonlinearSolve, ModelingToolkit
+
+@testset "VBS_typical" begin
+    Ci = [2.5, 1.8, 4.0, 4.0, 5.8, 4.8, 6.3, 8.0]
+    ns = VBS(Ci)
+    simplens = structural_simplify(ns)
+    @unpack C_OA = simplens
+    guess = [29]
+    params = []
+    prob = NonlinearProblem(simplens, guess, params)
+    sol = solve(prob, TrustRegion())
+    @test sol[C_OA] ≈ 10.6 atol=1e-2
+end
+
+@testset "different_temperature" begin
+    Ci = [2.5, 1.8, 4.0, 4.0, 5.8, 4.8, 6.3, 8.0]
+    ns = VBS(Ci)
+    simplens = structural_simplify(ns)
+    @unpack C_OA,T = simplens
+    guess = [29]
+    params = [T=>285]
+    prob = NonlinearProblem(simplens, guess, params)
+    sol = solve(prob, TrustRegion())
+    @test sol[C_OA] ≈ 15.93 atol=1e-2
+end

test/runtests.jl

@@ -4,4 +4,7 @@ using SafeTestsets
 
 @testset "Aerosol.jl" begin
     @safetestset "isorropia" begin include(joinpath(@__DIR__, "isorropia", "runtests.jl")) end
+    @safetestset "VBS" begin 
+        include("VBS_test.jl") 
+    end 
 end