From a70f197215ff0664b7dde0cd386af7efa8d42007 Mon Sep 17 00:00:00 2001 From: Rowan Orlijan-Rhyne <129808953+rorlija1@users.noreply.github.com> Date: Sat, 24 Aug 2024 03:03:28 +0200 Subject: [PATCH] P3 skeleton framework (#185) add p3 advection case --- .buildkite/pipeline.yml | 4 +- Project.toml | 2 +- src/Common/Common.jl | 1 + src/Common/equation_types.jl | 9 +- src/Common/helper_functions.jl | 9 + src/Common/initial_condition.jl | 114 ++++++++ src/Common/netcdf_io.jl | 10 + src/Common/ode_utils.jl | 42 ++- src/Common/tendency.jl | 82 +++++- src/K1DModel/tendency.jl | 198 ++++++++++++- test/Project.toml | 2 +- .../KiD_driver/parse_commandline.jl | 16 +- .../KiD_driver/run_KiD_simulation.jl | 46 ++- test/plotting_utils.jl | 264 +++++++++++++++--- test/unit_tests/common_unit_test.jl | 84 +++++- test/unit_tests/k1d_unit_test.jl | 110 +++++++- test/unit_tests/unit_test.jl | 14 + 17 files changed, 935 insertions(+), 72 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index c8b3fca3..5786f2c7 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -135,5 +135,7 @@ steps: command: "julia --color=yes --project=test test/experiments/KiD_driver/KiD_driver.jl --moisture_choice=CloudyMoisture --precipitation_choice=CloudyPrecip --n_elem=128 --num_moments=7" artifact_paths: "test/experiments/KiD_driver/Output_CloudyMoisture_CloudyPrecip_7/figures/*" - + - label: ":snowflake: Experiment: P3 Setup" + command: "julia --color=yes --project=test test/experiments/KiD_driver/KiD_driver.jl --moisture_choice=MoistureP3 --precipitation_choice=PrecipitationP3 --n_elem=24 --z_max=3000 --t_end=75 --w1=0 --dt=0.5" + artifact_paths: "test/experiments/KiD_driver/Output_MoistureP3_PrecipitationP3/figures/*" diff --git a/Project.toml b/Project.toml index 46b5253d..f3212eba 100644 --- a/Project.toml +++ b/Project.toml @@ -30,7 +30,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] ClimaCore = "0.14" ClimaParams = "0.10.7" -CloudMicrophysics = "0.22.1" +CloudMicrophysics = "0.22.2" Cloudy = "0.5.6" DiffEqBase = "6" Distributions = "0.25" diff --git a/src/Common/Common.jl b/src/Common/Common.jl index b1023cec..772185e1 100644 --- a/src/Common/Common.jl +++ b/src/Common/Common.jl @@ -17,6 +17,7 @@ const CM1 = CM.Microphysics1M const CM2 = CM.Microphysics2M const CMAM = CM.AerosolModel const CMAA = CM.AerosolActivation +const CMP3 = CM.P3Scheme include("parameters.jl") include("equation_types.jl") diff --git a/src/Common/equation_types.jl b/src/Common/equation_types.jl index 6a045392..b865d4e7 100644 --- a/src/Common/equation_types.jl +++ b/src/Common/equation_types.jl @@ -52,9 +52,10 @@ struct Precipitation2M{PT, ST} <: AbstractPrecipitationStyle rain_formation::PT sedimentation::ST end -struct PrecipitationP3{PT, ST, P3} <: AbstractPrecipitationStyle - rain_formation::PT - sedimentation::ST - p3_parameters::P3 +struct PrecipitationP3{P3, CH, PT, B} <: AbstractPrecipitationStyle + p3_params::P3 + Chen2022::CH + sb2006::PT + p3_boundary_condition::B end struct CloudyPrecip <: AbstractPrecipitationStyle end diff --git a/src/Common/helper_functions.jl b/src/Common/helper_functions.jl index e36da5e7..cef19ead 100644 --- a/src/Common/helper_functions.jl +++ b/src/Common/helper_functions.jl @@ -8,6 +8,8 @@ function get_moisture_type(moisture_choice::String, toml_dict) moisture = NonEquilibriumMoisture(CMP.CloudLiquid(toml_dict), CMP.CloudIce(toml_dict)) elseif moisture_choice == "CloudyMoisture" moisture = CloudyMoisture() + elseif moisture_choice == "MoistureP3" + moisture = MoistureP3() else error("Invalid moisture choice: $moisture_choice") end @@ -22,6 +24,7 @@ function get_precipitation_type( toml_dict; rain_formation_choice::Union{Nothing, String} = nothing, sedimentation_choice::Union{Nothing, String} = nothing, + boundary::Union{Nothing, NamedTuple} = nothing, ) if precipitation_choice == "NoPrecipitation" precip = NoPrecipitation() @@ -77,6 +80,12 @@ function get_precipitation_type( end elseif precipitation_choice == "CloudyPrecip" precip = CloudyPrecip() + elseif precipitation_choice == "PrecipitationP3" + FT = Float64 + p3_params = CMP.ParametersP3(FT) + Chen2022 = CMP.Chen2022VelType(FT) + sb2006 = CMP.SB2006(toml_dict) + precip = PrecipitationP3(p3_params, Chen2022, sb2006, boundary) else error("Invalid precipitation choice: $precipitation_choice") end diff --git a/src/Common/initial_condition.jl b/src/Common/initial_condition.jl index 1658b82f..ccbae713 100644 --- a/src/Common/initial_condition.jl +++ b/src/Common/initial_condition.jl @@ -281,3 +281,117 @@ function cloudy_initial_condition(pdists, ip, k = 1) return merge(ip, (; moments = moments, pdists = pdists, cloudy_moments_zero)) end + +function p3_initial_condition( + ::Type{FT}, + kid_params, + thermo_params, + z; + _q_init = FT(0), + _N_init = FT(0), + _F_rim = FT(0), + _F_liq = FT(0), + _ρ_r = FT(0), + z_top, + ice_start, + dry = false, +) where {FT} + + # initialize water vapor profile + # using same setup as initial_condition_1d + q_vap::FT = init_profile(FT, kid_params, thermo_params, z, dry = dry).qv + + # if ice_start, start with small ice bump + # otherwise, introduce particles solely through + # the boundary in advection_tendency!() + _z_band = z_top * FT(0.2) + has_ice(z) = ifelse(ice_start, z > (z_top - _z_band), false) + ice_bump(χ, z) = χ * cos(z - (z_top - 0.5 * _z_band)) + + # P3 state variables (minus B_rim -- computed below) + # (initialized as kg/kg before converting to kg/m3) + q_ice::FT = has_ice(z) ? ice_bump(_q_init, z) : FT(0) + N_ice_kg::FT = has_ice(z) ? ice_bump(_N_init, z) : FT(0) + q_liqonice::FT = _F_liq * q_ice + q_rim::FT = _F_rim * (q_ice - q_liqonice) + + # 2M warm rain state variables + q_liq::FT = FT(0) + N_liq::FT = FT(0) + q_rai::FT = FT(0) + N_rai::FT = FT(0) + + # q_tot - single p3 category + 2M categories + vapor + q_tot::FT = q_ice + q_liq + q_rai + q_vap + + # thermodynamics: + # for Case 1 of Cholette et al 2019 paper + # use approximate temperature values from + # sounding: kin1d/soundings/Tsfc2.txt + # (path in p3 fortran github repo) + T::FT = -0.004 * (z - 500) + 273.15 # temperature + p::FT = 990 - (0.1 * z) # pressure + _q = TD.PhasePartition(q_tot, q_liq, q_ice) + ρ::FT = TD.air_density(thermo_params, T, p, _q) # total density + ρ_dry::FT = TD.air_density(thermo_params, T, p) + ts = TD.PhaseNonEquil_ρTq(thermo_params, ρ, T, _q) # thermodynamic state + + # P3 scheme uses L (kg/m3) = ρ * q + # so we compute ρq and pass that to P3 + # also we convert N to (1/m3) + # since paper init values are given + # in (kg/kg), (1/kg) + ρq_tot::FT = ρ * q_tot + N_ice::FT = ρ * N_ice_kg + ρq_ice::FT = ρ * q_ice + ρq_rim::FT = ρ * q_rim + ρq_liqonice::FT = ρ * q_liqonice + ρq_rai::FT = ρ * q_rai + ρq_liq::FT = ρ * q_liq + ρq_vap::FT = ρ * q_vap + + # also compute B_rim from L_rim, ρ_r + B_rim::FT = ρq_rim / _ρ_r + + + # unused quantities: + q_sno::FT = FT(0) + ρq_sno::FT = FT(0) + θ_liq_ice::FT = TD.liquid_ice_pottemp(thermo_params, ts) + N_aer::FT = FT(0) + θ_dry::FT = TD.dry_pottemp(thermo_params, T, ρ) + + # zero: + zero::FT = FT(0) + + return (; + ρ, + ρ_dry, + T, + p, + θ_liq_ice, + θ_dry, + ρq_tot, + ρq_liq, + ρq_ice, + ρq_rai, + ρq_rim, + ρq_sno, + ρq_vap, + ρq_liqonice, + q_tot, + q_liq, + q_ice, + q_rai, + q_sno, + q_rim, + q_liqonice, + q_vap, + B_rim, + N_liq, + N_ice, + N_rai, + N_aer, + zero, + ) +end diff --git a/src/Common/netcdf_io.jl b/src/Common/netcdf_io.jl index 7ebe13db..272ebad6 100644 --- a/src/Common/netcdf_io.jl +++ b/src/Common/netcdf_io.jl @@ -23,11 +23,21 @@ function NetCDFIO_Stats( :q_tot => "q_tot", :q_liq => "q_liq", :q_ice => "q_ice", + :q_rim => "q_rim", + :q_liqonice => "q_liqonice", :q_rai => "q_rai", :q_sno => "q_sno", + :q_vap => "q_vap", :N_aer => "N_aer", :N_liq => "N_liq", :N_rai => "N_rai", + :N_ice => "N_ice", + :B_rim => "B_rim", + :ρq_ice => "ρq_ice", + :ρq_rim => "ρq_rim", + :ρq_liqonice => "ρq_liqonice", + :ρq_vap => "ρq_vap", + :ρq_rai => "ρq_rai", ), ) FT = Float64 diff --git a/src/Common/ode_utils.jl b/src/Common/ode_utils.jl index ab1fdcfd..0f4b4b0a 100644 --- a/src/Common/ode_utils.jl +++ b/src/Common/ode_utils.jl @@ -76,11 +76,15 @@ function initialise_state(::MoistureP3, ::PrecipitationP3, initial_profiles) ρq_rai = initial_profiles.ρq_rai, ρq_ice = initial_profiles.ρq_ice, ρq_rim = initial_profiles.ρq_rim, + ρq_liqonice = initial_profiles.ρq_liqonice, B_rim = initial_profiles.B_rim, N_liq = initial_profiles.N_liq, N_rai = initial_profiles.N_rai, N_ice = initial_profiles.N_ice, N_aer = initial_profiles.N_aer, + q_vap = initial_profiles.q_vap, + ρq_vap = initial_profiles.ρq_vap, + q_rai = initial_profiles.q_rai, ) end function initialise_state(::CloudyMoisture, ::CloudyPrecip, initial_profiles) @@ -116,7 +120,7 @@ function initialise_aux( if moisture isa EquilibriumMoisture ts = @. TD.PhaseEquil_ρθq(thermo_params, ip.ρ, ip.θ_liq_ice, ip.q_tot) cloud_sources = nothing - elseif moisture isa NonEquilibriumMoisture + elseif moisture isa NonEquilibriumMoisture || moisture isa MoistureP3 q = @. TD.PhasePartition(ip.q_tot, ip.q_liq, ip.q_ice) ts = @. TD.PhaseNonEquil_ρθq(thermo_params, ip.ρ, ip.θ_liq_ice, q) @@ -181,11 +185,20 @@ function initialise_aux( q_rai = ip.q_rai, q_ice = ip.q_ice, q_rim = ip.q_rim, + q_liqonice = ip.q_liqonice, B_rim = ip.B_rim, N_liq = ip.N_liq, N_rai = ip.N_rai, N_ice = ip.N_ice, N_aer = ip.N_aer, + ρq_tot = ip.ρq_tot, + ρq_liq = ip.ρq_liq, + ρq_rai = ip.ρq_rai, + ρq_ice = ip.ρq_ice, + ρq_rim = ip.ρq_rim, + ρq_liqonice = ip.ρq_liqonice, + q_vap = ip.q_vap, + ρq_vap = ip.ρq_vap, ) velocities = (; term_vel_rai = copy(ip.zero), @@ -199,6 +212,7 @@ function initialise_aux( q_rai::FT, q_ice::FT, q_rim::FT, + q_liqonice::FT, N_aer::FT, N_liq::FT, N_rai::FT, @@ -217,11 +231,31 @@ function initialise_aux( copy(ip.zero), copy(ip.zero), copy(ip.zero), + copy(ip.zero), + ), + ) + p3_boundary_condition_eltype = @NamedTuple{ + ice_start::Bool, + _magnitude::FT, + _q_flux::FT, + _N_flux::FT, + _F_rim::FT, + _F_liq::FT, + _ρ_r_init::FT, + } + p3_boundary_condition = @. p3_boundary_condition_eltype( + tuple( + copy(precip.p3_boundary_condition.ice_start), + copy(precip.p3_boundary_condition._magnitude), + copy(precip.p3_boundary_condition._q_flux), + copy(precip.p3_boundary_condition._N_flux), + copy(precip.p3_boundary_condition._F_rim), + copy(precip.p3_boundary_condition._F_liq), + copy(precip.p3_boundary_condition._ρ_r_init), ), ) activation_sources = nothing - @assert moisture isa NonEquilibrumMoisture elseif precip isa CloudyPrecip microph_variables = (; q_tot = ip.q_tot, @@ -241,7 +275,7 @@ function initialise_aux( cloudy_variables = (; nm_cloud = Val(cloudy_params.NProgMoms[1])) scratch = merge(scratch, (; tmp_cloudy = similar(ip.cloudy_moments_zero))) else - error("Wrong precipitation choise $precip") + error("Wrong precipitation choice $precip") end aux = (; @@ -262,6 +296,8 @@ function initialise_aux( if precip isa CloudyPrecip aux = merge(aux, (; cloudy_params, cloudy_variables)) + elseif precip isa PrecipitationP3 + aux = merge(aux, (; p3_boundary_condition)) end return aux diff --git a/src/Common/tendency.jl b/src/Common/tendency.jl index 0487fd99..9807ab81 100644 --- a/src/Common/tendency.jl +++ b/src/Common/tendency.jl @@ -47,8 +47,41 @@ end @. θ_dry = TD.dry_pottemp(thermo_params, T, ρ_dry) @. θ_liq_ice = TD.liquid_ice_pottemp(thermo_params, ts) end +@inline function precompute_aux_thermo!(::MoistureP3, Y, aux) + + FT = eltype(Y.ρq_liq) + + (; thermo_params) = aux + (; ts, ρ, ρ_dry, p, T, θ_dry, θ_liq_ice) = aux.thermo_variables + (; ρq_liq, ρq_ice, ρq_rim, ρq_liqonice, ρq_vap, ρq_rai, q_tot, q_liq, q_ice, q_rim, q_liqonice, q_vap, q_rai) = + aux.microph_variables + + @. ρ = ρ_dry + Y.ρq_tot + @. ρq_rim = Y.ρq_rim + @. ρq_ice = Y.ρq_ice + @. ρq_liqonice = Y.ρq_liqonice + @. ρq_liq = Y.ρq_liq + @. ρq_rai = Y.ρq_rai + @. ρq_vap = Y.ρq_vap + + @. q_liq = q_(Y.ρq_liq, ρ) + @. q_ice = q_(Y.ρq_ice, ρ) + @. q_rim = q_(Y.ρq_rim, ρ) + @. q_liqonice = q_(Y.ρq_liqonice, ρ) + @. q_vap = q_(Y.ρq_vap, ρ) + @. q_rai = q_(Y.ρq_rai, ρ) + @. q_tot = q_liq + q_ice + q_rai + q_vap + + @. ts = TD.PhaseNonEquil_ρTq(thermo_params, ρ, T, PP(q_tot, q_liq, q_ice)) + @. p = TD.air_pressure(thermo_params, ts) + @. θ_dry = TD.dry_pottemp(thermo_params, T, ρ_dry) + @. θ_liq_ice = TD.liquid_ice_pottemp(thermo_params, ts) + +end @inline function precompute_aux_thermo!(::NonEquilibriumMoisture, Y, aux) + FT = eltype(Y.ρq_liq) + (; thermo_params) = aux (; ts, ρ, ρ_dry, p, T, θ_dry, θ_liq_ice) = aux.thermo_variables (; q_tot, q_liq, q_ice) = aux.microph_variables @@ -143,19 +176,53 @@ end end @inline function precompute_aux_precip!(ps::PrecipitationP3, Y, aux) + # update state variables FT = eltype(Y.ρq_rai) - (; ρ) = aux.thermo_variables - (; q_rai, q_ice, q_rim, B_rim, N_rai, N_ice, N_liq) = aux.microph_variables + (; ρ, ρ_dry) = aux.thermo_variables + (; ρq_ice, ρq_rim, N_ice, ρq_liqonice, B_rim) = aux.microph_variables + @. ρq_ice = Y.ρq_ice + @. ρq_rim = Y.ρq_rim + @. ρq_liqonice = Y.ρq_liqonice + @. N_ice = max(FT(0), Y.N_ice) + @. B_rim = max(FT(0), Y.B_rim) + + # compute predicted quantities + ρ_r = aux.scratch.tmp + F_rim = aux.scratch.tmp2 + F_liq = aux.scratch.tmp3 + # prohibit ρ_r < 0 + @. ρ_r = ifelse(B_rim < eps(FT), eps(FT), max(FT(0), ρq_rim / B_rim)) + # keep F_r in range [0, 1) + @. F_rim = + ifelse((ρq_ice - ρq_liqonice) < eps(FT), FT(0), min(max(FT(0), ρq_rim / (ρq_ice - ρq_liqonice)), 1 - eps(FT))) + # prohibit F_liq < 0 + @. F_liq = ifelse(Y.ρq_ice < eps(FT), FT(0), Y.ρq_liqonice / Y.ρq_ice) + + p3 = ps.p3_params + Chen2022 = ps.Chen2022 + + (; term_vel_ice, term_vel_N_ice, term_vel_rai, term_vel_N_rai) = aux.velocities + use_aspect_ratio = true + @. term_vel_N_ice = + getindex(CMP3.ice_terminal_velocity(p3, Chen2022, ρq_ice, N_ice, ρ_r, F_rim, F_liq, ρ_dry, use_aspect_ratio), 1) + @. term_vel_ice = + getindex(CMP3.ice_terminal_velocity(p3, Chen2022, ρq_ice, N_ice, ρ_r, F_rim, F_liq, ρ_dry, use_aspect_ratio), 2) + + # adding 2M rain below: + + sb2006 = ps.sb2006 + + (; q_rai, q_liq, N_rai, N_liq) = aux.microph_variables + (; term_vel_rai, term_vel_N_rai) = aux.velocities @. q_rai = q_(Y.ρq_rai, ρ) - @. q_ice = q_(Y.ρq_ice, ρ) - @. q_rim = q_(Y.ρq_rim, ρ) - @. B_rim = q_(Y.ρq_rim, ρ) + @. q_liq = q_(Y.ρq_liq, ρ) @. N_rai = max(FT(0), Y.N_rai) @. N_liq = max(FT(0), Y.N_liq) - @. N_ice = max(FT(0), Y.N_ice) - # TODO... + @. term_vel_N_rai = getindex(CM2.rain_terminal_velocity(sb2006, Chen2022.rain, q_rai, ρ, N_rai), 1) + @. term_vel_rai = getindex(CM2.rain_terminal_velocity(sb2006, Chen2022.rain, q_rai, ρ, N_rai), 2) + end function get_dists_moments(moments, NProgMoms::NTuple{ND, Int}) where {ND} FT = eltype(moments) @@ -492,6 +559,7 @@ end end @inline function precompute_aux_precip_sources!(ps::PrecipitationP3, aux) + # TODO [P3] return nothing end diff --git a/src/K1DModel/tendency.jl b/src/K1DModel/tendency.jl index 7d86a47b..480ac5c1 100644 --- a/src/K1DModel/tendency.jl +++ b/src/K1DModel/tendency.jl @@ -95,7 +95,7 @@ Aerosol activation tendencies error("activation_tendency not implemented for a given $sp") end @inline function precompute_aux_activation!( - ::Union{CO.NoPrecipitation, CO.Precipitation0M, CO.Precipitation1M}, + ::Union{CO.NoPrecipitation, CO.Precipitation0M, CO.Precipitation1M, CO.PrecipitationP3}, dY, Y, aux, @@ -269,6 +269,25 @@ end return dY end +@inline function advection_tendency!(::CO.MoistureP3, dY, Y, aux, t) + FT = eltype(Y.ρq_tot) + + # advect aerosols + + If = CC.Operators.InterpolateC2F() + ∂ = CC.Operators.DivergenceF2C( + bottom = CC.Operators.SetValue(CC.Geometry.WVector(aux.prescribed_velocity.ρw0 * aux.q_surf)), + top = CC.Operators.Extrapolate(), + ) + @. dY.N_aer += -∂(aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ) * If(Y.N_aer)) + @. dY.ρq_vap += -∂(aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ) * If(Y.ρq_vap)) + + fcc = CC.Operators.FluxCorrectionC2C(bottom = CC.Operators.Extrapolate(), top = CC.Operators.Extrapolate()) + @. dY.N_aer += fcc(aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ), Y.N_aer) + @. dY.ρq_vap += fcc(aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ), Y.ρq_vap) + + return dY +end @inline function advection_tendency!(::Union{CO.NoPrecipitation, CO.Precipitation0M}, dY, Y, aux, t) end @inline function advection_tendency!(::CO.Precipitation1M, dY, Y, aux, t) FT = eltype(Y.ρq_tot) @@ -393,3 +412,180 @@ end return dY end + +@inline function advection_tendency!(::CO.PrecipitationP3, dY, Y, aux, t) + FT = eltype(Y.ρq_tot) + # TODO - flux magnitude at top should agree + # with ρ * q instead of being only based on q + + # P3 advection (introducing particles through boundary): + # TODO - change run_KiD_simulation so that temp. variables + # (_q_flux, _N_flux, etc) are not hard coded in the advection tendency + + # update aux with state variables + (; ρ) = aux.thermo_variables + (; ρq_tot, ρq_liq, ρq_rai, ρq_ice, ρq_rim, ρq_liqonice, N_ice, N_rai, N_aer, N_liq, B_rim, ρq_vap, q_rai) = + aux.microph_variables + @. ρq_ice = Y.ρq_ice + @. ρq_rim = Y.ρq_rim + @. ρq_liqonice = Y.ρq_liqonice + @. N_ice = Y.N_ice + @. B_rim = Y.B_rim + + # choose flux characteristics + (; ice_start, _magnitude, _q_flux, _N_flux, _F_rim, _F_liq, _ρ_r_init) = aux.p3_boundary_condition + # if we have initial ice signal then do not introduce boundary flux + if ice_start + _magnitude = FT(0) + end + + If = CC.Operators.InterpolateC2F() + + # define divergence operators with different boundary conditions + # for each state variable + ∂q_ice = CC.Operators.DivergenceF2C( + bottom = CC.Operators.Extrapolate(), + top = CC.Operators.SetValue(CC.Geometry.WVector(_magnitude * _q_flux)), + ) + ∂q_rim = CC.Operators.DivergenceF2C( + bottom = CC.Operators.Extrapolate(), + top = CC.Operators.SetValue(CC.Geometry.WVector(_magnitude * _F_rim * (1 - _F_liq) * _q_flux)), + ) + ∂q_liqonice = CC.Operators.DivergenceF2C( + bottom = CC.Operators.Extrapolate(), + top = CC.Operators.SetValue(CC.Geometry.WVector(_magnitude * _F_liq * _q_flux)), + ) + ∂B_rim = CC.Operators.DivergenceF2C( + bottom = CC.Operators.Extrapolate(), + top = CC.Operators.SetValue(CC.Geometry.WVector(_magnitude * _F_rim * (1 - _F_liq) * _q_flux / _ρ_r_init)), + ) + ∂N_ice = CC.Operators.DivergenceF2C( + bottom = CC.Operators.Extrapolate(), + top = CC.Operators.SetValue(CC.Geometry.WVector(_magnitude * _N_flux)), + ) + + # apply divergence + @. dY.ρq_ice += ∂q_ice( + FT(-1) * (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_ice))) * If(Y.ρq_ice), + ) + @. dY.ρq_rim += ∂q_rim( + FT(-1) * (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_ice))) * If(Y.ρq_rim), + ) + @. dY.B_rim += ∂B_rim( + FT(-1) * (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_ice))) * If(Y.B_rim), + ) + @. dY.ρq_liqonice += ∂q_liqonice( + FT(-1) * (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_ice))) * If(Y.ρq_liqonice), + ) + @. dY.N_ice += ∂N_ice( + FT(-1) * (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_N_ice))) * If(Y.N_ice), + ) + fcc = CC.Operators.FluxCorrectionC2C(bottom = CC.Operators.Extrapolate(), top = CC.Operators.Extrapolate()) + + # apply flux correction + @. dY.ρq_ice += fcc( + (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_ice) * FT(-1))), + Y.ρq_ice, + ) + @. dY.ρq_rim += fcc( + (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_ice) * FT(-1))), + Y.ρq_rim, + ) + @. dY.ρq_liqonice += fcc( + (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_ice) * FT(-1))), + Y.ρq_liqonice, + ) + @. dY.B_rim += fcc( + (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_ice) * FT(-1))), + Y.B_rim, + ) + @. dY.N_ice += fcc( + (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_N_ice) * FT(-1))), + Y.N_ice, + ) + + # 2M rain advection (zero boundary flux): + + @. ρq_rai = Y.ρq_rai + @. ρq_liq = Y.ρq_liq + @. N_rai = Y.N_rai + @. N_liq = Y.N_liq + @. N_aer = Y.N_aer + + ∂ = CC.Operators.DivergenceF2C( + bottom = CC.Operators.Extrapolate(), + top = CC.Operators.SetValue(CC.Geometry.WVector(0.0)), + ) + ∂ₐ = CC.Operators.DivergenceF2C(bottom = CC.Operators.Extrapolate(), top = CC.Operators.Extrapolate()) + + @. dY.N_aer += -∂ₐ((aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) * If(Y.N_aer)) + @. dY.N_liq += -∂((aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) * If(Y.N_liq)) + + @. dY.N_rai += + -∂( + ( + aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ) + + CC.Geometry.WVector(If(aux.velocities.term_vel_N_rai) * FT(-1)) + ) * If(Y.N_rai), + ) + @. dY.q_rai += + -∂( + ( + aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ) + + CC.Geometry.WVector(If(aux.velocities.term_vel_rai) * FT(-1)) + ) * If(Y.ρq_rai), + ) + + @. dY.N_aer += fcc((aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)), Y.N_aer) + @. dY.N_liq += fcc((aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)), Y.N_liq) + + @. dY.N_rai += fcc( + ( + aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ) + + CC.Geometry.WVector(If(aux.velocities.term_vel_N_rai) * FT(-1)) + ), + Y.N_rai, + ) + @. dY.ρq_rai += fcc( + ( + aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ) + + CC.Geometry.WVector(If(aux.velocities.term_vel_rai) * FT(-1)) + ), + Y.ρq_rai, + ) + + # advecting q_tot: + + @. ρq_vap = Y.ρq_vap + @. ρq_tot = ρq_ice + ρq_rai + ρq_vap + ρq_liq + + ∂q_tot = CC.Operators.DivergenceF2C( + bottom = CC.Operators.Extrapolate(), + top = CC.Operators.SetValue(CC.Geometry.WVector(_magnitude * _q_flux)), + ) + + @. dY.ρq_tot += ∂q_tot( + FT(-1) * (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_ice))) * If(Y.ρq_ice) + + CC.Geometry.WVector(If(aux.velocities.term_vel_rai)) * If(ρq_rai), + ) + + @. dY.ρq_tot += fcc( + (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) + + (CC.Geometry.WVector(If(aux.velocities.term_vel_ice) * FT(-1))) + + CC.Geometry.WVector(If(aux.velocities.term_vel_rai) * FT(-1)), + Y.ρq_tot, + ) + + return dY +end diff --git a/test/Project.toml b/test/Project.toml index f2c65081..f966ebcd 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -31,7 +31,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" ClimaCore = "0.14" ClimaCorePlots = "0.2" ClimaParams = "0.10.7" -CloudMicrophysics = "0.22.1" +CloudMicrophysics = "0.22.2" Cloudy = "0.5.6" DiffEqBase = "6.75" NCDatasets = "0.14" diff --git a/test/experiments/KiD_driver/parse_commandline.jl b/test/experiments/KiD_driver/parse_commandline.jl index 4027d7c5..0559b997 100644 --- a/test/experiments/KiD_driver/parse_commandline.jl +++ b/test/experiments/KiD_driver/parse_commandline.jl @@ -13,11 +13,11 @@ function parse_commandline() arg_type = String default = "Float64" "--moisture_choice" - help = "Mositure model choice: EquilibriumMoisture, NonEquilibriumMoisture, CloudyMoisture" + help = "Mositure model choice: EquilibriumMoisture, NonEquilibriumMoisture, CloudyMoisture, MoistureP3" arg_type = String default = "EquilibriumMoisture" "--precipitation_choice" - help = "Precipitation model choice: NoPrecipitation, Precipitation0M, Precipitation1M, Precipitation2M, CloudyPrecip" + help = "Precipitation model choice: NoPrecipitation, Precipitation0M, Precipitation1M, Precipitation2M, CloudyPrecip, PrecipitationP3" arg_type = String default = "Precipitation1M" "--num_moments" @@ -145,6 +145,18 @@ function parse_commandline() help = "Initial condition theta2 [K]" arg_type = Float64 default = Float64(312.66) + "--p3_boundary_condition" + help = "Characteristics of particle flux being introduced into P3 domain top (if ice_start = false) or of the initial ice signal (if ice_start = true)" + arg_type = NamedTuple + default = (; + ice_start = false, + _magnitude = Float64(0.5), + _q_flux = Float64(0.65e-4), + _N_flux = Float64(40000), + _F_rim = Float64(0.2), + _F_liq = Float64(0.2), + _ρ_r_init = Float64(900), + ) end return AP.parse_args(s) diff --git a/test/experiments/KiD_driver/run_KiD_simulation.jl b/test/experiments/KiD_driver/run_KiD_simulation.jl index 6e59c994..c897844b 100644 --- a/test/experiments/KiD_driver/run_KiD_simulation.jl +++ b/test/experiments/KiD_driver/run_KiD_simulation.jl @@ -36,6 +36,11 @@ function run_KiD_simulation(::Type{FT}, opts) where {FT} if opts["open_system_activation"] output_folder = output_folder * "_OSA" end + if precipitation_choice == "PrecipitationP3" + p3_boundary_condition = opts["p3_boundary_condition"] + else + p3_boundary_condition = nothing + end path = joinpath(@__DIR__, output_folder) mkpath(path) @@ -79,6 +84,7 @@ function run_KiD_simulation(::Type{FT}, opts) where {FT} toml_dict; rain_formation_choice = rain_formation_choice, sedimentation_choice = sedimentation_choice, + boundary = p3_boundary_condition, ) @info "Initialising" @@ -107,6 +113,25 @@ function run_KiD_simulation(::Type{FT}, opts) where {FT} ), coord, ) + elseif precipitation_choice == "PrecipitationP3" + cloudy_params = nothing + (; ice_start, _q_flux, _N_flux, _F_rim, _F_liq, _ρ_r_init) = precip.p3_boundary_condition + init = map( + coord -> CO.p3_initial_condition( + FT, + kid_params, + thermo_params, + coord.z; + _q_init = _q_flux, + _N_init = _N_flux, + _F_rim = _F_rim, + _F_liq = _F_liq, + _ρ_r = _ρ_r_init, + z_top = FT(opts["z_max"]), + ice_start = ice_start, + ), + coord, + ) else cloudy_params = nothing init = map( @@ -166,12 +191,21 @@ function run_KiD_simulation(::Type{FT}, opts) where {FT} z_centers = parent(CC.Fields.coordinate_field(space)) plot_final_aux_profiles(z_centers, aux, precip, output = plot_folder) - plot_animation(string("experiments/KiD_driver/", output_folder, "/Output.nc"), output = plot_folder) - plot_timeheight( - string("experiments/KiD_driver/", output_folder, "/Output.nc"), - output = plot_folder, - mixed_phase = false, - ) + if precip isa CO.PrecipitationP3 + plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, plot_folder) + plot_timeheight_p3( + string("experiments/KiD_driver/", output_folder, "/Output.nc"), + precip, + output = plot_folder, + ) + else + plot_animation(string("experiments/KiD_driver/", output_folder, "/Output.nc"), output = plot_folder) + plot_timeheight( + string("experiments/KiD_driver/", output_folder, "/Output.nc"), + output = plot_folder, + mixed_phase = false, + ) + end end return diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index c3b52abd..9c078ae2 100644 --- a/test/plotting_utils.jl +++ b/test/plotting_utils.jl @@ -66,6 +66,14 @@ function plot_final_aux_profiles(z_centers, aux, precip; output = "output") elseif precip isa CO.Precipitation2M q_rai_end = parent(aux.microph_variables.q_rai) q_sno_end = q_tot_end .* 0.0 + elseif precip isa CO.PrecipitationP3 + q_rai_end = parent(aux.microph_variables.q_rai) + q_liqonice_end = parent(aux.microph_variables.q_liqonice) + q_rim_end = parent(aux.microph_variables.q_rim) + B_rim_end = parent(aux.microph_variables.B_rim) + N_ice_end = parent(aux.microph_variables.N_ice) + N_liq_end = parent(aux.microph_variables.N_liq) + N_rai_end = parent(aux.microph_variables.N_rai) else q_rai_end = q_tot_end .* 0.0 q_sno_end = q_tot_end .* 0.0 @@ -86,54 +94,168 @@ function plot_final_aux_profiles(z_centers, aux, precip; output = "output") p3 = Plots.plot(q_ice_end .* 1e3, z_centers, xlabel = "q_ice [g/kg]", ylabel = "z [m]") p4 = Plots.plot(T_end, z_centers, xlabel = "T [K]", ylabel = "z [m]") p5 = Plots.plot(q_rai_end .* 1e3, z_centers, xlabel = "q_rai [g/kg]", ylabel = "z [m]") - p6 = Plots.plot(q_sno_end .* 1e3, z_centers, xlabel = "q_sno [g/kg]", ylabel = "z [m]") - p8 = Plots.plot(N_aer_end .* 1e-6, z_centers, xlabel = "N_aer [1/cm3]", ylabel = "z [m]") - p9 = Plots.plot(N_liq_end .* 1e-6, z_centers, xlabel = "N_liq [1/cm3]", ylabel = "z [m]") - p10 = Plots.plot(N_rai_end .* 1e-6, z_centers, xlabel = "N_rai [1/cm3]", ylabel = "z [m]") + if precip isa CO.PrecipitationP3 + p6 = Plots.plot(N_ice_end .* 1e-6, z_centers, xlabel = "N_ice [1/cm3]", ylabel = "z [m]") + p11 = Plots.plot(N_ice_end .* 1e-6, z_centers, xlabel = "N_ice [1/cm3]", ylabel = "z [m]") + p12 = Plots.plot(q_liqonice_end .* 1e3, z_centers, xlabel = "q_liqonice [g/kg]", ylabel = "z [m]") + p13 = Plots.plot(q_rim_end .* 1e3, z_centers, xlabel = "q_rim [g/kg]", ylabel = "z [m]") + p14 = Plots.plot(B_rim_end, z_centers, xlabel = "B_rim [-]", ylabel = "z [m]") + p = Plots.plot( + p1, + p2, + p3, + p4, + p5, + p6, + p11, + p12, + p13, + p14, + size = (1800.0, 1200.0), + bottom_margin = 40.0 * Plots.PlotMeasures.px, + left_margin = 80.0 * Plots.PlotMeasures.px, + layout = (2, 5), + ) + Plots.png(p, joinpath(path, "final_aux_profiles.png")) + else + p6 = Plots.plot(q_sno_end .* 1e3, z_centers, xlabel = "q_sno [g/kg]", ylabel = "z [m]") + p8 = Plots.plot(N_aer_end .* 1e-6, z_centers, xlabel = "N_aer [1/cm3]", ylabel = "z [m]") + p9 = Plots.plot(N_liq_end .* 1e-6, z_centers, xlabel = "N_liq [1/cm3]", ylabel = "z [m]") + p10 = Plots.plot(N_rai_end .* 1e-6, z_centers, xlabel = "N_rai [1/cm3]", ylabel = "z [m]") + + p7 = Plots.plot(xlabel = "precipitation susceptibility", ylabel = "z [m]") + if precip isa CO.Precipitation2M + N_liq_end = parent(aux.microph_variables.N_liq) + precip_sus_aut = + CMPS.precipitation_susceptibility_autoconversion.( + Ref(precip.rain_formation), + q_liq_end, + q_rai_end, + ρ_end, + N_liq_end, + ) + precip_sus_acc = + CMPS.precipitation_susceptibility_accretion.( + Ref(precip.rain_formation), + q_liq_end, + q_rai_end, + ρ_end, + N_liq_end, + ) + Plots.plot!([r.d_ln_pp_d_ln_q_liq for r in precip_sus_aut], z_centers, label = "aut, q_liq", color = :red) + Plots.plot!([r.d_ln_pp_d_ln_q_rai for r in precip_sus_aut], z_centers, label = "aut, q_rai", color = :brown) + Plots.plot!([r.d_ln_pp_d_ln_q_liq for r in precip_sus_acc], z_centers, label = "acc, q_liq", color = :blue) + Plots.plot!([r.d_ln_pp_d_ln_q_rai for r in precip_sus_acc], z_centers, label = "acc, q_rai", color = :green) + Plots.plot!(legend = :outerright) + end + p = Plots.plot( + p1, + p2, + p3, + p4, + p5, + p6, + p7, + p8, + p9, + p10, + size = (1800.0, 1200.0), + bottom_margin = 40.0 * Plots.PlotMeasures.px, + left_margin = 80.0 * Plots.PlotMeasures.px, + ) + Plots.png(p, joinpath(path, "final_aux_profiles.png")) + end +end + +function plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, output = plot_folder) + + path = joinpath(@__DIR__, output) + mkpath(path) + + ρ = parent(aux.thermo_variables.ρ) + nz = length(z_centers) + nt = length(solver.u) + q_tot = zeros(nz, nt) + q_liq = zeros(nz, nt) + q_ice = zeros(nz, nt) + q_rai = zeros(nz, nt) + q_sno = zeros(nz, nt) + N_rai = zeros(nz, nt) + N_liq = zeros(nz, nt) + N_ice = zeros(nz, nt) + ρq_tot = zeros(nz, nt) + ρq_liq = zeros(nz, nt) + ρq_ice = zeros(nz, nt) + ρq_liqonice = zeros(nz, nt) + ρq_rim = zeros(nz, nt) + ρq_rai = zeros(nz, nt) + B_rim = zeros(nz, nt) + + for (i, u) in enumerate(solver.u) + + + ρq_tot[:, i] = parent(u.ρq_tot) .* 1e3 + ρq_liq[:, i] = parent(u.ρq_liq) .* 1e3 + ρq_ice[:, i] = parent(u.ρq_ice) .* 1e3 + ρq_liqonice[:, i] = parent(u.ρq_liqonice) .* 1e3 + ρq_rim[:, i] = parent(u.ρq_rim) .* 1e3 + ρq_rai[:, i] = parent(u.ρq_rai) .* 1e3 + B_rim[:, i] = parent(u.B_rim) .* 1e3 + N_rai[:, i] = parent(u.N_rai) ./ 1e6 + N_liq[:, i] = parent(u.N_liq) ./ 1e6 + N_ice[:, i] = parent(u.N_ice) ./ 1e6 - p7 = Plots.plot(xlabel = "precipitation susceptibility", ylabel = "z [m]") - if precip isa CO.Precipitation2M - N_liq_end = parent(aux.microph_variables.N_liq) - precip_sus_aut = - CMPS.precipitation_susceptibility_autoconversion.( - Ref(precip.rain_formation), - q_liq_end, - q_rai_end, - ρ_end, - N_liq_end, - ) - precip_sus_acc = - CMPS.precipitation_susceptibility_accretion.( - Ref(precip.rain_formation), - q_liq_end, - q_rai_end, - ρ_end, - N_liq_end, - ) - Plots.plot!([r.d_ln_pp_d_ln_q_liq for r in precip_sus_aut], z_centers, label = "aut, q_liq", color = :red) - Plots.plot!([r.d_ln_pp_d_ln_q_rai for r in precip_sus_aut], z_centers, label = "aut, q_rai", color = :brown) - Plots.plot!([r.d_ln_pp_d_ln_q_liq for r in precip_sus_acc], z_centers, label = "acc, q_liq", color = :blue) - Plots.plot!([r.d_ln_pp_d_ln_q_rai for r in precip_sus_acc], z_centers, label = "acc, q_rai", color = :green) - Plots.plot!(legend = :outerright) end - p = Plots.plot( - p1, - p2, - p3, - p4, - p5, - p6, - p8, - p9, - p10, - p7, - size = (1800.0, 1200.0), - bottom_margin = 40.0 * Plots.PlotMeasures.px, - left_margin = 80.0 * Plots.PlotMeasures.px, - ) - Plots.png(p, joinpath(path, "final_aux_profiles.png")) + function plot_data(data, data_label, max_val, title = "") + return Plots.plot( + data, + z_centers, + xlabel = data_label, + ylabel = "z [m]", + legend = false, + title = title, + titlefontsize = 30, + xlim = [0, 1.1 * max_val], + ) + end + + anim = Plots.@animate for i in 1:nt + + title = "time = " * string(floor(Int, solver.t[i])) * " [s]" + + p1 = plot_data(ρq_tot[:, i], "ρq_tot [g/m3]", maximum(ρq_tot)) + p2 = plot_data(ρq_liq[:, i], "ρq_liq [g/m3]", maximum(ρq_liq)) + p3 = plot_data(ρq_ice[:, i], "ρq_ice [g/m3]", maximum(ρq_ice), title) + p4 = plot_data(ρq_liqonice[:, i], "ρq_liqonice [g/m3]", maximum(ρq_liqonice)) + p5 = plot_data(ρq_rim[:, i], "ρq_rim [g/m3]", maximum(ρq_rim)) + p6 = plot_data(ρq_rai[:, i], "ρq_rai [g/m3]", maximum(ρq_rai)) + p7 = plot_data(B_rim[:, i], "B_rim [-]", maximum(B_rim)) + p8 = plot_data(N_liq[:, i], "N_liq [1/cm^3]", maximum(N_liq)) + p9 = plot_data(N_ice[:, i], "N_ice [1/cm^3]", maximum(N_ice)) + p10 = plot_data(N_rai[:, i], "N_rai [1/cm^3]", maximum(N_rai)) + plot( + p1, + p2, + p3, + p4, + p5, + p6, + p7, + p8, + p9, + p10, + size = (1800.0, 1500.0), + bottom_margin = 30.0 * Plots.PlotMeasures.px, + left_margin = 30.0 * Plots.PlotMeasures.px, + top_margin = 30.0 * Plots.PlotMeasures.px, + right_margin = 30.0 * Plots.PlotMeasures.px, + layout = (5, 2), + ) + end + + Plots.mp4(anim, joinpath(path, "animation.mp4"), fps = 10) end function plot_animation(nc_data_file; output = "output") @@ -198,6 +320,60 @@ function plot_animation(nc_data_file; output = "output") Plots.mp4(anim, joinpath(path, "animation.mp4"), fps = 10) end +function plot_timeheight_p3(nc_data_file, precip; output = "output") + path = joinpath(@__DIR__, output) + mkpath(path) + ds = NC.NCDataset(joinpath(@__DIR__, nc_data_file)) + t_plt = Array(ds.group["profiles"]["t"]) + z_plt = Array(ds.group["profiles"]["zc"]) + q_tot_plt = Array(ds.group["profiles"]["q_tot"]) + q_liq_plt = Array(ds.group["profiles"]["q_liq"]) + ρq_ice_plt = Array(ds.group["profiles"]["ρq_ice"]) + ρq_rim_plt = Array(ds.group["profiles"]["ρq_rim"]) + ρq_liqonice_plt = Array(ds.group["profiles"]["ρq_liqonice"]) + q_rai_plt = Array(ds.group["profiles"]["q_rai"]) + q_vap_plt = Array(ds.group["profiles"]["q_vap"]) + N_aer_plt = Array(ds.group["profiles"]["N_aer"]) + N_liq_plt = Array(ds.group["profiles"]["N_liq"]) + N_rai_plt = Array(ds.group["profiles"]["N_rai"]) + N_ice_plt = Array(ds.group["profiles"]["N_ice"]) + B_rim_plt = Array(ds.group["profiles"]["B_rim"]) + #! format: off + p1 = Plots.heatmap(t_plt, z_plt, q_tot_plt .* 1e3, title = "q_tot [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + p2 = Plots.heatmap(t_plt, z_plt, q_liq_plt .* 1e3, title = "q_liq [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + p3 = Plots.heatmap(t_plt, z_plt, ρq_ice_plt .* 1e3, title = "ρq_ice [g/m3]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + p4 = Plots.heatmap(t_plt, z_plt, q_rai_plt .* 1e3, title = "q_rai [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + p5 = Plots.heatmap(t_plt, z_plt, q_vap_plt .* 1e3, title = "q_vap [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + p9 = Plots.heatmap(t_plt, z_plt, ρq_rim_plt .* 1e3, title = "ρq_rim [g/m3]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + p10 = Plots.heatmap(t_plt, z_plt, ρq_liqonice_plt .* 1e3, title = "ρq_liqonice [g/m3]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + p6 = Plots.heatmap(t_plt, z_plt, N_aer_plt .* 1e-6, title = "N_aer [1/cm3]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis, clims=(0, 100)) + p7 = Plots.heatmap(t_plt, z_plt, N_liq_plt .* 1e-6, title = "N_liq [1/cm3]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + p8 = Plots.heatmap(t_plt, z_plt, N_rai_plt .* 1e-6, title = "N_rai [1/cm3]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + p11 = Plots.heatmap(t_plt, z_plt, N_ice_plt .* 1e-6, title = "N_ice [1/cm3]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + p12 = Plots.heatmap(t_plt, z_plt, B_rim_plt, title = "B_rim [-]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + + #! format: on + p = Plots.plot( + p1, + p2, + p3, + p4, + p5, + p6, + p7, + p8, + p9, + p10, + p11, + p12, + size = (2000.0, 1500.0), + bottom_margin = 30.0 * Plots.PlotMeasures.px, + left_margin = 30.0 * Plots.PlotMeasures.px, + layout = (4, 3), + ) + Plots.png(p, joinpath(path, "timeheight.png")) +end + function plot_timeheight(nc_data_file; output = "output", mixed_phase = true, pysdm = false) path = joinpath(@__DIR__, output) mkpath(path) diff --git a/test/unit_tests/common_unit_test.jl b/test/unit_tests/common_unit_test.jl index a026d2cd..7e026e17 100644 --- a/test/unit_tests/common_unit_test.jl +++ b/test/unit_tests/common_unit_test.jl @@ -8,6 +8,7 @@ params = (common_params, thermo_params, air_params, activation_params) p0m = "Precipitation0M" p1m = "Precipitation1M" p2m = "Precipitation2M" + pp3 = "PrecipitationP3" rf_1 = "CliMA_1M" rf_2 = "KK2000" rf_3 = "B1994" @@ -29,6 +30,7 @@ params = (common_params, thermo_params, air_params, activation_params) precip_1m_6 = CO.get_precipitation_type(p1m, toml_dict, rain_formation_choice = rf_6, sedimentation_choice = st_1) precip_1m_7 = CO.get_precipitation_type(p1m, toml_dict, rain_formation_choice = rf_6, sedimentation_choice = st_2) precip_2m = CO.get_precipitation_type(p2m, toml_dict, rain_formation_choice = rf_7) + precip_p3 = CO.get_precipitation_type(pp3, toml_dict) #test @test_throws Exception CO.get_precipitation_type("_", toml_dict, rain_formation_choice = rf_1) @@ -44,6 +46,7 @@ params = (common_params, thermo_params, air_params, activation_params) @test precip_1m_6 isa CO.Precipitation1M @test precip_1m_7 isa CO.Precipitation1M @test precip_2m isa CO.Precipitation2M + @test precip_p3 isa CO.PrecipitationP3 end @@ -53,18 +56,22 @@ end toml_dict = CP.create_toml_dict(FT) eqm = "EquilibriumMoisture" neqm = "NonEquilibriumMoisture" + mp3 = "MoistureP3" #action moisture_eq = CO.get_moisture_type(eqm, toml_dict) moisture_neq = CO.get_moisture_type(neqm, toml_dict) + moisture_p3 = CO.get_moisture_type(mp3, toml_dict) #test @test CO.EquilibriumMoisture <: CO.AbstractMoistureStyle @test CO.NonEquilibriumMoisture <: CO.AbstractMoistureStyle + @test CO.MoistureP3 <: CO.AbstractMoistureStyle @test_throws Exception CO.get_moisture_type("_", toml_dict) @test moisture_eq isa CO.EquilibriumMoisture @test moisture_neq isa CO.NonEquilibriumMoisture + @test moisture_p3 isa CO.MoistureP3 end @@ -72,7 +79,23 @@ end @test_throws Exception CO.initialise_state(CO.AbstractMoistureStyle(), CO.AbstractPrecipitationStyle(), 0) - initial_profiles = (; ρq_tot = 0, ρq_liq = 0, ρq_ice = 0, ρq_rai = 0, ρq_sno = 0, N_liq = 0, N_rai = 0, N_aer = 0) + initial_profiles = (; + ρq_tot = 0, + ρq_liq = 0, + ρq_ice = 0, + ρq_rai = 0, + ρq_sno = 0, + N_liq = 0, + N_rai = 0, + N_aer = 0, + ρq_rim = 0, + ρq_liqonice = 0, + N_ice = 0, + B_rim = 0, + q_vap = 0, + ρq_vap = 0, + q_rai = 0, + ) state = CO.initialise_state(equil_moist, no_precip, initial_profiles) @test state isa CC.Fields.FieldVector @@ -125,6 +148,20 @@ end @test LA.norm(state.N_rai) == 0 @test LA.norm(state.N_aer) == 0 + state = CO.initialise_state(p3_moist, precip_p3, initial_profiles) + @test state isa CC.Fields.FieldVector + @test LA.norm(state.ρq_tot) == 0 + @test LA.norm(state.ρq_liq) == 0 + @test LA.norm(state.ρq_rai) == 0 + @test LA.norm(state.ρq_ice) == 0 + @test LA.norm(state.ρq_rim) == 0 + @test LA.norm(state.ρq_liqonice) == 0 + @test LA.norm(state.N_ice) == 0 + @test LA.norm(state.B_rim) == 0 + @test LA.norm(state.N_liq) == 0 + @test LA.norm(state.N_rai) == 0 + @test LA.norm(state.N_aer) == 0 + end @@ -139,11 +176,17 @@ end q_ice = 2e-3 / 1.2, q_rai = 1e-3 / 1.2, q_sno = 2e-3 / 1.2, + q_rim = 0.5e-3 / 1.2, + q_liqonice = 0.5e-3 / 1.2, + q_vap = 1e-2 / 1.2, ρq_tot = 15e-3, ρq_liq = 1e-3, ρq_ice = 2e-3, ρq_rai = 1e-3, ρq_sno = 2e-3, + ρq_rim = 0.5e-3, + ρq_liqonice = 0.5e-3, + ρq_vap = 1e-2, p = 101300.0, T = 280.0, θ_dry = 360.0, @@ -152,6 +195,7 @@ end N_rai = 1e4, N_sno = 1e4, N_aer = 1e10, + B_rim = 0.5e-3 / 900, zero = 0.0, ) domain = CC.Domains.IntervalDomain( @@ -203,6 +247,12 @@ end CO.zero_tendencies!(dY) @test LA.norm(dY) == 0 + aux = CO.initialise_aux(FT, ip, params..., 0.0, 0.0, p3_moist, precip_p3) + Y = CO.initialise_state(p3_moist, precip_p3, ip) + dY = similar(Y) + CO.zero_tendencies!(dY) + @test LA.norm(dY) == 0 + @test_throws Exception CO.cloud_sources_tendency!(CO.AbstractMoistureStyle(), dY, Y, aux, 1.0) @test_throws Exception CO.precip_sources_tendency!(CO.AbstractPrecipitationStyle(), dY, Y, aux, 1.0) @@ -240,6 +290,29 @@ end #@test get_value(θ_liq_ice)[1] == TD.liquid_ice_pottemp(thermo_params, get_value(ts)[1]) end + # test p3 precompute_aux_thermo + ms = p3_moist + Y = CO.initialise_state(ms, precip_p3, ip) + aux = CO.initialise_aux(FT, ip, params..., 0.0, 0.0, ms, precip_p3) + CO.precompute_aux_thermo!(ms, Y, aux) + (; ρ_dry, p, T, θ_dry, θ_liq_ice, ts, ρ, ρ_dry) = aux.thermo_variables + (; q_tot, q_liq, q_ice) = aux.microph_variables + for el in aux.thermo_variables + if el != aux.thermo_variables.ts + @test all(isfinite, get_value(el)) + end + end + @test get_value(q_tot)[1] >= 0.0 + @test get_value(q_liq)[1] >= 0.0 + @test get_value(q_ice)[1] >= 0.0 + + @test ρ_dry == ρ .- Y.ρq_tot + @test p == TD.air_pressure.(thermo_params, ts) + @test T == TD.air_temperature.(thermo_params, ts) + @test θ_dry == TD.dry_pottemp.(thermo_params, T, ρ_dry) + #TODO - check Thermodynamics? + #@test get_value(θ_liq_ice)[1] == TD.liquid_ice_pottemp(thermo_params, get_value(ts)[1]) + # test precompute_aux_precip for ps in (precip_1m, precip_2m) Y = CO.initialise_state(equil_moist, ps, ip) @@ -252,6 +325,15 @@ end end end + Y = CO.initialise_state(p3_moist, precip_p3, ip) + aux = CO.initialise_aux(FT, ip, params..., 0.0, 0.0, p3_moist, precip_p3) + CO.precompute_aux_thermo!(p3_moist, Y, aux) + CO.precompute_aux_precip!(precip_p3, Y, aux) + for el in merge(aux.velocities, aux.microph_variables) + @test all(isfinite, get_value(el)) + @test all(get_value(el) .>= FT(0)) + end + # test precompute_aux_moisture_sources Y = CO.initialise_state(nequil_moist, precip_1m, ip) aux = CO.initialise_aux(FT, ip, params..., 0.0, 0.0, nequil_moist, precip_1m) diff --git a/test/unit_tests/k1d_unit_test.jl b/test/unit_tests/k1d_unit_test.jl index 60fa5141..ada4db82 100644 --- a/test/unit_tests/k1d_unit_test.jl +++ b/test/unit_tests/k1d_unit_test.jl @@ -32,11 +32,17 @@ end q_ice = 2e-3 / 1.2, q_rai = 1e-3 / 1.2, q_sno = 2e-3 / 1.2, + q_rim = 0.5e-3 / 1.2, + q_liqonice = 0.5e-3 / 1.2, + q_vap = 1e-2 / 1.2, ρq_tot = 15e-3, ρq_liq = 1e-3, ρq_ice = 2e-3, ρq_rai = 1e-3, ρq_sno = 2e-3, + ρq_rim = 0.5e-3, + ρq_liqonice = 0.5e-3, + ρq_vap = 1e-2, p = 101300.0, T = 280.0, θ_dry = 360.0, @@ -44,7 +50,8 @@ end N_ice = 1e8, N_rai = 1e4, N_sno = 1e4, - N_aer = 5e7, + N_aer = 1e10, + B_rim = 0.5e-3 / 900, zero = 0.0, ) domain = CC.Domains.IntervalDomain( @@ -74,6 +81,11 @@ end @test aux.cloud_sources === nothing end end + ms = p3_moist + ps = precip_p3 + aux = K1D.initialise_aux(FT, ip, params..., 0.0, 0.0, face_space, ms, ps) + @test aux isa NamedTuple + @test LA.norm(aux.precip_sources) == 0 end @testset "advection_tendency" begin @@ -107,6 +119,102 @@ end end end +@testset "advection_tendency p3" begin + + space, face_space = K1D.make_function_space(FT, 0, 100, 5) + coord = CC.Fields.coordinate_field(space) + ice_start = false + _magnitude = Float64(0.5) + _q_flux = Float64(0.65e-4) + _N_flux = Float64(40000) + _F_rim = Float64(0.2) + _F_liq = Float64(0.2) + _ρ_r_init = Float64(900) + init = map( + coord -> CO.p3_initial_condition( + FT, + kid_params, + thermo_params, + coord.z; + _q_init = _q_flux, + _N_init = _N_flux, + _F_rim = _F_rim, + _F_liq = _F_liq, + _ρ_r = _ρ_r_init, + z_top = FT(3000), + ice_start = ice_start, + ), + coord, + ) + t = 13.0 + + # eq + aux = K1D.initialise_aux(FT, init, params..., 0.0, 0.0, face_space, p3_moist, precip_p3) + Y = CO.initialise_state(p3_moist, precip_p3, init) + dY = Y / 10 + @test_throws Exception K1D.advection_tendency!(CO.AbstractMoistureStyle(), dY, Y, aux, t) + @test_throws Exception K1D.advection_tendency!(CO.AbstractPrecipitationStyle(), dY, Y, aux, t) + + ms = p3_moist + ps = precip_p3 + aux = K1D.initialise_aux(FT, init, params..., 0.0, 0.0, face_space, ms, ps) + Y = CO.initialise_state(ms, ps, init) + dY = Y / 10 + ρw = 0.0 + + @. aux.prescribed_velocity.ρw = CC.Geometry.WVector.(ρw) + aux.prescribed_velocity.ρw0 = ρw + K1D.advection_tendency!(ms, dY, Y, aux, t) + K1D.advection_tendency!(ps, dY, Y, aux, t) + @test dY >= Y / 10 + + # same tests with ice_start = true + ice_start = true + _magnitude = Float64(0.5) + _q_flux = Float64(0.65e-4) + _N_flux = Float64(40000) + _F_rim = Float64(0.2) + _F_liq = Float64(0.2) + _ρ_r_init = Float64(900) + init = map( + coord -> CO.p3_initial_condition( + FT, + kid_params, + thermo_params, + coord.z; + _q_init = _q_flux, + _N_init = _N_flux, + _F_rim = _F_rim, + _F_liq = _F_liq, + _ρ_r = _ρ_r_init, + z_top = FT(3000), + ice_start = ice_start, + ), + coord, + ) + t = 13.0 + + # eq + aux = K1D.initialise_aux(FT, init, params..., 0.0, 0.0, face_space, p3_moist, precip_p3) + Y = CO.initialise_state(p3_moist, precip_p3, init) + dY = Y / 10 + @test_throws Exception K1D.advection_tendency!(CO.AbstractMoistureStyle(), dY, Y, aux, t) + @test_throws Exception K1D.advection_tendency!(CO.AbstractPrecipitationStyle(), dY, Y, aux, t) + + ms = p3_moist + ps = precip_p3 + aux = K1D.initialise_aux(FT, init, params..., 0.0, 0.0, face_space, ms, ps) + Y = CO.initialise_state(ms, ps, init) + dY = Y / 10 + ρw = 0.0 + + @. aux.prescribed_velocity.ρw = CC.Geometry.WVector.(ρw) + aux.prescribed_velocity.ρw0 = ρw + K1D.advection_tendency!(ms, dY, Y, aux, t) + K1D.advection_tendency!(ps, dY, Y, aux, t) + @test dY >= Y / 10 +end + @testset "aerosol activation for 2M schemes" begin #setup _ip = (; diff --git a/test/unit_tests/unit_test.jl b/test/unit_tests/unit_test.jl index 84e99d27..62f3f0fb 100644 --- a/test/unit_tests/unit_test.jl +++ b/test/unit_tests/unit_test.jl @@ -30,9 +30,22 @@ thermo_params = create_thermodynamics_parameters(toml_dict) kid_params = create_kid_parameters(toml_dict) air_params = CMP.AirProperties(toml_dict) activation_params = CMP.AerosolActivationParameters(toml_dict) +p3 = CMP.ParametersP3(FT) +Chen2022 = CMP.Chen2022VelType(FT) +# (use default boundary condition) +p3_boundary_condition = (; + ice_start = false, + _magnitude = Float64(0.5), + _q_flux = Float64(0.65e-4), + _N_flux = Float64(40000), + _F_rim = Float64(0.2), + _F_liq = Float64(0.2), + _ρ_r_init = Float64(900), +) # ... for cloud condensate options ... equil_moist = CO.EquilibriumMoisture() nequil_moist = CO.NonEquilibriumMoisture(CMP.CloudLiquid(toml_dict), CMP.CloudIce(toml_dict)) +p3_moist = CO.MoistureP3() # ... and precipitation options no_precip = CO.NoPrecipitation() precip_0m = CO.Precipitation0M(CMP.Parameters0M(toml_dict)) @@ -46,6 +59,7 @@ precip_1m = CO.Precipitation1M( CMP.Blk1MVelType(toml_dict), ) precip_2m = CO.Precipitation2M(CMP.SB2006(toml_dict), CMP.SB2006VelType(toml_dict)) +precip_p3 = CO.PrecipitationP3(p3, Chen2022, CMP.SB2006(toml_dict), p3_boundary_condition) # common unit tests include("./common_unit_test.jl")