Merge #134

134: Moving H2SO4_sat_vap_pressure & Delta_a_w r=amylu00 a=amylu00 Co-authored-by: amylu00 <[email protected]>
CliMA · Jun 29, 2023 · 51778ba · 51778ba
2 parents d27e1ac + fbd45d7
commit 51778ba
Show file tree

Hide file tree

Showing 8 changed files with 154 additions and 152 deletions.
diff --git a/docs/src/API.md b/docs/src/API.md
@@ -78,8 +78,6 @@ AerosolActivation.total_M_activated
 ```@docs
 HetIceNucleation
 HetIceNucleation.dust_activated_number_fraction
-HetIceNucleation.H2SO4_soln_saturation_vapor_pressure
-HetIceNucleation.ABIFM_Delta_a_w
 HetIceNucleation.ABIFM_J
 ```
 
@@ -90,6 +88,8 @@ Common
 Common.G_func
 Common.logistic_function
 Common.logistic_function_integral
+Common.H2SO4_soln_saturation_vapor_pressure
+Common.Delta_a_w
 ```
 
 # Common utility types

diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md
@@ -18,7 +18,6 @@ The parameterization for deposition on dust particles is an implementation of
     heterogeneous immersion freezing or homogeneous freezing
     and modeling the competition between them.
 
-
 ## Activated fraction for deposition freezing on dust
 The parameterization models the activated fraction
   as an empirical function of ice saturation ratio,
@@ -55,7 +54,6 @@ Water Activity-Based Immersion Freezing Model (ABFIM)
   aerosols are assumed to contain an insoluble and soluble material. When immersed in water,
   the soluble material diffuses into the liquid water to create a sulphuric acid solution.
 
-
 Using empirical coefficients, ``m`` and ``c``, from [KnopfAlpert2013](@cite), 
   the heterogeneous nucleation rate coefficient in units of ``cm^{-2}s^{-1}`` can be determined by the linear equation
 ```math
@@ -82,7 +80,7 @@ Using empirical coefficients, ``m`` and ``c``, from [KnopfAlpert2013](@cite),
 where ``p_{sol}`` is saturated vapor pressure of water above solution, ``p_{sat}``
   is saturated vapor pressure above pure liquid water, and ``p_{i,sat}`` is saturated
   vapor pressure above ice. ``p_{sol}`` is determined in mbar using a parameterization
-  for supercooled, binary ``H_2SO_4/H_2O`` solution from [Luo1995](@cite) which is valid for 185K < T < 235K:
+  for supercooled, binary ``H_2SO_4/H_2O`` solution from [Luo1995](@cite) which is valid for ``185K < T < 235K``:
 ```math
 \begin{equation}
   ln(p_{sol}) = 23.306 - 5.3465x + 12xw_h - 8.19xw_h^2 + \frac{1}{T}(-5814 + 928.9x - 1876.7xw_h)
@@ -98,7 +96,8 @@ Once ``J_{het}`` is calculated, it can be used to determine the ice production r
   P_{ice} = J_{het}A(N_{tot}-N_{ice})
 \end{equation}
 ```
-where ``A`` is surface area of an individual ice nuclei, ``N_{tot}`` is total number of ice nuclei, and ``N_{ice}`` is number of ice crystals already in the system. 
+where ``A`` is surface area of an individual ice nuclei, ``N_{tot}`` is total number 
+  of ice nuclei, and ``N_{ice}`` is number of ice crystals already in the system. 
 
 ## ABIFM Example Figures
 ```@example
@@ -112,6 +111,7 @@ const PL = Plots
 const IN = CloudMicrophysics.HetIceNucleation
 const CMP = CloudMicrophysics.Parameters
 const CT = CloudMicrophysics.CommonTypes
+const CO = CloudMicrophysics.Common
 const CP =  CLIMAParameters
 const TD = Thermodynamics
 
@@ -133,7 +133,7 @@ dust_type = CT.KaoliniteType()
 
 it = 1
 for T in temp
-        Delta_a[it] = IN.ABIFM_Delta_a_w(prs, x, T)
+        Delta_a[it] = CO.Delta_a_w(prs, x, T)
         J[it] = IN.ABIFM_J(dust_type, Delta_a[it])
         global it += 1
 end
@@ -153,4 +153,4 @@ PL.plot!(KA13_Delta_a_param, KA13_log10J_param, linecolor = :red, label="paper p
 
 PL.savefig("Knopf_Alpert_fig_1.svg")
 ```
-![](Knopf_Alpert_fig_1.svg)
+![](Knopf_Alpert_fig_1.svg)
diff --git a/src/Common.jl b/src/Common.jl
@@ -11,6 +11,8 @@ const CMP = Parameters
 const APS = Parameters.AbstractCloudMicrophysicsParameters
 
 export G_func
+export H2SO4_soln_saturation_vapor_pressure
+export ABIFM_Delta_a_w
 
 """
     G_func(param_set, T, Liquid())
@@ -116,4 +118,52 @@ function logistic_function_integral(x::FT, x_0::FT, k::FT) where {FT <: Real}
     return _result
 end
 
+"""
+    H2SO4_soln_saturation_vapor_pressure(x, T)
+
+ - `x` - wt percent sulphuric acid [unitless] 
+ - `T` - air temperature [K].
+
+Returns the saturation vapor pressure above a sulphuric acid solution droplet in Pa.
+`x` is, for example, 0.1 if droplets are 10 percent sulphuric acid by weight
+"""
+function H2SO4_soln_saturation_vapor_pressure(x::FT, T::FT) where {FT <: Real}
+
+    @assert T < FT(235)
+    @assert T > FT(185)
+
+    w_h = 1.4408 * x
+    p_sol =
+        exp(
+            23.306 - 5.3465 * x + 12 * x * w_h - 8.19 * x * w_h^2 +
+            (-5814 + 928.9 * x - 1876.7 * x * w_h) / T,
+        ) * 100 # * 100 converts mbar --> Pa
+    return p_sol
+end
+
+"""
+    Delta_a_w(prs, x, T)
+
+ - `prs` - set with model parameters
+ - `x` - wt percent sulphuric acid [unitless]
+ - `T` - air temperature [K].
+
+Returns the change in water activity when droplet undergoes immersion freezing.
+`x` is, for example, 0.1 if droplets are 10 percent sulphuric acid by weight.
+"""
+function Delta_a_w(prs::APS, x::FT, T::FT) where {FT <: Real}
+
+    thermo_params = CMP.thermodynamics_params(prs)
+
+    p_sol = H2SO4_soln_saturation_vapor_pressure(x, T)
+    p_sat = TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid())
+    p_ice = TD.saturation_vapor_pressure(thermo_params, T, TD.Ice())
+
+    a_w = p_sol / p_sat
+    a_w_ice = p_ice / p_sat
+    Δa_w = a_w - a_w_ice
+
+    return min(Δa_w, FT(1))
+end
+
 end
diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl
@@ -10,8 +10,6 @@ import Thermodynamics as TD
 const APS = CMP.AbstractCloudMicrophysicsParameters
 
 export dust_activated_number_fraction
-export H2SO4_soln_saturation_vapor_pressure
-export ABIFM_Delta_a_w
 export ABIFM_J
 
 S0_warm(prs::APS, ::CT.ArizonaTestDustType) = CMP.S0_warm_ATD_Mohler2006(prs)
@@ -67,74 +65,22 @@ function dust_activated_number_fraction(
     end
 end
 
-"""
-    H2SO4_soln_saturation_vapor_pressure(x, T)
-
- - `x` - wt percent sulphuric acid [unitless] 
- - `T` - air temperature [K].
-
-Returns the saturation vapor pressure above a sulphuric acid solution droplet in Pa.
-`x` is, for example, 0.1 if droplets are 10 percent sulphuric acid by weight
-"""
-function H2SO4_soln_saturation_vapor_pressure(x::FT, T::FT) where {FT <: Real}
-
-    @assert T < FT(235)
-    @assert T > FT(185)
-
-    w_h = 1.4408 * x
-    p_sol =
-        exp(
-            23.306 - 5.3465 * x + 12 * x * w_h - 8.19 * x * w_h^2 +
-            (-5814 + 928.9 * x - 1876.7 * x * w_h) / T,
-        ) * 100 # * 100 converts mbar --> Pa
-    return p_sol
-end
-
-"""
-    ABIFM_Delta_a_w(prs, x, T)
-
- - `prs` - set with model parameters
- - `x` - wt percent sulphuric acid [unitless]
- - `T` - air temperature [K].
-
-Returns the change in water activity when droplet undergoes immersion freezing.
-`x` is, for example, 0.1 if droplets are 10 percent sulphuric acid by weight.
-To be used in conjuction with ABIFM_J to obtain nucleation rate coefficient
-"""
-function ABIFM_Delta_a_w(prs::APS, x::FT, T::FT) where {FT <: Real}
-
-    thermo_params = CMP.thermodynamics_params(prs)
-
-    p_sol = H2SO4_soln_saturation_vapor_pressure(x, T)
-    p_sat = TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid())
-    p_ice = TD.saturation_vapor_pressure(thermo_params, T, TD.Ice())
-
-    a_w = p_sol / p_sat
-    a_w_ice = p_ice / p_sat
-    Delta_a_w = a_w - a_w_ice
-
-    return min(Delta_a_w, FT(1))
-end
-
 """
     ABIFM_J(dust_type, Delta_a_w)
 
  - `dust_type` - choosing aerosol type
  - `Delta_a_w` - change in water activity [unitless].
 
 Returns the immersion freezing nucleation rate coefficient, `J`, in m^-2 s^-1 for sulphuric acid containing solutions. 
-For other solutions, p_sol should be adjusted accordingly. Delta_a_w can be found using ABIFM_Delta_a_w. 
+For other solutions, p_sol should be adjusted accordingly. Delta_a_w can be found using the Delta_a_w function in Common.jl. 
 `m` and `c` constants are taken from Knopf & Alpert 2013.
 """
-function ABIFM_J(
-    dust_type::CT.AbstractAerosolType,
-    Delta_a_w::FT,
-) where {FT <: Real}
+function ABIFM_J(dust_type::CT.AbstractAerosolType, Δa_w::FT) where {FT <: Real}
 
     m = J_het_coeff_m(dust_type)
     c = J_het_coeff_c(dust_type)
 
-    logJ = m * Delta_a_w + c
+    logJ = m * Δa_w + c
 
     return max(0, 10^logJ * 100^2) # converts cm^-2 s^-1 to m^-2 s^-1
 end

diff --git a/test/common_functions_tests.jl b/test/common_functions_tests.jl
@@ -1,8 +1,14 @@
-import Test
+import Test as TT
 
 import CloudMicrophysics
 const CO = CloudMicrophysics.Common
 
+import Thermodynamics as TD
+import CLIMAParameters as CP
+
+include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
+
+@info "Common Functions Tests"
 
 TT.@testset "logistic_function unit tests" begin
 
@@ -43,3 +49,63 @@ TT.@testset "logistic_function_integral unit tests" begin
     )
 
 end
+
+function test_H2SO4_soln_saturation_vapor_pressure(FT)
+    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+    prs = cloud_microphysics_parameters(toml_dict)
+
+    TT.@testset "H2SO4 solution saturated vapor pressure" begin
+
+        T_warm = FT(225.0)
+        T_cold = FT(200.0)
+        T_too_warm = FT(240)
+        T_too_cold = FT(180)
+        x_sulph = FT(0.1)
+
+        # If T out of range
+        TT.@test_throws AssertionError("T < FT(235)") CO.H2SO4_soln_saturation_vapor_pressure(
+            x_sulph,
+            T_too_warm,
+        )
+        TT.@test_throws AssertionError("T > FT(185)") CO.H2SO4_soln_saturation_vapor_pressure(
+            x_sulph,
+            T_too_cold,
+        )
+
+        # p_sol should be higher at warmer temperatures
+        TT.@test CO.H2SO4_soln_saturation_vapor_pressure(x_sulph, T_warm) >
+                 CO.H2SO4_soln_saturation_vapor_pressure(x_sulph, T_cold)
+    end
+end
+
+function test_Delta_a_w(FT)
+    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+    prs = cloud_microphysics_parameters(toml_dict)
+
+    TT.@testset "Delta_a_w" begin
+
+        T_warm = FT(229.2)
+        T_cold = FT(228.8)
+        x_sulph = FT(0.1)
+
+        # Delta_a_w never greater than 1
+        for T in [T_warm, T_cold]
+            TT.@test CO.Delta_a_w(prs, x_sulph, T) <= FT(1)
+        end
+    end
+end
+
+println("Testing Float64")
+test_H2SO4_soln_saturation_vapor_pressure(Float64)
+
+
+println("Testing Float32")
+test_H2SO4_soln_saturation_vapor_pressure(Float32)
+
+
+println("Testing Float64")
+test_Delta_a_w(Float64)
+
+
+println("Testing Float32")
+test_Delta_a_w(Float32)
diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl
@@ -14,6 +14,7 @@ include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
 const AM = CM.AerosolModel
 const AA = CM.AerosolActivation
 const CMI = CM.HetIceNucleation
+const CO = CM.Common
 const CMT = CM.CommonTypes
 const CM0 = CM.Microphysics0M
 const CM1 = CM.Microphysics1M
@@ -267,7 +268,7 @@ end
     end
 end
 
-@kernel function test_IceNucleation_H2SO4_soln_saturation_vapor_pressure_kernel!(
+@kernel function test_Common_H2SO4_soln_saturation_vapor_pressure_kernel!(
     output::AbstractArray{FT},
     x_sulph,
     T,
@@ -276,11 +277,11 @@ end
     i = @index(Group, Linear)
 
     @inbounds begin
-        output[i] = CMI.H2SO4_soln_saturation_vapor_pressure(x_sulph[i], T[i])
+        output[i] = CO.H2SO4_soln_saturation_vapor_pressure(x_sulph[i], T[i])
     end
 end
 
-@kernel function test_IceNucleation_ABIFM_Delta_a_w_kernel!(
+@kernel function test_Common_Delta_a_w_kernel!(
     prs,
     output::AbstractArray{FT},
     x_sulph,
@@ -290,7 +291,7 @@ end
     i = @index(Group, Linear)
 
     @inbounds begin
-        output[i] = CMI.ABIFM_Delta_a_w(prs, x_sulph[i], T[i])
+        output[i] = CO.Delta_a_w(prs, x_sulph[i], T[i])
     end
 end
 
@@ -552,7 +553,7 @@ function test_gpu(FT)
         @test Array(output)[3] ≈ FT(0)
     end
 
-    @testset "Ice Nucleation kernels" begin
+    @testset "Common Kernels" begin
         data_length = 1
         output = ArrayType(Array{FT}(undef, 1, data_length))
         fill!(output, FT(-44))
@@ -564,11 +565,10 @@ function test_gpu(FT)
         T = ArrayType([FT(230)])
         x_sulph = ArrayType([FT(0.1)])
 
-        kernel! =
-            test_IceNucleation_H2SO4_soln_saturation_vapor_pressure_kernel!(
-                dev,
-                work_groups,
-            )
+        kernel! = test_Common_H2SO4_soln_saturation_vapor_pressure_kernel!(
+            dev,
+            work_groups,
+        )
         event = kernel!(output, x_sulph, T, ndrange = ndrange)
         wait(dev, event)
 
@@ -586,13 +586,15 @@ function test_gpu(FT)
         T = ArrayType([FT(230)])
         x_sulph = ArrayType([FT(0.1)])
 
-        kernel! = test_IceNucleation_ABIFM_Delta_a_w_kernel!(dev, work_groups)
+        kernel! = test_Common_Delta_a_w_kernel!(dev, work_groups)
         event = kernel!(make_prs(FT), output, x_sulph, T, ndrange = ndrange)
         wait(dev, event)
 
-        # test if ABIFM_Delta_a_w is callable and returns reasonable values
+        # test if Delta_a_w is callable and returns reasonable values
         @test Array(output)[1] ≈ FT(0.2750536615)
+    end
 
+    @testset "Ice Nucleation kernels" begin
         data_length = 2
         output = ArrayType(Array{FT}(undef, 1, data_length))
         fill!(output, FT(-44.0))