From 17a1572fe25699ff9cb4c0f73e41acaf9c8b196e Mon Sep 17 00:00:00 2001 From: amylu00 Date: Mon, 8 Apr 2024 14:35:27 -0700 Subject: [PATCH] Linear HOM J option --- docs/src/API.md | 3 +- docs/src/IceNucleation.md | 8 ++++ docs/src/plots/HomFreezingPlots.jl | 8 ++-- docs/src/plots/linear_HOM_J.jl | 41 ++++++++++++++++++++ parcel/ParcelTendencies.jl | 2 +- src/IceNucleation.jl | 24 ++++++++++-- src/parameters/IceNucleation.jl | 6 +++ test/gpu_tests.jl | 30 ++++++++++++--- test/homogeneous_ice_nucleation_tests.jl | 49 ++++++++++++++++++++---- test/performance_tests.jl | 3 +- 10 files changed, 151 insertions(+), 23 deletions(-) create mode 100644 docs/src/plots/linear_HOM_J.jl diff --git a/docs/src/API.md b/docs/src/API.md index e8b224cd0..2c7d208eb 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -109,7 +109,8 @@ HetIceNucleation.INP_concentration_mean # Homogeneous ice nucleation ```@docs HomIceNucleation -HomIceNucleation.homogeneous_J +HomIceNucleation.homogeneous_J_cubic +HomIceNucleation.homogeneous_J_linear ``` # Common utility functions diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md index 9eeffd9b1..4dcec044e 100644 --- a/docs/src/IceNucleation.md +++ b/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 diff --git a/docs/src/plots/HomFreezingPlots.jl b/docs/src/plots/HomFreezingPlots.jl index e202631a6..aaac07e47 100644 --- a/docs/src/plots/HomFreezingPlots.jl +++ b/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 diff --git a/docs/src/plots/linear_HOM_J.jl b/docs/src/plots/linear_HOM_J.jl new file mode 100644 index 000000000..997c272c2 --- /dev/null +++ b/docs/src/plots/linear_HOM_J.jl @@ -0,0 +1,41 @@ +import CairoMakie as MK +import CloudMicrophysics as CM +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.HomIceNucleation as CMH + +FT = Float32 +const ip = CMP.IceNucleationParameters(FT) + +# 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) + +# 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, 600), fontsize = 18) +ax1 = MK.Axis( + fig[1, 1], + xlabel = "Δa_w [-]", + ylabel = "log10(J), J = [cm^-3 s^-1]", + title = "Linear fit to Koop2000", +) + +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 = :rb) + +#! format: on + +MK.save("linear_HOM_J.svg", fig) diff --git a/parcel/ParcelTendencies.jl b/parcel/ParcelTendencies.jl index 84d4630ab..245854be5 100644 --- a/parcel/ParcelTendencies.jl +++ b/parcel/ParcelTendencies.jl @@ -180,7 +180,7 @@ function homogeneous_freezing(params::ABHOM, PSD, state) Δa_w = ips.homogeneous.Δa_w_max end - J = CMI_hom.homogeneous_J(ips.homogeneous, Δa_w) + J = CMI_hom.homogeneous_J_cubic(ips.homogeneous, Δa_w) return min(max(FT(0), J * Nₗ * PSD.Vₗ), Nₗ / const_dt) end diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl index 10564880f..8a1ac7f0b 100644 --- a/src/IceNucleation.jl +++ b/src/IceNucleation.jl @@ -214,10 +214,11 @@ module HomIceNucleation import ..Parameters as CMP import Thermodynamics as TD -export homogeneous_J +export homogeneous_J_cubic +export homogeneous_J_linear """ - homogeneous_J(ip, Δa_w) + homogeneous_J_cubic(ip, Δa_w) - `ip` - a struct with ice nucleation parameters - `Δa_w` - change in water activity [-]. @@ -226,7 +227,7 @@ Returns the homogeneous freezing nucleation rate coefficient, `J`, in m^-3 s^-1 for sulphuric acid solutions. Parameterization based on Koop 2000, DOI: 10.1038/35020537. """ -function homogeneous_J(ip::CMP.Koop2000, Δa_w::FT) where {FT} +function homogeneous_J_cubic(ip::CMP.Koop2000, Δa_w::FT) where {FT} @assert Δa_w >= ip.Δa_w_min @assert Δa_w <= ip.Δa_w_max @@ -236,4 +237,21 @@ function homogeneous_J(ip::CMP.Koop2000, Δa_w::FT) where {FT} return FT(10)^(logJ) * 1e6 end +""" + homogeneous_J_linear(ip, Δa_w) + + - `ip` - a struct with ice nucleation parameters + - `Δa_w` - change in water activity [-]. + +Returns the homogeneous freezing nucleation rate coefficient, +`J`, in m^-3 s^-1 for sulphuric acid solutions. +Parameterization derived from a linear fit of the Koop 2000 parameterization, DOI: 10.1038/35020537. +""" +function homogeneous_J_linear(ip::CMP.Koop2000, Δa_w::FT) where {FT} + + logJ::FT = ip.linear_c₂ * Δa_w + ip.linear_c₁ + + return FT(10)^(logJ) * 1e6 +end + end # end module diff --git a/src/parameters/IceNucleation.jl b/src/parameters/IceNucleation.jl index 02a97572a..00052bac4 100644 --- a/src/parameters/IceNucleation.jl +++ b/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_coeff2 => :linear_c₂, ) parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics") FT = CP.float_type(td) diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index ae2e392c8..06b203d45 100644 --- a/test/gpu_tests.jl +++ b/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 @@ -1088,14 +1101,21 @@ function test_gpu(FT) (; 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_kernel!(backend, work_groups) + kernel! = + IceNucleation_homogeneous_J_cubic_kernel!(backend, work_groups) kernel!(ip, output, Delta_a_w; ndrange) - # test homogeneous_J is callable and returns a reasonable value + # test homogeneous_J_cubic is callable and returns a reasonable value @test Array(output)[1] ≈ FT(2.66194650334444e12) + + kernel! = + IceNucleation_homogeneous_J_linear_kernel!(backend, work_groups) + kernel!(ip, output, Delta_a_w; ndrange) + + # test homogeneous_J_linear is callable and returns a reasonable value + @test Array(output)[1] ≈ FT(7.156568123338207e11) end @testset "Modal nucleation kernels" begin diff --git a/test/homogeneous_ice_nucleation_tests.jl b/test/homogeneous_ice_nucleation_tests.jl index 4d578b665..2b7b6bdfc 100644 --- a/test/homogeneous_ice_nucleation_tests.jl +++ b/test/homogeneous_ice_nucleation_tests.jl @@ -10,13 +10,13 @@ import CloudMicrophysics.HomIceNucleation as CMH @info "Homogeneous Ice Nucleation Tests" -function test_homogeneous_J(FT) +function test_homogeneous_J_cubic(FT) tps = TD.Parameters.ThermodynamicsParameters(FT) H2SO4_prs = CMP.H2SO4SolutionParameters(FT) ip = CMP.IceNucleationParameters(FT) - TT.@testset "Homogeneous J" begin + TT.@testset "Homogeneous J (cubic parameterization)" begin T_warm = FT(229.2) T_cold = FT(228.8) @@ -25,30 +25,63 @@ function test_homogeneous_J(FT) Δa_w_too_large = FT(0.35) # higher nucleation rate at colder temperatures - TT.@test CMH.homogeneous_J( + TT.@test CMH.homogeneous_J_cubic( ip.homogeneous, CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold) - CO.a_w_ice(tps, T_cold), - ) > CMH.homogeneous_J( + ) > CMH.homogeneous_J_cubic( ip.homogeneous, CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm) - CO.a_w_ice(tps, T_warm), ) # If Δa_w out of range - TT.@test_throws AssertionError("Δa_w >= ip.Δa_w_min") CMH.homogeneous_J( + TT.@test_throws AssertionError("Δa_w >= ip.Δa_w_min") CMH.homogeneous_J_cubic( ip.homogeneous, Δa_w_too_small, ) - TT.@test_throws AssertionError("Δa_w <= ip.Δa_w_max") CMH.homogeneous_J( + TT.@test_throws AssertionError("Δa_w <= ip.Δa_w_max") CMH.homogeneous_J_cubic( ip.homogeneous, Δa_w_too_large, ) 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(Float64) +test_homogeneous_J_linear(Float64) println("Testing Float32") -test_homogeneous_J(Float32) +test_homogeneous_J_linear(Float32) diff --git a/test/performance_tests.jl b/test/performance_tests.jl index c0ef080c7..cadef0413 100644 --- a/test/performance_tests.jl +++ b/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_linear, (ip.homogeneous, Delta_a_w), 230) # non-equilibrium bench_press(CMN.τ_relax, (liquid,), 10)