Add "Flexible Microphysics" option, based on Cloudy.jl (#333)

* Add preliminary cloudy condensation + coalescence, with unit tests

* Add sedimentation

* Make cloudy deps an extension instead

* Add docs

* Add condensation-only parcel model example for Cloudy

* Add compat entry for Cloudy

* Update Cloudy version #

* Rename sedimentation + fix condensation tests

* Put units into SI

* Julia Formatter

* Add normalization factors

* Fix documentation errors

Project.toml

@@ -12,14 +12,17 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
 [weakdeps]
+Cloudy = "9e3b23bb-e7cc-4b94-886c-65de2234ba87"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 MLJ = "add582a8-e3ab-11e8-2d5e-e98b27df1bc7"
 
 [extensions]
+CloudyExt = "Cloudy"
 EmulatorModelsExt = ["DataFrames", "MLJ"]
 
 [compat]
 ClimaParams = "0.10"
+Cloudy = "0.4"
 DataFrames = "1.6"
 DocStringExtensions = "0.8, 0.9"
 ForwardDiff = "0.10"

docs/make.jl

@@ -37,6 +37,7 @@ Parameterizations = [
     "0-moment precipitation microphysics" => "Microphysics0M.md",
     "1-moment precipitation microphysics" => "Microphysics1M.md",
     "2-moment precipitation microphysics" => "Microphysics2M.md",
+    "N-moment precipitation microphysics" => "MicrophysicsFlexible.md",
     "P3 Scheme" => "P3Scheme.md",
     "Non-equilibrium cloud formation" => "MicrophysicsNonEq.md",
     "Smooth transition at thresholds" => "ThresholdsTransition.md",

docs/src/API.md

@@ -53,6 +53,15 @@ Microphysics2M.rain_evaporation
 Microphysics2M.conv_q_liq_to_q_rai
 ```
 
+# Flexible N-moment precipitation microphysics
+```@docs
+MicrophysicsFlexible
+MicrophysicsFlexible.CLSetup
+MicrophysicsFlexible.coalescence
+MicrophysicsFlexible.condensation
+MicrophysicsFlexible.weighted_vt
+```
+
 # P3 scheme
 
 ```@docs

docs/src/MicrophysicsFlexible.md

@@ -0,0 +1,41 @@
+# Microphysics Flexible
+
+The `MicrophysicsFlexible.jl` module relies on the extension defined in `ext/CloudyExt.jl`, based on a flexible N-moment microphysics scheme built in the external package [Cloudy.jl](https://github.com/CliMA/Cloudy.jl). This option currently handles warm-rain processes including coalescence, condensation/evaporation, and sedimentation (terminal velocity). Unlike typical moment-based schemes which distinguish between categories such as rain and cloud, and which determine rates of conversion between categories (the canonical autoconversion, accretion, and self-collection), this option gives the user the flexibility to define as many or as few moments as they please, with these coalescence-based processes being solved directly without relying on conversion rates. Likewise, the rate of condensation/evaporation is defined through the rate of diffusion of water vapor to/from the surface of droplets defined by the subdistributions which underpin the method. The user has not only the flexibility to specify the number of moments (and therefore the complexity/accuracy) to use, but also the assumed size distributions corresponding to these moments. For instance, one might define a 5-moment implementation using an Exponential mode for smaller cloud droplets, plus a Gamma mode for larger rain droplets. Or, more creatively, perhaps a 12-moment implementation comprised of four Gamma modes.
+
+Options for dynamics and size distributions are under continuous development in the `Cloudy.jl` package, thus only the default and suggested use cases are described in detail here.
+
+## Moments and Sub-Distributions
+
+The prognostic variables of this parameterization are a set of N moments, which can be further divided into P sets of moments, each of which correponds to a subdistribution p. By design these moments begin at order 0 and increase as integers up to the maximum number of parameters for the chosen subdistribution. The first three such default moments have interpretable meanings:
+  - ``M_0`` - the number density of droplets [1/m^3]
+  - ``M_1`` - the mass density of droplets [kg/m^3]
+  - ``M_2`` - proportional to the radar reflectivity [kg^2/m^3]
+and can be converted to more canonical definitions of `q_liq` and `q_rai` through numerical integration.
+
+When the user wishes to use more than 2 or 3 total variables to represent the system, these moments must be divided between ``P > 1`` sub-distributions, each of which assumes the form of a particular mathematical distribution, such as an Exponential, Lognormal, or Monodisperse (each of which has two parameters), or a Gamma distribution (which takes 3 parameters). 
+
+## Loading the extension
+The package `Cloudy.jl` and its dependencies are not loaded by default when using `CloudMicrophysics.jl`. Rather, one must specify:
+```
+using CloudMicrophysics
+using Cloudy
+```
+from the Julia REPL. Upon recognizing that `Cloudy.jl` is being loaded, the extension `CloudyExt.jl` will then be loaded and overwrite the function stubs defined in `src/MicrophysicsFlexible.jl`.
+
+## Setting up a system
+All the details from the number of moments and type of subdistributions, to the parameterizations of coalescence, condensation, and sedimentation are defined through the `CLSetup` (CLoudySetup) mutable struct. This struct is mutable specifically because certain of its components, such as backend-computed coalescence tendencies, are updated prior to being passed to the timestepper. The components of a `CLSetup` object and their defaults are further described below.
+
+|   component         |   description                              |   default                  |
+|---------------------|--------------------------------------------|----------------------------|
+| ``pdists``          | Vector of subdistributions corresponding   | ``[Exponential, Gamma]``   |
+|                     | to the moments                             |                            |
+| ``mom``             | Prognostic mass moments, in the same order |``[1e8 / m^3, 1e-2 kg/m^3,``|
+|                     | as the corresponding subdistributions;     |``1e6/m^3, 1e-3 kg/m^3,``   |
+|                     | first 2 for Exp, next 3 for Gamma          |``2e-12 kg^2/m^3]``         |
+| ``KernelFunc``      | Form of the coalescence kernel function    | ``LongKernelFunction``     |
+| ``mass_thresholds`` | Particle size thresholds for coalescence   | ``[10.0, Inf]``            |
+|                     | integration                                |                            |
+| ``kernel order``    | Polynomial order for the approx. kernel    | ``1``                      |
+| ``kernel_limit``    | Size threshold for approx. kernel          | ``500``                    |
+| ``vel``             | Power-series coefficients for velocity     | ``[2.0, 1/6]``             |
+| ``norms``           | Normalizing number density & mass          | ``[1e6/m^3, 1e-9 kg]``     |

ext/CloudyExt.jl

@@ -0,0 +1,158 @@
+"""
+Flexible N-moment microphysics representation, including: 
+  - Generalized collisional-coalescence described by a kernel K(x,y), with
+  default parameterization based on Long collision kernel 
+  - Power-series representation of fall-speed 
+  - Power-series representation of condensational growth
+  TODO: no representation of ventilation effects
+"""
+module CloudyExt
+
+import Thermodynamics as TD
+import Thermodynamics.Parameters as TDP
+
+import CloudMicrophysics.Common as CO
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.MicrophysicsFlexible:
+    CLSetup, coalescence, condensation, weighted_vt
+
+import Cloudy as CL
+import Cloudy.ParticleDistributions as CPD
+import Cloudy.KernelFunctions as CLK
+
+"""
+A structure containing the subdistributions, their moments, and additional
+dynamical parameters corresponding to rates of collision, sedimentation, and
+condensation/evaporation
+"""
+function CLSetup{FT}(;
+    pdists::Vector{<:CPD.PrimitiveParticleDistribution{FT}} = Vector([
+        CPD.ExponentialPrimitiveParticleDistribution(
+            FT(100 * 1e6),
+            FT(1e5 * 1e-18 * 1e3),
+        ), # 100/cm^3; 10^5 µm^3 = 1e-10 kg
+        CPD.GammaPrimitiveParticleDistribution(
+            FT(1 * 1e6),
+            FT(1e6 * 1e-18 * 1e3),
+            FT(1),
+        ),   # 1/cm^3; 10^6 µm^3 = 1e-9 kg; k=1
+    ]),
+    mom::Vector{FT} = FT.([100.0 * 1e6, 1e-2, 1.0 * 1e6, 1e-3, 2e-12]),
+    NProgMoms::Vector{Int} = [Integer(CPD.nparams(dist)) for dist in pdists],
+    KernelFunc::CLK.KernelFunction{FT} = CLK.LongKernelFunction(
+        FT(5.236e-10),  # 5.236e-10 kg;  
+        FT(9.44e9),     # 9.44e9 m^3/kg^2/s;
+        FT(5.78),       # 5.78 m^3/kg/s
+    ),
+    mass_thresholds::Vector{FT} = [FT(1e-9), FT(Inf)],
+    kernel_order::Int = 1,
+    kernel_limit::FT = FT(1 * 1e-9 * 1e3),   # ~1mm Diam ~ 1 mm^3 volume
+    coal_data = nothing,
+    vel::Vector{Tuple{FT, FT}} = [(FT(50.0), FT(1.0 / 6))], # 50 m/s/kg^(1/6),
+    norms::Vector{FT} = [1e6, 1e-9], # 1e6 / m^3; 1e-9 kg
+) where {FT <: AbstractFloat}
+
+    CLSetup{FT}(
+        pdists,
+        mom,
+        NProgMoms,
+        KernelFunc,
+        mass_thresholds,
+        kernel_order,
+        kernel_limit,
+        coal_data,
+        vel,
+        norms,
+    )
+end
+
+"""
+    coalescence(clinfo)
+
+ - `clinfo` - kwarg structure containing pdists, moments, and coalescence parameters
+ TODO: currently implemented only for analytical coalescence style
+
+Returns a vector of moment tendencies due to collisional coalescence
+"""
+function coalescence(clinfo::CLSetup{FT}) where {FT}
+    kernel_tensor = CL.KernelTensors.CoalescenceTensor(
+        clinfo.KernelFunc,
+        clinfo.kernel_order,
+        clinfo.kernel_limit,
+    )
+    if isnothing(clinfo.coal_data)
+        clinfo.coal_data = CL.Coalescence.initialize_coalescence_data(
+            CL.EquationTypes.AnalyticalCoalStyle(),
+            kernel_tensor,
+            clinfo.NProgMoms,
+            norms = clinfo.norms,
+            dist_thresholds = clinfo.mass_thresholds,
+        )
+    end
+    mom_norms =
+        CL.get_moments_normalizing_factors(clinfo.NProgMoms, clinfo.norms)
+    mom_normalized = clinfo.mom ./ mom_norms
+    # first, update the particle distributions
+    for (i, dist) in enumerate(clinfo.pdists)
+        ind_rng = CL.get_dist_moments_ind_range(clinfo.NProgMoms, i)
+        CPD.update_dist_from_moments!(dist, mom_normalized[ind_rng]) #clinfo.mom[ind_rng])
+    end
+    CL.Coalescence.update_coal_ints!(
+        CL.EquationTypes.AnalyticalCoalStyle(),
+        clinfo.pdists,
+        clinfo.coal_data,
+    )
+    return clinfo.coal_data.coal_ints .* mom_norms
+end
+
+"""
+    condensation(clinfo, aps, tps, q, ρ, T)
+
+ - `clinfo` - kwarg structure containing pdists, moments, and coalescence parameters
+ - `aps` - air properties
+ - `tps` - thermodynamics parameters
+ - `T` - air temperature
+ - `S` - saturation ratio (supersaturation = S - 1)
+Returns a vector of moment tendencies due to condensation/evaporation
+"""
+function condensation(
+    clinfo::CLSetup{FT},
+    aps::CMP.AirProperties{FT},
+    tps::TDP.ThermodynamicsParameters{FT},
+    T::FT,
+    S::FT,
+) where {FT}
+    ξ = CO.G_func(aps, tps, T, TD.Liquid())
+
+    mom_norms =
+        CL.get_moments_normalizing_factors(clinfo.NProgMoms, clinfo.norms)
+    mom_normalized = clinfo.mom ./ mom_norms
+    # first, update the particle distributions
+    for (i, dist) in enumerate(clinfo.pdists)
+        ind_rng = CL.get_dist_moments_ind_range(clinfo.NProgMoms, i)
+        CPD.update_dist_from_moments!(dist, mom_normalized[ind_rng])
+    end
+    return CL.Condensation.get_cond_evap(
+        S - 1,
+        (; ξ = ξ, pdists = clinfo.pdists),
+    ) .* mom_norms
+end
+
+"""
+    weighted_vt(clinfo)
+
+ - `clinfo` - kwarg structure containing pdists, moments, and coalescence parameters
+Returns the integrated fall speeds corresponding to the rate of change of prognostic moments
+"""
+function weighted_vt(clinfo::CLSetup{FT}) where {FT}
+    for (i, dist) in enumerate(clinfo.pdists)
+        ind_rng = CL.get_dist_moments_ind_range(clinfo.NProgMoms, i)
+        CPD.update_dist_from_moments!(dist, clinfo.mom[ind_rng])
+    end
+    sed_flux =
+        CL.Sedimentation.get_sedimentation_flux(clinfo.pdists, clinfo.vel)
+    return -sed_flux ./ clinfo.mom
+end
+
+
+end #module