Add condensation dynamic and simple test cases (#138)

* Add Condensation module

* Add unit tests

* Add simple condensation examples

* Add new examples to tests and buildkite

* CliMA Formatter

* Fix example and update version #

* Small changes, e.g. don't compute negative moments

* Get rid of array partition; update version number

* Increment subversion

* rm stray .x

* Remove array partition

* Update plotting helper

.buildkite/pipeline.yml

@@ -64,6 +64,14 @@ steps:
         command: "julia --color=yes --project=test test/examples/Analytical/box_gamma_mixture_long.jl"
         artifact_paths: "test/outputs/*"
 
+      - label: ":crystal_ball: Examples: analytical, box, single gamma, condensation"
+        command: "julia --color=yes --project=test test/examples/Analytical/condensation_single_gamma.jl"
+        artifact_paths: "test/outputs/*"
+
+      - label: ":crystal_ball: Examples: analytical, box, exp-gamma mixture, condensation"
+        command: "julia --color=yes --project=test test/examples/Analytical/condensation_exp_gamma.jl"
+        artifact_paths: "test/outputs/*"
+
       - label: ":crystal_ball: Examples: analytical, rainshaft, single gamma"
         command: "julia --color=yes --project=test test/examples/Analytical/rainshaft_single_gamma.jl"
         artifact_paths: "test/outputs/*"

Project.toml

@@ -1,7 +1,7 @@
 name = "Cloudy"
 uuid = "9e3b23bb-e7cc-4b94-886c-65de2234ba87"
 authors = ["Climate Modeling Alliance"]
-version = "0.3.1"
+version = "0.4.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

src/Cloudy.jl

@@ -6,6 +6,7 @@ include(joinpath("Kernels", "KernelFunctions.jl"))
 include(joinpath("ParticleDistributions", "ParticleDistributions.jl"))
 include(joinpath("Sources", "EquationTypes.jl"))
 include(joinpath("Sources", "Coalescence.jl"))
+include(joinpath("Sources", "Condensation.jl"))
 include(joinpath("Sources", "Sedimentation.jl"))
 
 end

src/Sources/Condensation.jl

@@ -0,0 +1,48 @@
+"""
+  multi-moment bulk microphysics implementation of condensation/evaporation
+
+  Includes only a single variations that involves analytical integration
+"""
+module Condensation
+
+using Cloudy
+using Cloudy.ParticleDistributions
+
+export get_cond_evap
+
+
+"""
+    get_cond_evap(s::FT, par::NamedTuple)
+
+    's' - supersaturation
+    `par` - ODE parameters, a NamedTuple containing a list of ParticleDistributions and the condensation coefficient ξ
+Returns the rate of change of all prognostic moments due to condensation and evaporation (without ventilation effects)
+based on the equation dg(x) / dt = -3ξs d/dx(x^{1/3} g(x))
+"""
+function get_cond_evap(s::FT, par::NamedTuple) where {FT <: Real}
+    ξ = par.ξ
+    n_dist = length(par.pdists)
+    NProgMoms = [nparams(dist) for dist in par.pdists]
+
+    # build diagnostic moments
+    mom = zeros(FT, sum(NProgMoms))
+    for i in 1:n_dist
+        for j in 2:NProgMoms[i]
+            ind = get_dist_moment_ind(NProgMoms, i, j)
+            mom[ind] = moment(par.pdists[i], FT(j - 1 - 2 / 3))
+        end
+    end
+
+    # calculate condensation/evaporation flux for prognostic moments
+    cond_evap_int = zeros(FT, sum(NProgMoms))
+    for i in 1:n_dist
+        for j in 2:NProgMoms[i]
+            ind = get_dist_moment_ind(NProgMoms, i, j)
+            cond_evap_int[ind] = 3 * FT(ξ) * s * (j - 1) * mom[ind]
+        end
+    end
+
+    return cond_evap_int
+end
+
+end #module Condensation.jl

src/Sources/Sedimentation.jl

@@ -11,10 +11,9 @@ export get_sedimentation_flux
 
 
 """
-  get_sedimentation_flux(mom_p::Array{Real}, par::Dict)
+  get_sedimentation_flux(par::NamedTuple)
 
-  `mom_p` - prognostic moments
-  `par` - ODE parameters, a NamedTuple containing a list of ParticleDistributions and terminal celocity coefficients.
+  `par` - ODE parameters, a NamedTuple containing a list of ParticleDistributions and terminal velocity coefficients.
 Returns sedimentation flux of all prognostic moments, which is the integral of terminal velocity times prognostic moments. The
 terminal velocity of particles is assumed to be expressed as: ∑ vel[i][1] * x^(vel[i][2]) where vel is a vector of tuples.
 """

test/examples/Analytical/condensation_exp_gamma.jl

@@ -0,0 +1,30 @@
+"Box model with a single gamma mode"
+
+using DifferentialEquations
+
+include("../utils/box_model_helpers.jl")
+include("../utils/plotting_helpers.jl")
+
+FT = Float64
+
+# Initial condition
+moments_init = [100.0, 10.0, 10.0, 10.0, 20.0]
+dist_init = [
+    ExponentialPrimitiveParticleDistribution(FT(100), FT(0.1)), # 100/cm^3; 10^5 µm^3
+    GammaPrimitiveParticleDistribution(FT(10), FT(1), FT(1)),   # 10/cm^3; 10^6 µm^3; k=1
+]
+
+# Solver
+s = 0.05
+ξ = 1e-2
+tspan = (FT(0), FT(120))
+rhs!(dm, m, par, t) = rhs_condensation!(dm, m, par, s)
+NProgMoms = [nparams(dist) for dist in dist_init]
+ODE_parameters = (; ξ = ξ, pdists = dist_init, NProgMoms = NProgMoms, dt = FT(10))
+prob = ODEProblem(rhs!, moments_init, tspan, ODE_parameters)
+sol = solve(prob, SSPRK33(), dt = ODE_parameters.dt)
+
+plot_params!(sol, ODE_parameters; file_name = "condensation_expgam_params.pdf")
+plot_moments!(sol, ODE_parameters; file_name = "condensation_expgam_moments.pdf")
+plot_spectra!(sol, ODE_parameters; file_name = "condensation_expgam_spectra.pdf", logxrange = (-3, 6))
+print_box_results!(sol, ODE_parameters)

test/examples/Analytical/condensation_single_gamma.jl

@@ -0,0 +1,27 @@
+"Box model with a single gamma mode"
+
+using DifferentialEquations
+
+include("../utils/box_model_helpers.jl")
+include("../utils/plotting_helpers.jl")
+
+FT = Float64
+
+# Initial condition
+moments_init = [100.0, 10.0, 2.0]
+dist_init = [GammaPrimitiveParticleDistribution(FT(100), FT(0.1), FT(1))]
+
+# Solver
+s = 0.05
+ξ = 1e-2
+tspan = (FT(0), FT(120))
+NProgMoms = [nparams(dist) for dist in dist_init]
+rhs!(dm, m, par, t) = rhs_condensation!(dm, m, par, s)
+ODE_parameters = (; ξ = ξ, pdists = dist_init, NProgMoms = NProgMoms, dt = FT(10))
+prob = ODEProblem(rhs!, moments_init, tspan, ODE_parameters)
+sol = solve(prob, SSPRK33(), dt = ODE_parameters.dt)
+
+plot_params!(sol, ODE_parameters; file_name = "condensation_single_gamma_params.pdf")
+plot_moments!(sol, ODE_parameters; file_name = "condensation_single_gamma_moments.pdf")
+plot_spectra!(sol, ODE_parameters; file_name = "condensation_single_gamma_spectra.pdf", logxrange = (-3, 6))
+print_box_results!(sol, ODE_parameters)

test/examples/utils/box_model_helpers.jl

@@ -6,6 +6,7 @@ using Cloudy.KernelFunctions
 using Cloudy.KernelTensors
 using Cloudy.ParticleDistributions
 using Cloudy.Coalescence
+using Cloudy.Condensation
 using Cloudy.EquationTypes
 
 const CPD = Cloudy.ParticleDistributions
@@ -30,6 +31,14 @@ function rhs_coal!(coal_type::CoalescenceStyle, dmom, mom, p)
     dmom .= p.coal_data.coal_ints
 end
 
+function rhs_condensation!(dmom, mom, p, s)
+    for (i, dist) in enumerate(p.pdists)
+        ind_rng = get_dist_moments_ind_range(p.NProgMoms, i)
+        update_dist_from_moments!(dist, mom[ind_rng])
+    end
+    dmom .= get_cond_evap(s, p)
+end
+
 """
   golovin_analytical_solution(x, x0, t; b = 1.5e-3, n = 1)
 

test/examples/utils/plotting_helpers.jl

@@ -195,9 +195,8 @@ function print_box_results!(sol, p)
     moments = vcat(sol.u'...)
 
     Ndist = length(p.pdists)
-    n_params = [nparams(p.pdists[i]) for i in 1:Ndist]
-    Nmom_min = minimum(n_params)
-    Nmom_max = maximum(n_params)
+    Nmom_min = minimum(p.NProgMoms)
+    Nmom_max = maximum(p.NProgMoms)
     moments_sum = zeros(length(time), Nmom_min)
 
     for i in 1:Ndist
@@ -208,18 +207,17 @@ function print_box_results!(sol, p)
         end
     end
     @show time
-    @show moments_sum[:, 1]
-    @show moments_sum[:, 2]
-    @show moments_sum[:, 3]
+    for j in 1:Nmom_min
+        @show moments_sum[:, j]
+    end
 
     t_ind = [1, ceil(Int, length(sol.t) / 2), length(sol.t)]
-    n_params = [nparams(p.pdists[i]) for i in 1:Ndist]
-    params = zeros(length(t_ind), Ndist, n_params[1])
+    params = zeros(length(t_ind), Ndist, Nmom_max)
     for i in 1:3
         for j in 1:Ndist
             ind_rng = get_dist_moments_ind_range(p.NProgMoms, j)
             update_dist_from_moments!(p.pdists[j], moments[t_ind[i], ind_rng])
-            params[i, j, :] = vcat(CPD.get_params(p.pdists[j])[2]...)
+            params[i, j, 1:p.NProgMoms[j]] = vcat(CPD.get_params(p.pdists[j])[2]...)
         end
     end
     @show t_ind

test/runtests.jl

@@ -11,6 +11,8 @@ include("./examples/Analytical/box_lognorm_mixture.jl")
 include("./examples/Analytical/box_mono_gamma_mixture.jl")
 include("./examples/Analytical/box_gamma_mixture_hydro.jl")
 include("./examples/Analytical/box_gamma_mixture_long.jl")
+include("./examples/Analytical/condensation_single_gamma.jl")
+include("./examples/Analytical/condensation_exp_gamma.jl")
 include("./examples/Analytical/rainshaft_single_gamma.jl")
 include("./examples/Analytical/rainshaft_gamma_mixture.jl")
 include("./examples/Numerical/single_particle_exp.jl")

test/unit_tests/test_Sources_correctness.jl

@@ -19,6 +19,7 @@ using Cloudy.Coalescence:
     initialize_coalescence_data,
     update_coal_ints!
 using Cloudy.Sedimentation
+using Cloudy.Condensation
 using Cloudy.KernelTensors
 using Cloudy.KernelFunctions
 using SpecialFunctions: gamma, gamma_inc
@@ -254,5 +255,21 @@ end
 ## Sedimentation.jl
 # Sedimentation moment flux tests
 par = (; pdists = [ExponentialPrimitiveParticleDistribution(1.0, 1.0)], vel = [(1.0, 0.0), (-1.0, 1.0 / 6)])
-mom = [1.0, 1.0]
 @test get_sedimentation_flux(par) ≈ [-1.0 + gamma(1.0 + 1.0 / 6), -1.0 + gamma(2.0 + 1.0 / 6)] rtol = rtol
+
+## Condensation.jl
+# Condensation moment tests
+par = (; pdists = [ExponentialPrimitiveParticleDistribution(1.0, 1.0)], ξ = 1e-6)
+@test get_cond_evap(0.01, par) ≈ [0.0, 3 * 1e-6 * 0.01 * moment(par.pdists[1], 1 - 2 / 3)] rtol = rtol
+
+par = (;
+    pdists = [ExponentialPrimitiveParticleDistribution(1.0, 1.0), GammaPrimitiveParticleDistribution(1.0, 2.0, 3.0)],
+    ξ = 1e-6,
+)
+@test get_cond_evap(0.01, par) ≈ [
+    0.0,
+    3 * 1e-6 * 0.01 * moment(par.pdists[1], 1 - 2 / 3),
+    0.0,
+    3 * 1e-6 * 0.01 * moment(par.pdists[2], 1 - 2 / 3),
+    3 * 2 * 1e-6 * 0.01 * moment(par.pdists[2], 2 - 2 / 3),
+] rtol = rtol

test/unit_tests/test_Sources_opt.jl

@@ -15,6 +15,8 @@ using Cloudy.Coalescence:
     initialize_coalescence_data,
     get_coalescence_integral_moment_qrs!,
     update_coal_ints!
+using Cloudy.Sedimentation
+using Cloudy.Condensation
 using JET: @test_opt
 using QuadGK
 
@@ -84,3 +86,13 @@ for pdists in ([dist1b], [dist1b, dist2b])
     @test_opt get_coalescence_integral_moment_qrs!(NumericalCoalStyle(), moment_order, pdists, cd)
     @test_opt update_coal_ints!(NumericalCoalStyle(), pdists, cd)
 end
+
+## Sedimentation.jl
+# Sedimentation moment flux tests
+par = (; pdists = [ExponentialPrimitiveParticleDistribution(1.0, 1.0)], vel = [(1.0, 0.0), (-1.0, 1.0 / 6)])
+@test_opt get_sedimentation_flux(par)
+
+## Condensation.jl
+# Condensation moment tests
+par = (; pdists = [ExponentialPrimitiveParticleDistribution(1.0, 1.0)], ξ = 1e-6)
+@test_opt get_cond_evap(0.01, par)