diff --git a/Project.toml b/Project.toml index 52ca40a09..249ac8029 100644 --- a/Project.toml +++ b/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" diff --git a/docs/make.jl b/docs/make.jl index f937a1a4d..876fccbe6 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/API.md b/docs/src/API.md index 6b7e209f9..3fbf1c802 100644 --- a/docs/src/API.md +++ b/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 diff --git a/docs/src/MicrophysicsFlexible.md b/docs/src/MicrophysicsFlexible.md new file mode 100644 index 000000000..28374db7d --- /dev/null +++ b/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]`` | \ No newline at end of file diff --git a/ext/CloudyExt.jl b/ext/CloudyExt.jl new file mode 100644 index 000000000..65aa4d572 --- /dev/null +++ b/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 diff --git a/parcel/Example_NMoment_condensation.jl b/parcel/Example_NMoment_condensation.jl new file mode 100644 index 000000000..8326d508a --- /dev/null +++ b/parcel/Example_NMoment_condensation.jl @@ -0,0 +1,178 @@ +import OrdinaryDiffEq as ODE +import CairoMakie as MK +import Thermodynamics as TD +import CloudMicrophysics as CM +import Cloudy as CL +import Cloudy.ParticleDistributions as CPD +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.MicrophysicsFlexible as CMF +import ClimaParams as CP + +FT = Float64 + +""" + ODE problem definitions +""" +function parcel_model_cloudy(dY, Y, p, t) + # Numerical precision used in the simulation + FT = eltype(Y[5:end]) + # Simulation parameters + (; wps, tps, aps, w, clinfo) = p + # Y values + Sₗ = Y[1] + p_air = Y[2] + T = Y[3] + qᵥ = Y[4] + moments = Y[5:end] + + # Constants + Rᵥ = TD.Parameters.R_v(tps) + grav = TD.Parameters.grav(tps) + ρₗ = wps.ρw + + # Get thermodynamic parameters, phase partition and create thermo state. + q = TD.PhasePartition(qᵥ, FT(0), FT(0)) # use dry air density to compute ql + ts = TD.PhaseNonEquil_pTq(tps, p_air, T, q) + ρ_air = TD.air_density(tps, ts) + + qᵢ = FT(0.0) + qₗ = CPD.get_standard_N_q(clinfo.pdists).M_liq * ρₗ / ρ_air + q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) + ts = TD.PhaseNonEquil_pTq(tps, p_air, T, q) + + # Constants and variables that depend on the moisture content + R_air = TD.gas_constant_air(tps, q) + cp_air = TD.cp_m(tps, q) + L_subl = TD.latent_heat_sublim(tps, T) + L_fus = TD.latent_heat_fusion(tps, T) + L_vap = TD.latent_heat_vapor(tps, T) + ρ_air = TD.air_density(tps, ts) + + # Adiabatic parcel coefficients + a1 = L_vap * grav / cp_air / T^2 / Rᵥ - grav / R_air / T + a2 = 1 / qᵥ + a3 = L_vap^2 / Rᵥ / T^2 / cp_air + a4 = L_vap * L_subl / Rᵥ / T^2 / cp_air + a5 = L_vap * L_fus / Rᵥ / cp_air / (T^2) + + # Cloudy condnsational growth / evaporation + dmom_ce = CMF.condensation(clinfo, aps, tps, T, Sₗ) + @show dmom_ce + + # ... water mass budget + dqₗ_dt_v2l = dmom_ce[2] + dmom_ce[4] + + # Update the tendecies + dqₗ_dt = dqₗ_dt_v2l + dqᵥ_dt = -dqₗ_dt + + dSₗ_dt = a1 * w * Sₗ - (a2 + a3) * Sₗ * dqₗ_dt_v2l + + dp_air_dt = -p_air * grav / R_air / T * w + + dT_dt = -grav / cp_air * w + L_vap / cp_air * dqₗ_dt_v2l + + # Set tendencies + dY[1] = dSₗ_dt # saturation ratio over liquid water + dY[2] = dp_air_dt # pressure + dY[3] = dT_dt # temperature + dY[4] = dqᵥ_dt # vapor specific humidity + dY[5:end] = dmom_ce +end + +function run_parcel_cloudy(Yinit, clinfo, t_0, t_end, pp) + FT = typeof(t_0) + + println(" ") + println("Size distributions: ", clinfo.pdists) + println("Condensation growth only ") + + # Parameters for the ODE solver + p = (wps = pp.wps, aps = pp.aps, tps = pp.tps, w = pp.w, clinfo = clinfo) + + problem = + ODE.ODEProblem(parcel_model_cloudy, Yinit, (FT(t_0), FT(t_end)), p) + return ODE.solve( + problem, + ODE.Euler(), + dt = pp.const_dt, + reltol = 100 * eps(FT), + abstol = 100 * eps(FT), + ) +end + +# Get free parameters +tps = TD.Parameters.ThermodynamicsParameters(FT) +wps = CMP.WaterProperties(FT) +aps = CMP.AirProperties(FT) +# Constants +ρₗ = wps.ρw +ρᵢ = wps.ρi +R_v = TD.Parameters.R_v(tps) +R_d = TD.Parameters.R_d(tps) + +# Initial conditions +dist_init = [ + CPD.ExponentialPrimitiveParticleDistribution( + FT(100 * 1e6), + FT(1e5 * 1e-18 * 1e3), + ), # 100/cm^3; 10^5 µm^3 + CPD.GammaPrimitiveParticleDistribution( + FT(1 * 1e6), + FT(1e6 * 1e-18 * 1e3), + FT(1), + ), # 1/cm^3; 10^6 µm^3; k=1 +] +moments_init = FT.([100.0 * 1e6, 1e-2, 1.0 * 1e6, 1e-3, 2e-12]) +p₀ = FT(800 * 1e2) +T₀ = FT(273.15 + 7.0) +e = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) +Sₗ = FT(1) +md_v = (p₀ - e) / R_d / T₀ +mv_v = e / R_v / T₀ +ml_v = moments_init[2] + moments_init[4] +qᵥ = mv_v / (md_v + mv_v + ml_v) +qₗ = ml_v / (md_v + mv_v + ml_v) +qᵢ = FT(0) + +# Simulation parameters passed into ODE solver +w = FT(1) # updraft speed +const_dt = FT(0.5) # model timestep +t_max = FT(20) +clinfo = CMF.CLSetup{FT}(pdists = dist_init, mom = moments_init) + +Y0 = [Sₗ, p₀, T₀, qᵥ, moments_init...] +dY = zeros(FT, 9) +p = ( + wps = wps, + aps = aps, + tps = tps, + w = w, + clinfo = clinfo, + const_dt = const_dt, +) + +# Setup the plots +fig = MK.Figure(size = (800, 600)) +ax1 = MK.Axis(fig[1, 1], ylabel = "Supersaturation [%]") +ax2 = MK.Axis(fig[3, 1], xlabel = "Time [s]", ylabel = "Temperature [K]") +ax3 = MK.Axis(fig[2, 1], ylabel = "q_vap [g/kg]") +ax4 = MK.Axis(fig[2, 2], xlabel = "Time [s]", ylabel = "q_liq [g/kg]") +ax5 = MK.Axis(fig[1, 2], ylabel = "radius [μm]") + +# solve ODE +sol = run_parcel_cloudy(Y0, clinfo, FT(0), t_max, p) +@show sol.u + +# plot results - time series +MK.lines!(ax1, sol.t, (sol[1, :] .- 1) * 100.0) +MK.lines!(ax2, sol.t, sol[3, :]) +MK.lines!(ax3, sol.t, sol[4, :] * 1e3) +ρ_air = sol[2, :] / R_v ./ sol[3, :] +M_l = sol[6, :] + sol[8, :] # kg / m^3 air +N_l = sol[5, :] + sol[7, :] # number / m^3 air +r_l = (M_l ./ N_l / ρₗ / 4 / π * 3) .^ (1 / 3) * 1e6 +MK.lines!(ax4, sol.t, M_l ./ ρ_air * 1e3) +MK.lines!(ax5, sol.t, r_l) + +MK.save("cloudy_parcel.svg", fig) diff --git a/parcel/Project.toml b/parcel/Project.toml index 0129945f4..0d19bba40 100644 --- a/parcel/Project.toml +++ b/parcel/Project.toml @@ -1,7 +1,8 @@ [deps] -ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" +Cloudy = "9e3b23bb-e7cc-4b94-886c-65de2234ba87" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" diff --git a/src/CloudMicrophysics.jl b/src/CloudMicrophysics.jl index 16987404b..298bf26d0 100644 --- a/src/CloudMicrophysics.jl +++ b/src/CloudMicrophysics.jl @@ -21,5 +21,6 @@ include("AerosolActivation.jl") include("Nucleation.jl") include("IceNucleation.jl") include("PrecipitationSusceptibility.jl") +include("MicrophysicsFlexible.jl") end # module diff --git a/src/MicrophysicsFlexible.jl b/src/MicrophysicsFlexible.jl new file mode 100644 index 000000000..519622179 --- /dev/null +++ b/src/MicrophysicsFlexible.jl @@ -0,0 +1,73 @@ +""" +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 + TODO: conversion back to N_rai, N_liq, q_rai, q_liq +""" +module MicrophysicsFlexible + +""" +A structure containing the subdistributions, their moments, and additional +dynamical parameters corresponding to rates of collision, sedimentation, and +condensation/evaporation +""" +mutable struct CLSetup{FT} + "Subdistributions of the hydrometeor size distribution; defaults to + exponential cloud mode and gamma rain mode" + pdists::Vector + "Moments of the subdistributions; should correspond with pdists" + mom::Vector{FT} + "Total number of prognostic moments for each subdistribution" + NProgMoms::Vector{Int} + "Kernel function for collisional coalescence" + KernelFunc::Any + "Particle mass thresholds for analytical integration style" + mass_thresholds::Vector{FT} + "Polynomial order of the kernel approximation for coalescence" + kernel_order::Int + "Upper limit for evaluation of the polynomial kernel approximation" + kernel_limit::FT + "Coalescence data" + coal_data::Any + "Sedimentation rate parameters" + vel::Vector{Tuple{FT, FT}} + "Normalizing factors for number density and particle mass" + norms::Vector{FT} +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 end + +""" + condensation(clinfo, aps, tps, q, ρ, T) + + - `clinfo` - kwarg structure containing pdists, moments, and coalescence parameters + - `aps` - air properties + - `tps` - thermodynamics parameters + - `q` - phase partition + - `ρ` - air density + - `T` - air temperature +Returns a vector of moment tendencies due to condensation/evaporation +""" +function condensation end + +""" + weighted_vt(clinfo) + + - `clinfo` - kwarg structure containing pdists, moments, and coalescence parameters +Returns the moment-weighted terminal velocities corresponding to the rate of change of +prognostic moments +""" +function weighted_vt end + +end diff --git a/test/Project.toml b/test/Project.toml index cf9d03732..453d1d3c6 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -6,6 +6,7 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d" ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" +Cloudy = "9e3b23bb-e7cc-4b94-886c-65de2234ba87" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964" EvoTrees = "f6006082-12f8-11e9-0c9c-0d5d367ab1e5" @@ -22,4 +23,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] -KernelAbstractions = "0.9" +KernelAbstractions = "0.9" \ No newline at end of file diff --git a/test/microphysicsflexible_tests.jl b/test/microphysicsflexible_tests.jl new file mode 100644 index 000000000..7a92fb527 --- /dev/null +++ b/test/microphysicsflexible_tests.jl @@ -0,0 +1,83 @@ +import Test as TT + +import Thermodynamics as TD +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.MicrophysicsFlexible as CMF +import Cloudy.ParticleDistributions as CPD +import Cloudy.KernelFunctions as CPK + +@info "Microphysics Tests" + +function test_microphysics_flexible(FT) + # Thermodynamics and air properties parameters + aps = CMP.AirProperties(FT) + tps = TD.Parameters.ThermodynamicsParameters(FT) + + TT.@testset "Flexible microphysics - unit tests" begin + # Create Cloudy struct with different distributions + clinfo = CMF.CLSetup{FT}(; + pdists = [ + CPD.LognormalPrimitiveParticleDistribution( + FT(10.0), + FT(1.0), + FT(1.0), + ), + ], + mom = CPD.get_moments( + CPD.LognormalPrimitiveParticleDistribution( + FT(10.0), + FT(1.0), + FT(1.0), + ), + ), + NProgMoms = [3], + KernelFunc = CPK.ConstantKernelFunction(FT(1)), + mass_thresholds = [FT(100.0)], + kernel_order = 4, + kernel_limit = FT(100), + ) + + # Create Cloudy struct with defaults + clinfo = CMF.CLSetup{FT}() + TT.@test length(clinfo.pdists) == 2 + TT.@test length(clinfo.mom) == 5 + TT.@test length(clinfo.NProgMoms) == 2 + + # Test coalescence + TT.@test isnothing(clinfo.coal_data) + dmom = CMF.coalescence(clinfo) + TT.@test ~isnothing(clinfo.coal_data) + TT.@test all(dmom[1:clinfo.NProgMoms[1]] .< FT(0)) + TT.@test all(dmom[(end - 1):end] .> FT(0)) + + # Test evaporation + T = FT(288.15) + S = 0.99 + cond_evap = CMF.condensation(clinfo, aps, tps, T, S) + # conserve number density + TT.@test cond_evap[1] == FT(0) + TT.@test cond_evap[clinfo.NProgMoms[1] + 1] == FT(0) + # decrease mass density + TT.@test cond_evap[2] < FT(0) + TT.@test cond_evap[clinfo.NProgMoms[1] + 2] < FT(0) + + # Test condensation + S = 1.01 + cond_evap = CMF.condensation(clinfo, aps, tps, T, S) + # conserve number density + TT.@test cond_evap[1] == FT(0) + TT.@test cond_evap[clinfo.NProgMoms[1] + 1] == FT(0) + # increase cloud mass density + TT.@test cond_evap[2] > FT(0) + + # Test sedimentation + weighted_vt = CMF.weighted_vt(clinfo) + TT.@test all(weighted_vt .> FT(0)) + end +end + +println("Testing Float64") +test_microphysics_flexible(Float64) + +# println("Testing Float32") +# test_microphysics_flexible(Float32)