create homogeneous_linear, add tests

docs/src/API.md

@@ -108,7 +108,8 @@ HetIceNucleation.INP_concentration_mean
 # Homogeneous ice nucleation
 ```@docs
 HomIceNucleation
-HomIceNucleation.homogeneous_J
+HomIceNucleation.homogeneous_J_cubic
+HomIceNucleation.homogeneous_J_linear
 ```
 
 # Common utility functions

docs/src/IceNucleation.md

@@ -243,6 +243,14 @@ Multiple sulphuric acid concentrations, ``x``,
     values that are obtained from pure water droplets. Though CliMA lines look far
     from the Spichtinger 2023 line, the lines seem to move closer as x approaches 0.
 
+### Linear fit homogeneous freezing
+To avoid complications in extrapolating the Koop 2000 parameterization outside of the valide temperature range, a linear fit was created.
+
+```@example
+include("plots/linear_HOM_J.jl")
+```
+![](linear_HOM_J.svg)
+
 ## Water Activity Based vs P3 Ice Nucleation Parameterizations
 The water activity based models are compared with P3 ice nucleation parameterizations as described
   in [MorrisonMilbrandt2015](@cite) using an adiabatic parcel model with depositional growth. The

docs/src/plots/HomFreezingPlots.jl

@@ -31,16 +31,16 @@ x_sulph = Vector{FT}([0.03, 0.04, 0.06])           # wt% sulphuric acid in dropl
     T in T_range
 ]
 
-J1 = @. CMI.homogeneous_J(ip.homogeneous, Δa1)
-J2 = @. CMI.homogeneous_J(ip.homogeneous, Δa2)
-J3 = @. CMI.homogeneous_J(ip.homogeneous, Δa3)
+J1 = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa1)
+J2 = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa2)
+J3 = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa3)
 
 log10J_1 = @. log10(J1)
 log10J_2 = @. log10(J2)
 log10J_3 = @. log10(J3)
 
 Δa_range = range(0.27, stop = 0.32, length = 50)
-J_given_Δa = @. CMI.homogeneous_J(ip.homogeneous, Δa_range)
+J_given_Δa = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa_range)
 
 #! format: off
 # Spichtinger et al 2023

docs/src/plots/linear_HOM_J.jl

@@ -1,64 +1,33 @@
-import OrdinaryDiffEq as ODE
 import CairoMakie as MK
-import Thermodynamics as TD
 import CloudMicrophysics as CM
-import ClimaParams as CP
-import Distributions as DS
 import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.HomIceNucleation as CMH
 
-const CMO = CM.Common
-FT = Float32
-const tps = TD.Parameters.ThermodynamicsParameters(FT)
+const ip = CMP.IceNucleationParameters(FT)
 
-# From Atkinson, Murray, and O'Sullivan 2016
-# (DOI:10.1021/acs.jpca.6b03843)
-# Using overall volume fit
-function Atkinson_J(T)
-    lnJ = -4.0106 * T + 963.706
-    return exp(lnJ)
-end
+# Koop2000 parameterization
+Koop_Δa = collect(0.26:0.0025:0.34)
+KoopJ = @. CMH.homogeneous_J_cubic(ip.homogeneous, Koop_Δa) # SI units
+log10KoopJ = @. log10(KoopJ .* 1e-6)
 
-# Initializing
-T_range = range(FT(234.0), stop = 237, length = 50)  # air temperature
-Sₗ = FT(1.1)
-eₛ = [TD.saturation_vapor_pressure(tps, T, TD.Liquid()) for T in T_range]
-e_range = Sₗ .* eₛ
-
-Δa = zeros(FT, length(T_range))
-
-# Solving for Δa and J values
-J = [Atkinson_J(T) for T in T_range] # [cm^-3 s^-1]
-log10J = @. log10(J)
-for (i, T) in enumerate(T_range)
-    Δa[i] = CMO.a_w_eT(tps, e_range[i], T) - CMO.a_w_ice(tps, T)
-end
-
-# Obtaining a linear parameterization for J in terms of Δa
-slope = (log10J[end] - log10J[1]) / (Δa[end] - Δa[1])
-intercept = slope * (Δa[1]) + log10J[1]
+# Linear fit on Koop 2000
+M = [ones(length(Koop_Δa)) Koop_Δa]
+linear_coeffs = M \ log10KoopJ
+new_log10J = [linear_coeffs[2] * Delta_a + linear_coeffs[1] for Delta_a in Koop_Δa]
 
 # Plotting J vs Δa
-fig = MK.Figure(resolution = (800, 700))
+fig = MK.Figure(resolution = (800, 600), fontsize = 18)
 ax1 = MK.Axis(
     fig[1, 1],
-    ylabel = "log10(J), J = [cm^-3 s^-1]",
-    xlabel = "Temperature [K]",
-    title = "log10(J) vs T",
-)
-ax2 = MK.Axis(
-    fig[2, 1],
     xlabel = "Δa_w [-]",
     ylabel = "log10(J), J = [cm^-3 s^-1]",
-    title = "log10(J) vs Δa_w",
+    title = "Linear fit to Koop2000",
 )
 
-MK.lines!(ax1, T_range, log10J, label = "ln(J) = -4.0106 T + 963.706",)
-MK.lines!(ax2, Δa, log10J,
-    label = "log10(J) = $slope Δa + $intercept",
-)
+MK.lines!(ax1, Koop_Δa, new_log10J, label = "log10(J) = $(linear_coeffs[2]) Δa + $(linear_coeffs[1])")
+MK.lines!(ax1, Koop_Δa, log10KoopJ, label = "Koop 2000 Parameterization")
 
-MK.axislegend(ax1, position = :rt, labelsize = 18.0f0)
-MK.axislegend(ax2, position = :rb, labelsize = 18.0f0)
+MK.axislegend(ax1, position = :rb)
 
 #! format: on
 

src/IceNucleation.jl

@@ -245,11 +245,11 @@ end
 
 Returns the homogeneous freezing nucleation rate coefficient,
 `J`, in m^-3 s^-1 for sulphuric acid solutions.
-Parameterization based on Atkinson et al 2016, DOI:10.1021/acs.jpca.6b03847.
+Parameterization derived from a linear fit of the Koop 2000 parameterization, DOI: 10.1038/35020537.
 """
-function homogeneous_J_linear(ip::CMP.Atkinson2016, Δa_w::FT) where {FT}
+function homogeneous_J_linear(ip::CMP.Koop2000, Δa_w::FT) where {FT}
 
-    logJ::FT = ip.c₁ * Δa_w + ip.c₂
+    logJ::FT = ip.linear_c₁ * Δa_w + ip.linear_c₂
 
     return FT(10)^(logJ) * 1e6
 end

src/parameters/IceNucleation.jl

@@ -49,6 +49,10 @@ Base.@kwdef struct Koop2000{FT} <: ParametersType{FT}
     c₃::FT
     "coefficient [-]"
     c₄::FT
+    "coefficient [-]"
+    linear_c₁::FT
+    "coefficient [-]"
+    linear_c₂::FT
 end
 
 function Koop2000(td::CP.AbstractTOMLDict)
@@ -59,6 +63,8 @@ function Koop2000(td::CP.AbstractTOMLDict)
         :Koop2000_J_hom_coeff2 => :c₂,
         :Koop2000_J_hom_coeff3 => :c₃,
         :Koop2000_J_hom_coeff4 => :c₄,
+        :Linear_J_hom_coeff1 => :linear_c₁,
+        :Linear_J_hom_coeff1 => :linear_c₂,
     )
     parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics")
     FT = CP.float_type(td)

test/gpu_tests.jl

@@ -629,7 +629,7 @@ end
     end
 end
 
-@kernel function IceNucleation_homogeneous_J_kernel!(
+@kernel function IceNucleation_homogeneous_J_cubic_kernel!(
     ip,
     output::AbstractArray{FT},
     Delta_a_w,
@@ -638,7 +638,20 @@ end
     i = @index(Group, Linear)
 
     @inbounds begin
-        output[1] = CMI_hom.homogeneous_J(ip.homogeneous, Delta_a_w[1])
+        output[1] = CMI_hom.homogeneous_J_cubic(ip.homogeneous, Delta_a_w[1])
+    end
+end
+
+@kernel function IceNucleation_homogeneous_J_linear_kernel!(
+    ip,
+    output::AbstractArray{FT},
+    Delta_a_w,
+) where {FT}
+
+    i = @index(Group, Linear)
+
+    @inbounds begin
+        output[1] = CMI_hom.homogeneous_J_linear(ip.homogeneous, Delta_a_w[1])
     end
 end
 
@@ -1091,10 +1104,23 @@ function test_gpu(FT)
         x_sulph = ArrayType([FT(0.15)])
         Delta_a_w = ArrayType([FT(0.2907389666103033)])
 
-        kernel! = IceNucleation_homogeneous_J_kernel!(backend, work_groups)
+        kernel! = IceNucleation_homogeneous_J_cubic_kernel!(backend, work_groups)
+        kernel!(ip, output, Delta_a_w; ndrange)
+
+        # test homogeneous_J_cubic is callable and returns a reasonable value
+        @test Array(output)[1] ≈ FT(2.66194650334444e12)
+
+        dims = (1, 1)
+        (; output, ndrange) = setup_output(dims, FT)
+
+        T = ArrayType([FT(220)])
+        x_sulph = ArrayType([FT(0.15)])
+        Delta_a_w = ArrayType([FT(0.2907389666103033)])
+
+        kernel! = IceNucleation_homogeneous_J_linear_kernel!(backend, work_groups)
         kernel!(ip, output, Delta_a_w; ndrange)
 
-        # test homogeneous_J is callable and returns a reasonable value
+        # test homogeneous_J_linear is callable and returns a reasonable value
         @test Array(output)[1] ≈ FT(2.66194650334444e12)
     end
 

test/homogeneous_ice_nucleation_tests.jl

@@ -47,8 +47,41 @@ function test_homogeneous_J_cubic(FT)
     end
 end
 
+function test_homogeneous_J_linear(FT)
+
+    tps = TD.Parameters.ThermodynamicsParameters(FT)
+    H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
+    ip = CMP.IceNucleationParameters(FT)
+
+    TT.@testset "Homogeneous J (linear parameterization)" begin
+
+        T_warm = FT(229.2)
+        T_cold = FT(228.8)
+        x_sulph = FT(0.1)
+
+        # higher nucleation rate at colder temperatures
+        TT.@test CMH.homogeneous_J_linear(
+            ip.homogeneous,
+            CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold) -
+            CO.a_w_ice(tps, T_cold),
+        ) > CMH.homogeneous_J_linear(
+            ip.homogeneous,
+            CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm) -
+            CO.a_w_ice(tps, T_warm),
+        )
+
+    end
+end
+
+
 println("Testing Float64")
 test_homogeneous_J_cubic(Float64)
 
 println("Testing Float32")
 test_homogeneous_J_cubic(Float32)
+
+println("Testing Float64")
+test_homogeneous_J_linear(Float64)
+
+println("Testing Float32")
+test_homogeneous_J_linear(Float32)

test/performance_tests.jl

@@ -159,7 +159,8 @@ function benchmark_test(FT)
         (ip_frostenberg, INPC, T_air_cold),
         150,
     )
-    bench_press(CMI_hom.homogeneous_J, (ip.homogeneous, Delta_a_w), 230)
+    bench_press(CMI_hom.homogeneous_J_cubic, (ip.homogeneous, Delta_a_w), 230)
+    bench_press(CMI_hom.homogeneous_J_lienar, (ip.homogeneous, Delta_a_w), 230)
 
     # non-equilibrium
     bench_press(CMN.τ_relax, (liquid,), 10)