From 18ae71df5bcdc0ea790bea06339dc794e0a780dc Mon Sep 17 00:00:00 2001 From: Rowan Orlijan-Rhyne Date: Tue, 30 Jul 2024 16:33:22 -0700 Subject: [PATCH 01/22] add p3 advection example --- src/Common/Common.jl | 1 + src/Common/equation_types.jl | 8 +- src/Common/helper_functions.jl | 8 + src/Common/initial_condition.jl | 97 ++++++ src/Common/netcdf_io.jl | 7 + src/Common/ode_utils.jl | 15 +- src/Common/tendency.jl | 78 ++++- src/K1DModel/tendency.jl | 197 ++++++++++- .../KiD_driver/parse_commandline.jl | 4 +- .../KiD_driver/run_KiD_simulation.jl | 29 +- test/plotting_utils.jl | 329 ++++++++++++------ test/unit_tests/common_unit_test.jl | 48 ++- test/unit_tests/unit_test.jl | 4 + 13 files changed, 705 insertions(+), 120 deletions(-) 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..cb392f4d 100644 --- a/src/Common/equation_types.jl +++ b/src/Common/equation_types.jl @@ -52,9 +52,9 @@ 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} <: AbstractPrecipitationStyle + p3_params::P3 + Chen2022::CH + sb2006::PT end struct CloudyPrecip <: AbstractPrecipitationStyle end diff --git a/src/Common/helper_functions.jl b/src/Common/helper_functions.jl index e36da5e7..1cb397d5 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 @@ -77,6 +79,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) else error("Invalid precipitation choice: $precipitation_choice") end diff --git a/src/Common/initial_condition.jl b/src/Common/initial_condition.jl index 1658b82f..bab74671 100644 --- a/src/Common/initial_condition.jl +++ b/src/Common/initial_condition.jl @@ -281,3 +281,100 @@ 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}, thermo_params, q_init, N_init, z; F_rim, F_liq, z_top, ice_start) where {FT} + + _ρ_r::FT = FT(900) # fixed initial ρ_r + + # 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) + 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 + q_tot::FT = q_ice + q_liq + q_rai + + # 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 + + # 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) + ρq_vap::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, + 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..324720dd 100644 --- a/src/Common/netcdf_io.jl +++ b/src/Common/netcdf_io.jl @@ -23,11 +23,18 @@ 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", :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", ), ) FT = Float64 diff --git a/src/Common/ode_utils.jl b/src/Common/ode_utils.jl index ab1fdcfd..95f5a79a 100644 --- a/src/Common/ode_utils.jl +++ b/src/Common/ode_utils.jl @@ -76,6 +76,7 @@ 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, @@ -116,7 +117,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 +182,18 @@ 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, ) velocities = (; term_vel_rai = copy(ip.zero), @@ -199,6 +207,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 +226,11 @@ function initialise_aux( copy(ip.zero), copy(ip.zero), copy(ip.zero), + copy(ip.zero), ), ) activation_sources = nothing - @assert moisture isa NonEquilibrumMoisture elseif precip isa CloudyPrecip microph_variables = (; q_tot = ip.q_tot, @@ -241,7 +250,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 = (; diff --git a/src/Common/tendency.jl b/src/Common/tendency.jl index 73aee173..c5dbe667 100644 --- a/src/Common/tendency.jl +++ b/src/Common/tendency.jl @@ -47,8 +47,35 @@ 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_tot, q_liq, q_ice, q_rim, q_liqonice, ρq_tot, ρq_liq, ρq_ice, ρq_rim, ρq_liqonice) = aux.microph_variables + + @. ρ = ρ_dry + Y.ρq_tot + @. ρq_tot = Y.ρq_tot + @. ρq_rim = Y.ρq_rim + @. ρq_ice = Y.ρq_ice + @. ρq_liqonice = Y.ρq_liqonice + + @. q_tot = q_(Y.ρq_tot, ρ) + @. q_liq = q_(Y.ρq_liq, ρ) + @. q_ice = q_(Y.ρq_ice, ρ) + @. q_rim = q_(Y.ρq_rim, ρ) + @. q_liqonice = q_(Y.ρq_liqonice, ρ) + + @. 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 = θ_dry #TD.liquid_ice_pottemp(thermo_params, ts) - prevents error... +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 @@ -62,7 +89,7 @@ end @. 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) + @. θ_liq_ice = θ_dry #TD.liquid_ice_pottemp(thermo_params, ts) - prevents error... end @inline function precompute_aux_thermo!(::CloudyMoisture, Y, aux) @@ -115,19 +142,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 # helper function for precomputing aux precip for cloudy @@ -477,6 +538,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..b453a1d9 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)) + + if Bool(aux.kid_params.qtot_flux_correction) + 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) + end + + 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,179 @@ 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) = aux.microph_variables + @. ρq_ice = Y.ρq_ice + @. ρq_liq = Y.ρq_liq + @. ρq_rim = Y.ρq_rim + @. ρq_liqonice = Y.ρq_liqonice + @. N_ice = Y.N_ice + @. B_rim = Y.B_rim + + # choose flux characteristics + _q_flux = FT(0.65e-4) + _N_flux = FT(40000) + _F_rim = FT(0.2) + _F_liq = FT(0.2) + _ρ_r_init = FT(900) + + _flux = FT(0.5) + + 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(_flux * _q_flux)), + ) + ∂q_rim = CC.Operators.DivergenceF2C( + bottom = CC.Operators.Extrapolate(), + top = CC.Operators.SetValue(CC.Geometry.WVector(_flux * _F_rim * (1 - _F_liq) * _q_flux)), + ) + ∂q_liqonice = CC.Operators.DivergenceF2C( + bottom = CC.Operators.Extrapolate(), + top = CC.Operators.SetValue(CC.Geometry.WVector(_flux * _F_liq * _q_flux)), + ) + ∂B_rim = CC.Operators.DivergenceF2C( + bottom = CC.Operators.Extrapolate(), + top = CC.Operators.SetValue(CC.Geometry.WVector(_flux * _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(_flux * _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_tot = CC.Operators.DivergenceF2C( + bottom = CC.Operators.Extrapolate(), + top = CC.Operators.SetValue(CC.Geometry.WVector(_flux * _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/experiments/KiD_driver/parse_commandline.jl b/test/experiments/KiD_driver/parse_commandline.jl index e973a5d2..4051d3f2 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 = "NonEquilibriumMoisture" "--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" "--rain_formation_scheme_choice" diff --git a/test/experiments/KiD_driver/run_KiD_simulation.jl b/test/experiments/KiD_driver/run_KiD_simulation.jl index f9bdd044..0cb3fdb5 100644 --- a/test/experiments/KiD_driver/run_KiD_simulation.jl +++ b/test/experiments/KiD_driver/run_KiD_simulation.jl @@ -104,6 +104,33 @@ function run_KiD_simulation(::Type{FT}, opts) where {FT} ), coord, ) + elseif precipitation_choice == "PrecipitationP3" + # commenting here the configuration which has been used for p3: + # config["precipitation_choice"] = "PrecipitationP3" + # config["moisture_choice"] = "MoistureP3" + # config["n_elem"] = 24 + # config["z_max"] = 3000 + # config["t_end"] = 75 + # config["w1"] = 0 + # config["rv_0"] = 0 + # config["rv_1"] = 0 + # config["rv_2"] = 0 + # config["dt"] = Float64(0.5) + cloudy_params = nothing + init = map( + coord -> CO.p3_initial_condition( + FT, + thermo_params, + FT(0.65e-4), + FT(40000), + coord.z, + F_rim = FT(0.2), + F_liq = FT(0.2), + z_top = FT(opts["z_max"]), + ice_start = false, + ), + coord, + ) else cloudy_params = nothing init = map( @@ -164,7 +191,7 @@ 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(z_centers, solver, aux, moisture, precip, K1D, output = plot_folder) - plot_timeheight(string("experiments/KiD_driver/", output_folder, "/Output.nc"), output = plot_folder) + plot_timeheight(string("experiments/KiD_driver/", output_folder, "/Output.nc"), precip, output = plot_folder) end return diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index 5872f5cd..c89024b4 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,53 +94,78 @@ 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), + ) + 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 - 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) + p = Plots.plot( + p1, + p2, + p3, + p4, + p5, + p6, + p7, + p8, + p9, + p10, + p11, + size = (1800.0, 1200.0), + bottom_margin = 40.0 * Plots.PlotMeasures.px, + left_margin = 80.0 * Plots.PlotMeasures.px, + ) 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")) end @@ -151,6 +184,15 @@ function plot_animation(z_centers, solver, aux, moisture, precip, KiD; output = 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) if moisture isa CO.EquilibriumMoisture @@ -162,6 +204,17 @@ function plot_animation(z_centers, solver, aux, moisture, precip, KiD; output = elseif moisture isa CO.CloudyMoisture q_tot[:, i] = (parent(u.moments.:2) .+ parent(u.ρq_vap)) ./ ρ .* 1e3 q_liq[:, i] = parent(u.moments.:2) ./ ρ .* 1e3 + elseif precip isa CO.PrecipitationP3 + ρ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 end if precip isa CO.Precipitation1M @@ -195,15 +248,137 @@ function plot_animation(z_centers, solver, aux, moisture, precip, KiD; output = 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/kg]", maximum(q_tot)) - p2 = plot_data(q_liq[:, i], "q_liq [g/kg]", maximum(q_liq), title) - p3 = plot_data(N_liq[:, i], "N_liq [1/cm^3]", maximum(N_liq)) - p4 = plot_data(q_rai[:, i], "q_rai [g/kg]", maximum(q_rai)) - p5 = plot_data(N_rai[:, i], "N_rai [1/cm^3]", maximum(N_rai)) - p6 = plot_data(q_ice[:, i], "q_ice [g/kg]", maximum(q_ice)) - p7 = plot_data(q_sno[:, i], "q_sno [g/kg]", maximum(q_sno)) - - plot( + + if precip isa CO.PrecipitationP3 + 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), + ) + else + p1 = plot_data(q_tot[:, i], "q_tot [g/kg]", maximum(q_tot)) + p2 = plot_data(q_liq[:, i], "q_liq [g/kg]", maximum(q_liq), title) + p3 = plot_data(N_liq[:, i], "N_liq [1/cm^3]", maximum(N_liq)) + p4 = plot_data(q_rai[:, i], "q_rai [g/kg]", maximum(q_rai)) + p5 = plot_data(N_rai[:, i], "N_rai [1/cm^3]", maximum(N_rai)) + p6 = plot_data(q_ice[:, i], "q_ice [g/kg]", maximum(q_ice)) + p7 = plot_data(q_sno[:, i], "q_sno [g/kg]", maximum(q_sno)) + plot( + p1, + p2, + p3, + p4, + p5, + p6, + p7, + 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 = (3, 3), + ) + end + end + + Plots.mp4(anim, joinpath(path, "animation.mp4"), fps = 10) +end + +function plot_timeheight(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"]) + if precip isa CO.PrecipitationP3 + 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_sno_plt = Array(ds.group["profiles"]["q_sno"]) + 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_sno_plt .* 1e3, title = "q_sno [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), + ) + else + 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_rai_plt = Array(ds.group["profiles"]["q_rai"]) + q_sno_plt = Array(ds.group["profiles"]["q_sno"]) + 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"]) + #! 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/kg]", 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_sno_plt .* 1e3, title = "q_sno [g/kg]", 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) + #! format: on + p = Plots.plot( p1, p2, p3, @@ -211,57 +386,13 @@ function plot_animation(z_centers, solver, aux, moisture, precip, KiD; output = p5, p6, p7, - size = (1800.0, 1500.0), + p8, + size = (2000.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 = (3, 3), ) end - - Plots.mp4(anim, joinpath(path, "animation.mp4"), fps = 10) -end - -function plot_timeheight(nc_data_file; 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_rai_plt = Array(ds.group["profiles"]["q_rai"]) - q_sno_plt = Array(ds.group["profiles"]["q_sno"]) - 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"]) - #! 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/kg]", 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_sno_plt .* 1e3, title = "q_sno [g/kg]", 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) - #! format: on - p = Plots.plot( - p1, - p2, - p3, - p4, - p5, - p6, - p7, - p8, - size = (2000.0, 1500.0), - bottom_margin = 30.0 * Plots.PlotMeasures.px, - left_margin = 30.0 * Plots.PlotMeasures.px, - layout = (3, 3), - ) Plots.png(p, joinpath(path, "timeheight.png")) end function plot_timeheight_no_ice_snow(nc_data_file; output = "output") diff --git a/test/unit_tests/common_unit_test.jl b/test/unit_tests/common_unit_test.jl index a026d2cd..8cf1c24c 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,20 @@ 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, + ) state = CO.initialise_state(equil_moist, no_precip, initial_profiles) @test state isa CC.Fields.FieldVector @@ -125,6 +145,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 +173,15 @@ 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_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, p = 101300.0, T = 280.0, θ_dry = 360.0, @@ -152,7 +190,7 @@ end N_rai = 1e4, N_sno = 1e4, N_aer = 1e10, - zero = 0.0, + B_rim = zero = 0.0, ) domain = CC.Domains.IntervalDomain( CC.Geometry.ZPoint{FT}(0), @@ -203,6 +241,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) diff --git a/test/unit_tests/unit_test.jl b/test/unit_tests/unit_test.jl index 84e99d27..8172511d 100644 --- a/test/unit_tests/unit_test.jl +++ b/test/unit_tests/unit_test.jl @@ -30,9 +30,12 @@ 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) # ... 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 +49,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)) # common unit tests include("./common_unit_test.jl") From 0dd332d2f18daf3f41b84412c241429b63c9aec1 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Wed, 21 Aug 2024 13:41:37 -0700 Subject: [PATCH 02/22] fixes + test --- .buildkite/pipeline.yml | 4 ++ src/Common/equation_types.jl | 3 +- src/Common/helper_functions.jl | 3 +- src/Common/initial_condition.jl | 34 +++++++--- src/Common/netcdf_io.jl | 3 + src/Common/ode_utils.jl | 27 ++++++++ src/Common/tendency.jl | 16 +++-- src/K1DModel/tendency.jl | 41 ++++++------ .../KiD_driver/parse_commandline.jl | 12 ++++ .../KiD_driver/run_KiD_simulation.jl | 32 +++++----- test/plotting_utils.jl | 4 +- test/unit_tests/common_unit_test.jl | 31 ++++++++- test/unit_tests/k1d_unit_test.jl | 63 ++++++++++++++++++- test/unit_tests/unit_test.jl | 12 +++- 14 files changed, 227 insertions(+), 58 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index de0eeaee..43e2e7ef 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -126,3 +126,7 @@ steps: - label: ":partly_sunny: Experiment: 1D with Cloudy" command: "julia --color=yes --project=test test/experiments/KiD_driver/KiD_driver.jl --moisture_choice=CloudyMoisture --precipitation_choice=CloudyPrecip" artifact_paths: "test/experiments/KiD_driver/Output_CloudyMoisture_CloudyPrecip/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=Float64(0.5)" + artifact_paths: "test/experiments/KiD_driver/Output_MoistureP3_PrecipitationP3/figures/*" \ No newline at end of file diff --git a/src/Common/equation_types.jl b/src/Common/equation_types.jl index cb392f4d..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{P3, CH, PT} <: AbstractPrecipitationStyle +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 1cb397d5..cef19ead 100644 --- a/src/Common/helper_functions.jl +++ b/src/Common/helper_functions.jl @@ -24,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() @@ -84,7 +85,7 @@ function get_precipitation_type( p3_params = CMP.ParametersP3(FT) Chen2022 = CMP.Chen2022VelType(FT) sb2006 = CMP.SB2006(toml_dict) - precip = PrecipitationP3(p3_params, Chen2022, sb2006) + 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 bab74671..81c7edcf 100644 --- a/src/Common/initial_condition.jl +++ b/src/Common/initial_condition.jl @@ -282,9 +282,24 @@ 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}, thermo_params, q_init, N_init, z; F_rim, F_liq, z_top, ice_start) where {FT} +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} - _ρ_r::FT = FT(900) # fixed initial ρ_r + # 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 @@ -294,10 +309,10 @@ function p3_initial_condition(::Type{FT}, thermo_params, q_init, N_init, z; F_ri ice_bump(χ, z) = χ * cos(z - (z_top - 0.5 * _z_band)) # P3 state variables (minus B_rim -- computed below) - 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) + 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) @@ -305,8 +320,8 @@ function p3_initial_condition(::Type{FT}, thermo_params, q_init, N_init, z; F_ri q_rai::FT = FT(0) N_rai::FT = FT(0) - # q_tot - single p3 category + 2M categories - q_tot::FT = q_ice + q_liq + q_rai + # 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 @@ -332,6 +347,7 @@ function p3_initial_condition(::Type{FT}, thermo_params, q_init, N_init, z; F_ri ρ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 @@ -340,7 +356,6 @@ function p3_initial_condition(::Type{FT}, thermo_params, q_init, N_init, z; F_ri # unused quantities: q_sno::FT = FT(0) ρq_sno::FT = FT(0) - ρq_vap::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, ρ) @@ -370,6 +385,7 @@ function p3_initial_condition(::Type{FT}, thermo_params, q_init, N_init, z; F_ri q_sno, q_rim, q_liqonice, + q_vap, B_rim, N_liq, N_ice, diff --git a/src/Common/netcdf_io.jl b/src/Common/netcdf_io.jl index 324720dd..272ebad6 100644 --- a/src/Common/netcdf_io.jl +++ b/src/Common/netcdf_io.jl @@ -27,6 +27,7 @@ function NetCDFIO_Stats( :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", @@ -35,6 +36,8 @@ function NetCDFIO_Stats( :ρ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 95f5a79a..0f4b4b0a 100644 --- a/src/Common/ode_utils.jl +++ b/src/Common/ode_utils.jl @@ -82,6 +82,9 @@ function initialise_state(::MoistureP3, ::PrecipitationP3, initial_profiles) 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) @@ -194,6 +197,8 @@ function initialise_aux( ρ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), @@ -229,6 +234,26 @@ function initialise_aux( 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 elseif precip isa CloudyPrecip @@ -271,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 c5dbe667..46025869 100644 --- a/src/Common/tendency.jl +++ b/src/Common/tendency.jl @@ -53,24 +53,30 @@ end (; thermo_params) = aux (; ts, ρ, ρ_dry, p, T, θ_dry, θ_liq_ice) = aux.thermo_variables - (; q_tot, q_liq, q_ice, q_rim, q_liqonice, ρq_tot, ρq_liq, ρq_ice, ρq_rim, ρq_liqonice) = aux.microph_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_tot = 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_tot = q_(Y.ρq_tot, ρ) @. 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 = θ_dry #TD.liquid_ice_pottemp(thermo_params, ts) - prevents error... + @. θ_liq_ice = TD.liquid_ice_pottemp(thermo_params, ts) + end @inline function precompute_aux_thermo!(::NonEquilibriumMoisture, Y, aux) @@ -89,7 +95,7 @@ end @. 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 = θ_dry #TD.liquid_ice_pottemp(thermo_params, ts) - prevents error... + @. θ_liq_ice = TD.liquid_ice_pottemp(thermo_params, ts) end @inline function precompute_aux_thermo!(::CloudyMoisture, Y, aux) diff --git a/src/K1DModel/tendency.jl b/src/K1DModel/tendency.jl index b453a1d9..480ac5c1 100644 --- a/src/K1DModel/tendency.jl +++ b/src/K1DModel/tendency.jl @@ -280,11 +280,11 @@ end 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)) - if Bool(aux.kid_params.qtot_flux_correction) - 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) - end + 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 @@ -424,22 +424,20 @@ end # 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) = aux.microph_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_liq = Y.ρq_liq @. ρq_rim = Y.ρq_rim @. ρq_liqonice = Y.ρq_liqonice @. N_ice = Y.N_ice @. B_rim = Y.B_rim # choose flux characteristics - _q_flux = FT(0.65e-4) - _N_flux = FT(40000) - _F_rim = FT(0.2) - _F_liq = FT(0.2) - _ρ_r_init = FT(900) - - _flux = FT(0.5) + (; 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() @@ -447,23 +445,23 @@ end # for each state variable ∂q_ice = CC.Operators.DivergenceF2C( bottom = CC.Operators.Extrapolate(), - top = CC.Operators.SetValue(CC.Geometry.WVector(_flux * _q_flux)), + 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(_flux * _F_rim * (1 - _F_liq) * _q_flux)), + 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(_flux * _F_liq * _q_flux)), + 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(_flux * _F_rim * (1 - _F_liq) * _q_flux / _ρ_r_init)), + 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(_flux * _N_flux)), + top = CC.Operators.SetValue(CC.Geometry.WVector(_magnitude * _N_flux)), ) # apply divergence @@ -540,7 +538,7 @@ end CC.Geometry.WVector(If(aux.velocities.term_vel_N_rai) * FT(-1)) ) * If(Y.N_rai), ) - @. dY.ρq_rai += + @. dY.q_rai += -∂( ( aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ) + @@ -568,9 +566,12 @@ end # 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(_flux * _q_flux)), + top = CC.Operators.SetValue(CC.Geometry.WVector(_magnitude * _q_flux)), ) @. dY.ρq_tot += ∂q_tot( diff --git a/test/experiments/KiD_driver/parse_commandline.jl b/test/experiments/KiD_driver/parse_commandline.jl index 4051d3f2..16537983 100644 --- a/test/experiments/KiD_driver/parse_commandline.jl +++ b/test/experiments/KiD_driver/parse_commandline.jl @@ -141,6 +141,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 0cb3fdb5..25b5725e 100644 --- a/test/experiments/KiD_driver/run_KiD_simulation.jl +++ b/test/experiments/KiD_driver/run_KiD_simulation.jl @@ -34,6 +34,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) @@ -77,6 +82,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" @@ -105,29 +111,21 @@ function run_KiD_simulation(::Type{FT}, opts) where {FT} coord, ) elseif precipitation_choice == "PrecipitationP3" - # commenting here the configuration which has been used for p3: - # config["precipitation_choice"] = "PrecipitationP3" - # config["moisture_choice"] = "MoistureP3" - # config["n_elem"] = 24 - # config["z_max"] = 3000 - # config["t_end"] = 75 - # config["w1"] = 0 - # config["rv_0"] = 0 - # config["rv_1"] = 0 - # config["rv_2"] = 0 - # config["dt"] = Float64(0.5) 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, - FT(0.65e-4), - FT(40000), - coord.z, - F_rim = FT(0.2), - F_liq = FT(0.2), + 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 = false, + ice_start = ice_start, ), coord, ) diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index c89024b4..847bebc6 100644 --- a/test/plotting_utils.jl +++ b/test/plotting_utils.jl @@ -320,7 +320,7 @@ function plot_timeheight(nc_data_file, precip; output = "output") ρ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_sno_plt = Array(ds.group["profiles"]["q_sno"]) + 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"]) @@ -331,7 +331,7 @@ function plot_timeheight(nc_data_file, precip; output = "output") 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_sno_plt .* 1e3, title = "q_sno [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)) diff --git a/test/unit_tests/common_unit_test.jl b/test/unit_tests/common_unit_test.jl index 8cf1c24c..a2a3f413 100644 --- a/test/unit_tests/common_unit_test.jl +++ b/test/unit_tests/common_unit_test.jl @@ -92,6 +92,9 @@ end ρ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) @@ -175,6 +178,7 @@ end 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, @@ -182,6 +186,7 @@ end ρ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, @@ -190,7 +195,8 @@ end N_rai = 1e4, N_sno = 1e4, N_aer = 1e10, - B_rim = zero = 0.0, + B_rim = 0.5e-3 / 900, + zero = 0.0, ) domain = CC.Domains.IntervalDomain( CC.Geometry.ZPoint{FT}(0), @@ -284,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) diff --git a/test/unit_tests/k1d_unit_test.jl b/test/unit_tests/k1d_unit_test.jl index 60fa5141..ea4d7c44 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,55 @@ 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 +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 8172511d..62f3f0fb 100644 --- a/test/unit_tests/unit_test.jl +++ b/test/unit_tests/unit_test.jl @@ -32,6 +32,16 @@ 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)) @@ -49,7 +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)) +precip_p3 = CO.PrecipitationP3(p3, Chen2022, CMP.SB2006(toml_dict), p3_boundary_condition) # common unit tests include("./common_unit_test.jl") From f27107af2f1a68226e788c939488fbb401efe00c Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Wed, 21 Aug 2024 14:16:34 -0700 Subject: [PATCH 03/22] parsing fix --- .../KiD_driver/run_KiD_simulation.jl | 20 +++++++++++-------- test/plotting_utils.jl | 3 ++- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/test/experiments/KiD_driver/run_KiD_simulation.jl b/test/experiments/KiD_driver/run_KiD_simulation.jl index 4e9fc001..81d58daa 100644 --- a/test/experiments/KiD_driver/run_KiD_simulation.jl +++ b/test/experiments/KiD_driver/run_KiD_simulation.jl @@ -192,15 +192,19 @@ 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) if precip isa PrecipitationP3 - plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, output = plot_folder) - plot_timeheight_p3(string("experiments/KiD_driver/", output_folder, "/Output.nc"), precip, output = plot_folder) + plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, output = 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, - ) + 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 diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index 91cf0139..99a675be 100644 --- a/test/plotting_utils.jl +++ b/test/plotting_utils.jl @@ -229,6 +229,7 @@ function plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, output N_liq[:, i] = parent(u.moments.:1) ./ 1e6 N_rai[:, i] = parent(u.moments.:4) ./ 1e6 end + end function plot_data(data, data_label, max_val, title = "") return Plots.plot( @@ -304,7 +305,7 @@ function plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, output Plots.mp4(anim, joinpath(path, "animation.mp4"), fps = 10) end - + function plot_animation(nc_data_file; output = "output") path = joinpath(@__DIR__, output) From becf4de035cb6b8a06b634c6e7e394803dc3b214 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Wed, 21 Aug 2024 14:27:35 -0700 Subject: [PATCH 04/22] test fix --- test/plotting_utils.jl | 33 ++++++++++++++++----------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index 99a675be..fd9492a6 100644 --- a/test/plotting_utils.jl +++ b/test/plotting_utils.jl @@ -147,24 +147,23 @@ function plot_final_aux_profiles(z_centers, aux, precip; output = "output") 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) + p = Plots.plot( + p1, + p2, + p3, + p4, + p5, + p6, + p7, + p8, + p9, + p10, + p11, + size = (1800.0, 1200.0), + bottom_margin = 40.0 * Plots.PlotMeasures.px, + left_margin = 80.0 * Plots.PlotMeasures.px, + ) end - - p = Plots.plot( - p1, - p2, - p3, - p4, - p5, - p6, - p7, - p8, - p9, - p10, - p11, - size = (1800.0, 1200.0), - bottom_margin = 40.0 * Plots.PlotMeasures.px, - left_margin = 80.0 * Plots.PlotMeasures.px, - ) end Plots.png(p, joinpath(path, "final_aux_profiles.png")) end From 2b15a389839201b948a409e5d702bad248035df1 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Wed, 21 Aug 2024 14:47:54 -0700 Subject: [PATCH 05/22] more test fixes --- test/experiments/KiD_driver/run_KiD_simulation.jl | 2 +- test/plotting_utils.jl | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/test/experiments/KiD_driver/run_KiD_simulation.jl b/test/experiments/KiD_driver/run_KiD_simulation.jl index 81d58daa..dfce4b51 100644 --- a/test/experiments/KiD_driver/run_KiD_simulation.jl +++ b/test/experiments/KiD_driver/run_KiD_simulation.jl @@ -191,7 +191,7 @@ 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) - if precip isa PrecipitationP3 + if precip isa CO.PrecipitationP3 plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, output = plot_folder) plot_timeheight_p3( string("experiments/KiD_driver/", output_folder, "/Output.nc"), diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index fd9492a6..1f60e590 100644 --- a/test/plotting_utils.jl +++ b/test/plotting_utils.jl @@ -117,6 +117,7 @@ function plot_final_aux_profiles(z_centers, aux, precip; output = "output") 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]") @@ -163,9 +164,9 @@ function plot_final_aux_profiles(z_centers, aux, precip; output = "output") 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 - Plots.png(p, joinpath(path, "final_aux_profiles.png")) end function plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, output = plot_folder) From 54eb2222ca1b8c66c60c15cb02c3a40bdcfc5ab4 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Wed, 21 Aug 2024 14:56:10 -0700 Subject: [PATCH 06/22] test fix again... oops --- test/plotting_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index 1f60e590..9dec062e 100644 --- a/test/plotting_utils.jl +++ b/test/plotting_utils.jl @@ -148,6 +148,7 @@ function plot_final_aux_profiles(z_centers, aux, precip; output = "output") 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, @@ -165,7 +166,6 @@ function plot_final_aux_profiles(z_centers, aux, precip; output = "output") left_margin = 80.0 * Plots.PlotMeasures.px, ) Plots.png(p, joinpath(path, "final_aux_profiles.png")) - end end end From 62525224c57c9b8251c76a434fe8b89c6bd36047 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Wed, 21 Aug 2024 15:11:33 -0700 Subject: [PATCH 07/22] bad day for my tests --- test/plotting_utils.jl | 33 ++++++++++++++++----------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index 9dec062e..d414d8f6 100644 --- a/test/plotting_utils.jl +++ b/test/plotting_utils.jl @@ -149,23 +149,22 @@ function plot_final_aux_profiles(z_centers, aux, precip; output = "output") 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, - p11, - 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")) + 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 From 279bcb4bdea3f3bcc0dbc4b5692e4cff3bdae7c5 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Wed, 21 Aug 2024 16:10:53 -0700 Subject: [PATCH 08/22] buildkite command line fix --- .buildkite/pipeline.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 31d5fd72..689d8b74 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -136,6 +136,6 @@ steps: 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=Float64(0.5)" + 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/*" From 26d694b26312cd385203a953f4701bf4ada529eb Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Thu, 22 Aug 2024 09:20:50 -0700 Subject: [PATCH 09/22] and again --- .buildkite/pipeline.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 689d8b74..5786f2c7 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -136,6 +136,6 @@ steps: 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" + 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/*" From 91ae4464153fb066c5128b586bd9b05801c417c6 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Thu, 22 Aug 2024 20:05:51 -0700 Subject: [PATCH 10/22] fingers crossed --- src/Common/initial_condition.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Common/initial_condition.jl b/src/Common/initial_condition.jl index 81c7edcf..e739875c 100644 --- a/src/Common/initial_condition.jl +++ b/src/Common/initial_condition.jl @@ -309,6 +309,7 @@ function p3_initial_condition( ice_bump(χ, z) = χ * cos(z - (z_top - 0.5 * _z_band)) # P3 state variables (minus B_rim -- computed below) + # (computed 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 From 70a5b360cf63efa7d9aedd3bbc089642d27e518f Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Fri, 23 Aug 2024 11:28:00 -0700 Subject: [PATCH 11/22] CM release --- Project.toml | 2 +- test/Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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/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" From 507bb6666a893045cb92e04917ba4503c796ac48 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Fri, 23 Aug 2024 13:41:27 -0700 Subject: [PATCH 12/22] buildkite --- src/Common/initial_condition.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Common/initial_condition.jl b/src/Common/initial_condition.jl index e739875c..ccbae713 100644 --- a/src/Common/initial_condition.jl +++ b/src/Common/initial_condition.jl @@ -309,7 +309,7 @@ function p3_initial_condition( ice_bump(χ, z) = χ * cos(z - (z_top - 0.5 * _z_band)) # P3 state variables (minus B_rim -- computed below) - # (computed as kg/kg before converting to kg/m3) + # (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 From d7c3338b9493c882af61d7e824b90cf72fc22906 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Fri, 23 Aug 2024 15:13:27 -0700 Subject: [PATCH 13/22] oops --- test/experiments/KiD_driver/run_KiD_simulation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/experiments/KiD_driver/run_KiD_simulation.jl b/test/experiments/KiD_driver/run_KiD_simulation.jl index dfce4b51..c897844b 100644 --- a/test/experiments/KiD_driver/run_KiD_simulation.jl +++ b/test/experiments/KiD_driver/run_KiD_simulation.jl @@ -192,7 +192,7 @@ 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) if precip isa CO.PrecipitationP3 - plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, output = plot_folder) + plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, plot_folder) plot_timeheight_p3( string("experiments/KiD_driver/", output_folder, "/Output.nc"), precip, From ea3a326e15a7793f475c353ece2c08691b9f7cdc Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Fri, 23 Aug 2024 15:48:08 -0700 Subject: [PATCH 14/22] oops again --- test/plotting_utils.jl | 129 +++++++++++++---------------------------- 1 file changed, 41 insertions(+), 88 deletions(-) diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index d414d8f6..9c078ae2 100644 --- a/test/plotting_utils.jl +++ b/test/plotting_utils.jl @@ -194,46 +194,24 @@ function plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, output for (i, u) in enumerate(solver.u) - if moisture isa CO.EquilibriumMoisture - q_tot[:, i] = parent(u.ρq_tot) ./ ρ .* 1e3 - elseif moisture isa CO.NonEquilibriumMoisture - q_tot[:, i] = parent(u.ρq_tot) ./ ρ .* 1e3 - q_liq[:, i] = parent(u.ρq_liq) ./ ρ .* 1e3 - q_ice[:, i] = parent(u.ρq_ice) ./ ρ .* 1e3 - elseif moisture isa CO.CloudyMoisture - q_tot[:, i] = (parent(u.moments.:2) .+ parent(u.ρq_vap)) ./ ρ .* 1e3 - q_liq[:, i] = parent(u.moments.:2) ./ ρ .* 1e3 - elseif precip isa CO.PrecipitationP3 - ρ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 - end - if precip isa CO.Precipitation1M - q_rai[:, i] = parent(u.ρq_rai) ./ ρ .* 1e3 - q_sno[:, i] = parent(u.ρq_sno) ./ ρ .* 1e3 - elseif precip isa CO.Precipitation2M - q_rai[:, i] = parent(u.ρq_rai) ./ ρ .* 1e3 - N_rai[:, i] = parent(u.N_rai) ./ 1e6 - N_liq[:, i] = parent(u.N_liq) ./ 1e6 - elseif precip isa CO.CloudyPrecip - q_rai[:, i] = parent(u.moments.:5) ./ ρ .* 1e3 - N_liq[:, i] = parent(u.moments.:1) ./ 1e6 - N_rai[:, i] = parent(u.moments.:4) ./ 1e6 - end + ρ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 + end function plot_data(data, data_label, max_val, title = "") return Plots.plot( data, - z_plt, + z_centers, xlabel = data_label, ylabel = "z [m]", legend = false, @@ -243,63 +221,38 @@ function plot_animation_p3(z_centers, solver, aux, moisture, precip, K1D, output ) end - anim = Plots.@animate for i in 1:length(t_plt) + anim = Plots.@animate for i in 1:nt title = "time = " * string(floor(Int, solver.t[i])) * " [s]" - if precip isa CO.PrecipitationP3 - 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), - ) - else - p1 = plot_data(q_tot[:, i], "q_tot [g/kg]", maximum(q_tot)) - p2 = plot_data(q_liq[:, i], "q_liq [g/kg]", maximum(q_liq), title) - p3 = plot_data(N_liq[:, i], "N_liq [1/cm^3]", maximum(N_liq)) - p4 = plot_data(q_rai[:, i], "q_rai [g/kg]", maximum(q_rai)) - p5 = plot_data(N_rai[:, i], "N_rai [1/cm^3]", maximum(N_rai)) - p6 = plot_data(q_ice[:, i], "q_ice [g/kg]", maximum(q_ice)) - p7 = plot_data(q_sno[:, i], "q_sno [g/kg]", maximum(q_sno)) - plot( - p1, - p2, - p3, - p4, - p5, - p6, - p7, - 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 = (3, 3), - ) - end + 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) From 0212e0aa11cf573118624205063424b0de208445 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Fri, 23 Aug 2024 16:46:55 -0700 Subject: [PATCH 15/22] codecov --- test/unit_tests/common_unit_test.jl | 9 ++++++ test/unit_tests/k1d_unit_test.jl | 47 +++++++++++++++++++++++++++++ 2 files changed, 56 insertions(+) diff --git a/test/unit_tests/common_unit_test.jl b/test/unit_tests/common_unit_test.jl index a2a3f413..bd47670b 100644 --- a/test/unit_tests/common_unit_test.jl +++ b/test/unit_tests/common_unit_test.jl @@ -325,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, equil_moist, ps) + 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 ea4d7c44..173674ed 100644 --- a/test/unit_tests/k1d_unit_test.jl +++ b/test/unit_tests/k1d_unit_test.jl @@ -161,6 +161,53 @@ end 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) From 1d363d74aaf55ce7f18a71b4b8be3a1e50c2a502 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Fri, 23 Aug 2024 16:49:38 -0700 Subject: [PATCH 16/22] frmt --- test/unit_tests/k1d_unit_test.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit_tests/k1d_unit_test.jl b/test/unit_tests/k1d_unit_test.jl index 173674ed..ada4db82 100644 --- a/test/unit_tests/k1d_unit_test.jl +++ b/test/unit_tests/k1d_unit_test.jl @@ -207,7 +207,7 @@ end 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) From 5c035b9dfd2a3755f9953ad9a3694695e9b70a2c Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Fri, 23 Aug 2024 16:59:27 -0700 Subject: [PATCH 17/22] test --- test/unit_tests/common_unit_test.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit_tests/common_unit_test.jl b/test/unit_tests/common_unit_test.jl index bd47670b..7e026e17 100644 --- a/test/unit_tests/common_unit_test.jl +++ b/test/unit_tests/common_unit_test.jl @@ -326,7 +326,7 @@ end end Y = CO.initialise_state(p3_moist, precip_p3, ip) - aux = CO.initialise_aux(FT, ip, params..., 0.0, 0.0, equil_moist, ps) + 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) From dbbef3a3d588d2fe4d5487d2f4fc3f9f27019a86 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Sat, 24 Aug 2024 16:31:39 -0700 Subject: [PATCH 18/22] wrong branch? --- Project.toml | 3 +- src/Common/ode_utils.jl | 14 +-- src/Common/tendency.jl | 114 +++++++++++++++++++++- test/experiments/KiD_driver/KiD_driver.jl | 22 ++++- 4 files changed, 140 insertions(+), 13 deletions(-) diff --git a/Project.toml b/Project.toml index f3212eba..fcaa1bf9 100644 --- a/Project.toml +++ b/Project.toml @@ -22,9 +22,11 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed" +TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] @@ -47,7 +49,6 @@ Plots = "1.29" ProgressLogging = "0.1" Random = "1" SpecialFunctions = "2.3" -Statistics = "1.10" TerminalLoggers = "0.1" Thermodynamics = "0.12.6" julia = "1.8" diff --git a/src/Common/ode_utils.jl b/src/Common/ode_utils.jl index 0f4b4b0a..0fc0bd16 100644 --- a/src/Common/ode_utils.jl +++ b/src/Common/ode_utils.jl @@ -141,6 +141,8 @@ function initialise_aux( tmp = similar(ip.q_tot), tmp2 = similar(ip.q_tot), tmp3 = similar(ip.q_tot), + tmp4 = similar(ip.q_tot), + tmp5 = similar(ip.q_tot), tmp_surface = similar(CC.Fields.level(ip.q_tot, 1), Tuple{FT, FT}), ) @@ -207,12 +209,12 @@ function initialise_aux( term_vel_N_ice = copy(ip.zero), ) precip_sources_eltype = @NamedTuple{ - q_tot::FT, - q_liq::FT, - q_rai::FT, - q_ice::FT, - q_rim::FT, - q_liqonice::FT, + ρq_tot::FT, + ρq_liq::FT, + ρq_rai::FT, + ρq_ice::FT, + ρq_rim::FT, + ρq_liqonice::FT, N_aer::FT, N_liq::FT, N_rai::FT, diff --git a/src/Common/tendency.jl b/src/Common/tendency.jl index 9807ab81..2339c247 100644 --- a/src/Common/tendency.jl +++ b/src/Common/tendency.jl @@ -191,13 +191,16 @@ end 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)) + @. ρ_r = ifelse(B_rim < eps(FT), FT(100), max(FT(100), ρ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) - + for x in [F_rim, F_liq] + @. x = ifelse(isnan(x), FT(0), x) + end + @. ρ_r = ifelse(isnan(ρ_r), FT(100), ρ_r) p3 = ps.p3_params Chen2022 = ps.Chen2022 @@ -559,8 +562,99 @@ end end @inline function precompute_aux_precip_sources!(ps::PrecipitationP3, aux) - # TODO [P3] - return nothing + (; dt) = aux.TS + (; thermo_params, air_params, common_params) = aux + FT = eltype(thermo_params) + (; ts, ρ, T) = aux.thermo_variables + (; ρq_tot, ρq_liq, ρq_rai, ρq_ice, ρq_rim, ρq_liqonice, N_ice, N_rai, B_rim) = aux.microph_variables + + dLdt = aux.scratch.tmp + dNdt = aux.scratch.tmp2 + _F_rim = aux.scratch.tmp3 + _ρ_rim = aux.scratch.tmp4 + _F_liq = aux.scratch.tmp5 + + p3 = ps.p3_params + vel = ps.Chen2022 + + # helper type and wrapper to populate the tuple with sources + precip_sources_eltype = @NamedTuple{ + ρq_tot::FT, + ρq_liq::FT, + ρq_rai::FT, + ρq_ice::FT, + ρq_rim::FT, + ρq_liqonice::FT, + N_aer::FT, + N_liq::FT, + N_rai::FT, + N_ice::FT, + B_rim::FT, + } + to_sources(args...) = @. precip_sources_eltype(tuple(args...)) + + # zero out the aux sources + @. aux.precip_sources = to_sources(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) + + # calculate temporary variables + @. _F_rim = max(FT(0), min(ρq_rim / (ρq_ice - ρq_liqonice), 1 - eps(FT))) + @. _ρ_rim = max(FT(100), min(ρq_rim / B_rim, p3.ρ_l)) + @. _F_liq = max(FT(0), min(ρq_liq / ρq_ice, FT(1))) + for x in [_F_rim, _F_liq] + @. x = ifelse(isnan(x), FT(0), x) + end + @. _ρ_rim = ifelse(isnan(_ρ_rim), FT(100), _ρ_rim) + if Bool(common_params.precip_sources) + # # Step 1: Broadcast the function call to obtain the results + # melt_results = CMP3.ice_melt.(p3, vel.snow_ice, air_params, thermo_params, ρq_ice, N_ice, T, ρ, _F_rim, _ρ_rim, dt) + + # # Step 2: Broadcast over the extracted fields and apply the limit function + # dLdt = limit.(ρq_ice, dt, getfield.(melt_results, :dLdt)) + # dNdt = limit.(N_ice, dt, getfield.(melt_results, :dNdt)) + + + @. dLdt = limit(ρq_ice, dt, CMP3.ice_melt(p3, vel.snow_ice, air_params, thermo_params, ρq_ice, N_ice, T, ρ, _F_rim, _ρ_rim, dt).dLdt) + @. dNdt = limit(N_ice, dt, CMP3.ice_melt(p3, vel.snow_ice, air_params, thermo_params, ρq_ice, N_ice, T, ρ, _F_rim, _ρ_rim, dt).dNdt) + # TODO - add in CM.jl what portion of melted ice + # is a source for ρq_liqonice and what portion goes to ρq_rai. + # TODO - also, calculate somehow what portion of the melted ice + # comes as a sink from ρq_rim and what portion as a sink + # from nonrimed ice. + # for now for testing purposes we assume: + # - source of ρq_liqonice = 0.7 * dLdt + # - source of ρq_rai = 0.3 * dLdt + # - sink of ρq_rim = _F_rim * dLdt + # - sink of ρq_ice = (1 - _F_rim) * (1 - _F_liq) * dLdt + # (this is because currently the notation in KiD + # is such that ρq_ice = ρq_liqonice + ρq_rim + # + other ice -- different than current P3 + # notation!!!) + # - sink of B_rim = _F_rim * dLdt / _ρ_rim + # which is good since it should be + # = [sink of L_rim] / _ρ_rim + # if we assume d[ρ_rim]/dt = 0 + # TODO - also, when we add the parameterization of melting + # as a source for ρq_liqonice, we need to keep in mind + # that the associated dNdt is zero (ice melts to liquid + # but the entire particle doesn't melt)... + # for now we assume dN_icedt = 0.3 * dNdt since that is how much + # we're assuming goes to rain at the moment + @. aux.precip_sources += to_sources( + 0, + 0, + 0.3 * dLdt, + -1 * (1 - _F_rim) * (1 - _F_liq) * dLdt, + -1 * _F_rim * dLdt, + 0.7 * dLdt, + 0, + 0, + 0.3 * dNdt, + -0.3 * dNdt, + _F_rim * dLdt / _ρ_rim, + ) + end + + end # helper function for precomputing aux sources for cloudy @@ -713,6 +807,18 @@ end return dY end @inline function precip_sources_tendency!(ms::MoistureP3, ps::PrecipitationP3, dY, Y, aux, t) + precompute_aux_precip_sources!(ps, aux) + + @. dY.ρq_tot += aux.precip_sources.ρq_tot + @. dY.ρq_liq += aux.precip_sources.ρq_liq + @. dY.ρq_rai += aux.precip_sources.ρq_rai + @. dY.ρq_ice += aux.precip_sources.ρq_ice + @. dY.ρq_rim += aux.precip_sources.ρq_rim + @. dY.ρq_liqonice += aux.precip_sources.ρq_liqonice + @. dY.B_rim += aux.precip_sources.B_rim + @. dY.N_ice += aux.precip_sources.N_ice + @. dY.N_rai += aux.precip_sources.N_rai + return dY end @inline function precip_sources_tendency!(ms::CloudyMoisture, ps::CloudyPrecip, dY, Y, aux, t) diff --git a/test/experiments/KiD_driver/KiD_driver.jl b/test/experiments/KiD_driver/KiD_driver.jl index e061dd02..2cd0b8c1 100644 --- a/test/experiments/KiD_driver/KiD_driver.jl +++ b/test/experiments/KiD_driver/KiD_driver.jl @@ -9,9 +9,27 @@ # Get the parameter values for the simulation include("parse_commandline.jl") -if !(@isdefined config) +bound = (; + ice_start = false, + _magnitude = Float64(0.5), + _q_flux = Float64(0.65e-4), + _N_flux = Float64(40000), + _F_rim = Float64(0.5), + _F_liq = Float64(0.0), + _ρ_r_init = Float64(900), + ) + +# if !(@isdefined config) config = parse_commandline() -end + config["precipitation_choice"] = "PrecipitationP3" + config["moisture_choice"] = "MoistureP3" + config["n_elem"] = 24 + config["z_max"] = 3000 + config["t_end"] = 75 + config["w1"] = 0 + config["dt"] = Float64(0.5) + config["p3_boundary_condition"] = bound +# end ft_choice = config["FLOAT_TYPE"] From 4aeada07d0f8253445e80308994cd0ac8907d3c5 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Sat, 24 Aug 2024 17:21:04 -0700 Subject: [PATCH 19/22] melting test code --- src/Common/tendency.jl | 71 +++++++++++++++++------ test/experiments/KiD_driver/KiD_driver.jl | 34 +++++------ 2 files changed, 70 insertions(+), 35 deletions(-) diff --git a/src/Common/tendency.jl b/src/Common/tendency.jl index 2339c247..97b87179 100644 --- a/src/Common/tendency.jl +++ b/src/Common/tendency.jl @@ -579,18 +579,18 @@ end # helper type and wrapper to populate the tuple with sources precip_sources_eltype = @NamedTuple{ - ρq_tot::FT, - ρq_liq::FT, - ρq_rai::FT, - ρq_ice::FT, - ρq_rim::FT, - ρq_liqonice::FT, - N_aer::FT, - N_liq::FT, - N_rai::FT, - N_ice::FT, - B_rim::FT, - } + ρq_tot::FT, + ρq_liq::FT, + ρq_rai::FT, + ρq_ice::FT, + ρq_rim::FT, + ρq_liqonice::FT, + N_aer::FT, + N_liq::FT, + N_rai::FT, + N_ice::FT, + B_rim::FT, + } to_sources(args...) = @. precip_sources_eltype(tuple(args...)) # zero out the aux sources @@ -613,8 +613,40 @@ end # dNdt = limit.(N_ice, dt, getfield.(melt_results, :dNdt)) - @. dLdt = limit(ρq_ice, dt, CMP3.ice_melt(p3, vel.snow_ice, air_params, thermo_params, ρq_ice, N_ice, T, ρ, _F_rim, _ρ_rim, dt).dLdt) - @. dNdt = limit(N_ice, dt, CMP3.ice_melt(p3, vel.snow_ice, air_params, thermo_params, ρq_ice, N_ice, T, ρ, _F_rim, _ρ_rim, dt).dNdt) + @. dLdt = limit( + ρq_ice, + dt, + CMP3.ice_melt( + p3, + vel.snow_ice, + air_params, + thermo_params, + (ρq_ice - ρq_liqonice), + N_ice, + T, + ρ, + _F_rim, + _ρ_rim, + dt, + ).dLdt, + ) + @. dNdt = limit( + N_ice, + dt, + CMP3.ice_melt( + p3, + vel.snow_ice, + air_params, + thermo_params, + (ρq_ice - ρq_liqonice), + N_ice, + T, + ρ, + _F_rim, + _ρ_rim, + dt, + ).dNdt, + ) # TODO - add in CM.jl what portion of melted ice # is a source for ρq_liqonice and what portion goes to ρq_rai. # TODO - also, calculate somehow what portion of the melted ice @@ -624,11 +656,14 @@ end # - source of ρq_liqonice = 0.7 * dLdt # - source of ρq_rai = 0.3 * dLdt # - sink of ρq_rim = _F_rim * dLdt - # - sink of ρq_ice = (1 - _F_rim) * (1 - _F_liq) * dLdt + # - change of ρq_ice = F_liq * (0.7 * dLdt) - (1 - F_liq) dLdt # (this is because currently the notation in KiD # is such that ρq_ice = ρq_liqonice + ρq_rim # + other ice -- different than current P3 - # notation!!!) + # notation!!! And so the change of ρq_ice should + # should be the source of ρq_liqonice minus the + # melted mass) + # # - sink of B_rim = _F_rim * dLdt / _ρ_rim # which is good since it should be # = [sink of L_rim] / _ρ_rim @@ -643,14 +678,14 @@ end 0, 0, 0.3 * dLdt, - -1 * (1 - _F_rim) * (1 - _F_liq) * dLdt, + _F_liq * (0.7 * dLdt) - (1 - _F_liq) * dLdt, -1 * _F_rim * dLdt, 0.7 * dLdt, 0, 0, 0.3 * dNdt, -0.3 * dNdt, - _F_rim * dLdt / _ρ_rim, + -1 * _F_rim * dLdt / _ρ_rim, ) end diff --git a/test/experiments/KiD_driver/KiD_driver.jl b/test/experiments/KiD_driver/KiD_driver.jl index 2cd0b8c1..1bdc49c4 100644 --- a/test/experiments/KiD_driver/KiD_driver.jl +++ b/test/experiments/KiD_driver/KiD_driver.jl @@ -10,25 +10,25 @@ include("parse_commandline.jl") bound = (; - ice_start = false, - _magnitude = Float64(0.5), - _q_flux = Float64(0.65e-4), - _N_flux = Float64(40000), - _F_rim = Float64(0.5), - _F_liq = Float64(0.0), - _ρ_r_init = Float64(900), - ) + ice_start = false, + _magnitude = Float64(0.5), + _q_flux = Float64(0.65e-4), + _N_flux = Float64(40000), + _F_rim = Float64(0.5), + _F_liq = Float64(0.0), + _ρ_r_init = Float64(900), +) # if !(@isdefined config) - config = parse_commandline() - config["precipitation_choice"] = "PrecipitationP3" - config["moisture_choice"] = "MoistureP3" - config["n_elem"] = 24 - config["z_max"] = 3000 - config["t_end"] = 75 - config["w1"] = 0 - config["dt"] = Float64(0.5) - config["p3_boundary_condition"] = bound +config = parse_commandline() +config["precipitation_choice"] = "PrecipitationP3" +config["moisture_choice"] = "MoistureP3" +config["n_elem"] = 24 +config["z_max"] = 3000 +config["t_end"] = 500 +config["w1"] = 0 +config["dt"] = Float64(0.5) +config["p3_boundary_condition"] = bound # end ft_choice = config["FLOAT_TYPE"] From f2f1f90ca3b7307bc25d1d3a82a5416295e0a8a7 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Sat, 24 Aug 2024 17:57:58 -0700 Subject: [PATCH 20/22] corrections --- src/Common/tendency.jl | 8 ++++---- test/experiments/KiD_driver/KiD_driver.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Common/tendency.jl b/src/Common/tendency.jl index 6762d6d8..c998ede3 100644 --- a/src/Common/tendency.jl +++ b/src/Common/tendency.jl @@ -657,13 +657,13 @@ end # - source of ρq_liqonice = 0.7 * dLdt # - source of ρq_rai = 0.3 * dLdt # - sink of ρq_rim = _F_rim * dLdt - # - change of ρq_ice = F_liq * (0.7 * dLdt) - (1 - F_liq) dLdt + # - change of ρq_ice = -0.3 * dLdt # (this is because currently the notation in KiD # is such that ρq_ice = ρq_liqonice + ρq_rim # + other ice -- different than current P3 # notation!!! And so the change of ρq_ice should - # should be the source of ρq_liqonice minus the - # melted mass) + # should be the mass that's lost to rain + # (solid ice that's melted into liquid is kept)) # # - sink of B_rim = _F_rim * dLdt / _ρ_rim # which is good since it should be @@ -679,7 +679,7 @@ end 0, 0, 0.3 * dLdt, - _F_liq * (0.7 * dLdt) - (1 - _F_liq) * dLdt, + -1 * (0.7 - _F_rim + (_F_rim + (0.3 - 0.7))) * dLdt, -1 * _F_rim * dLdt, 0.7 * dLdt, 0, diff --git a/test/experiments/KiD_driver/KiD_driver.jl b/test/experiments/KiD_driver/KiD_driver.jl index 1bdc49c4..68646985 100644 --- a/test/experiments/KiD_driver/KiD_driver.jl +++ b/test/experiments/KiD_driver/KiD_driver.jl @@ -14,7 +14,7 @@ bound = (; _magnitude = Float64(0.5), _q_flux = Float64(0.65e-4), _N_flux = Float64(40000), - _F_rim = Float64(0.5), + _F_rim = Float64(0.1), _F_liq = Float64(0.0), _ρ_r_init = Float64(900), ) From 8972cf97af305df789a292ee8a23df50b54ca9bf Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Wed, 28 Aug 2024 16:26:39 -0700 Subject: [PATCH 21/22] progress --- src/Common/initial_condition.jl | 4 + src/Common/ode_utils.jl | 14 +- src/Common/tendency.jl | 183 ++++++++++++---------- test/experiments/KiD_driver/KiD_driver.jl | 6 +- test/plotting_utils.jl | 4 +- 5 files changed, 124 insertions(+), 87 deletions(-) diff --git a/src/Common/initial_condition.jl b/src/Common/initial_condition.jl index ccbae713..b4240e3c 100644 --- a/src/Common/initial_condition.jl +++ b/src/Common/initial_condition.jl @@ -329,6 +329,10 @@ function p3_initial_condition( # use approximate temperature values from # sounding: kin1d/soundings/Tsfc2.txt # (path in p3 fortran github repo) + # although I think it could also be grabbed from + # this website: https://weather.uwyo.edu/upperair/sounding.html + # temperature is a constant gradient with a melting layer at 500 m, + # but of course that could easily be changed 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) diff --git a/src/Common/ode_utils.jl b/src/Common/ode_utils.jl index 0fc0bd16..54fbf0bc 100644 --- a/src/Common/ode_utils.jl +++ b/src/Common/ode_utils.jl @@ -202,6 +202,18 @@ function initialise_aux( q_vap = ip.q_vap, ρq_vap = ip.ρq_vap, ) + p3_predicted_eltype = @NamedTuple{ + F_rim::FT, + ρ_rim::FT, + F_liq::FT, + } + p3_predicted = @. p3_predicted_eltype( + tuple( + copy(ip.ρq_rim / (ip.ρq_ice - ip.ρq_liqonice)), + copy(ip.ρq_rim / (ip.ρq_ice - ip.ρq_liqonice) / ip.B_rim), + copy(ip.ρq_liqonice / ip.ρq_ice), + ), + ) velocities = (; term_vel_rai = copy(ip.zero), term_vel_ice = copy(ip.zero), @@ -299,7 +311,7 @@ function initialise_aux( if precip isa CloudyPrecip aux = merge(aux, (; cloudy_params, cloudy_variables)) elseif precip isa PrecipitationP3 - aux = merge(aux, (; p3_boundary_condition)) + aux = merge(aux, (; p3_boundary_condition, p3_predicted)) end return aux diff --git a/src/Common/tendency.jl b/src/Common/tendency.jl index c998ede3..02081f21 100644 --- a/src/Common/tendency.jl +++ b/src/Common/tendency.jl @@ -187,30 +187,26 @@ end @. 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), FT(100), max(FT(100), ρ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) - for x in [F_rim, F_liq] - @. x = ifelse(isnan(x), FT(0), x) - end - @. ρ_r = ifelse(isnan(ρ_r), FT(100), ρ_r) + + p3 = ps.p3_params + (; F_rim, ρ_rim, F_liq) = aux.p3_predicted + F_r(ρq_ice, ρq_rim, ρq_liqonice) = ifelse(ρq_ice > eps(FT), max(FT(0), min(ρq_rim / (ρq_ice - ρq_liqonice), 1 - eps(FT))), FT(0)) + ρ_r(ρq_rim, B_rim) = ifelse(B_rim > eps(FT), max(FT(0), min(ρq_rim / B_rim, p3.ρ_l)), FT(0)) + F_l(ρq_ice, ρq_liqonice) = ifelse(ρq_ice > eps(FT), max(FT(0), min(ρq_liqonice / ρq_ice, FT(1 - eps(FT)))), FT(0)) + @. F_rim = F_r(ρq_ice, ρq_rim, ρq_liqonice) + @. ρ_rim = ρ_r(ρq_rim, B_rim) + @. F_liq = F_l(ρq_ice, ρq_liqonice) 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 + # TODO - should we be using ρ or ρ_dry in the velocity calculation? @. 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) + getindex(CMP3.ice_terminal_velocity(p3, Chen2022, ρq_ice, N_ice, ρ_rim, F_rim, F_liq, ρ, 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) + getindex(CMP3.ice_terminal_velocity(p3, Chen2022, ρq_ice, N_ice, ρ_rim, F_rim, F_liq, ρ, use_aspect_ratio), 2) # adding 2M rain below: @@ -567,15 +563,17 @@ end (; thermo_params, air_params, common_params) = aux FT = eltype(thermo_params) (; ts, ρ, T) = aux.thermo_variables - (; ρq_tot, ρq_liq, ρq_rai, ρq_ice, ρq_rim, ρq_liqonice, N_ice, N_rai, B_rim) = aux.microph_variables + (; ρq_liq, ρq_ice, ρq_rim, ρq_liqonice, N_ice, B_rim) = aux.microph_variables + (; F_rim, ρ_rim, F_liq) = aux.p3_predicted - dLdt = aux.scratch.tmp - dNdt = aux.scratch.tmp2 - _F_rim = aux.scratch.tmp3 - _ρ_rim = aux.scratch.tmp4 - _F_liq = aux.scratch.tmp5 + dLdt_p3 = aux.scratch.tmp + dNdt_ice = aux.scratch.tmp2 + dLdt_rim = aux.scratch.tmp3 + dLdt_liq = aux.scratch.tmp4 + ddt_B_rim = aux.scratch.tmp5 p3 = ps.p3_params + type = typeof(p3) vel = ps.Chen2022 # helper type and wrapper to populate the tuple with sources @@ -598,95 +596,118 @@ end @. aux.precip_sources = to_sources(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) # calculate temporary variables - @. _F_rim = max(FT(0), min(ρq_rim / (ρq_ice - ρq_liqonice), 1 - eps(FT))) - @. _ρ_rim = max(FT(100), min(ρq_rim / B_rim, p3.ρ_l)) - @. _F_liq = max(FT(0), min(ρq_liq / ρq_ice, FT(1))) - for x in [_F_rim, _F_liq] - @. x = ifelse(isnan(x), FT(0), x) - end - @. _ρ_rim = ifelse(isnan(_ρ_rim), FT(100), _ρ_rim) - if Bool(common_params.precip_sources) - # # Step 1: Broadcast the function call to obtain the results - # melt_results = CMP3.ice_melt.(p3, vel.snow_ice, air_params, thermo_params, ρq_ice, N_ice, T, ρ, _F_rim, _ρ_rim, dt) - # # Step 2: Broadcast over the extracted fields and apply the limit function - # dLdt = limit.(ρq_ice, dt, getfield.(melt_results, :dLdt)) - # dNdt = limit.(N_ice, dt, getfield.(melt_results, :dNdt)) + F_r(ρq_ice, ρq_rim, ρq_liqonice) = ifelse(ρq_ice > eps(FT), max(FT(0), min(ρq_rim / (ρq_ice - ρq_liqonice), 1 - eps(FT))), FT(0)) + ρ_r(ρq_rim, B_rim) = ifelse(B_rim > 0, max(FT(0), min(ρq_rim / B_rim, p3.ρ_l)), FT(0)) + F_l(ρq_ice, ρq_liqonice) = ifelse(ρq_ice > eps(FT), max(FT(0), min(ρq_liqonice / ρq_ice, FT(1 - eps(FT)))), FT(0)) + @. F_rim = F_r(ρq_ice, ρq_rim, ρq_liqonice) + @. ρ_rim = ρ_r(ρq_rim, B_rim) + @. F_liq = F_l(ρq_ice, ρq_liqonice) + if Bool(common_params.precip_sources) - @. dLdt = limit( + @. dLdt_p3 = limit( + ρq_ice, + dt, + CMP3.ice_melt( + p3, + vel, + air_params, + thermo_params, + ρq_ice, + N_ice, + T, + ρ, + F_rim, + ρ_rim, + F_liq, + dt, + ).dLdt_p3_tot, + ) + @. dLdt_rim = limit( ρq_ice, dt, CMP3.ice_melt( p3, - vel.snow_ice, + vel, + air_params, + thermo_params, + ρq_ice, + N_ice, + T, + ρ, + F_rim, + ρ_rim, + F_liq, + dt, + ).dLdt_rim, + ) + @. dLdt_liq = limit( + ρq_ice, + dt, + CMP3.ice_melt( + p3, + vel, + air_params, + thermo_params, + ρq_ice, + N_ice, + T, + ρ, + F_rim, + ρ_rim, + F_liq, + dt, + ).dLdt_liq, + ) + @. dNdt_ice = limit( + N_ice, + dt, + CMP3.ice_melt( + p3, + vel, air_params, thermo_params, - (ρq_ice - ρq_liqonice), + ρq_ice, N_ice, T, ρ, - _F_rim, - _ρ_rim, + F_rim, + ρ_rim, + F_liq, dt, - ).dLdt, + ).dNdt_ice, ) - @. dNdt = limit( + @. ddt_B_rim= limit( N_ice, dt, CMP3.ice_melt( p3, - vel.snow_ice, + vel, air_params, thermo_params, - (ρq_ice - ρq_liqonice), + ρq_ice, N_ice, T, ρ, - _F_rim, - _ρ_rim, + F_rim, + ρ_rim, + F_liq, dt, - ).dNdt, + ).ddtB_rim, ) - # TODO - add in CM.jl what portion of melted ice - # is a source for ρq_liqonice and what portion goes to ρq_rai. - # TODO - also, calculate somehow what portion of the melted ice - # comes as a sink from ρq_rim and what portion as a sink - # from nonrimed ice. - # for now for testing purposes we assume: - # - source of ρq_liqonice = 0.7 * dLdt - # - source of ρq_rai = 0.3 * dLdt - # - sink of ρq_rim = _F_rim * dLdt - # - change of ρq_ice = -0.3 * dLdt - # (this is because currently the notation in KiD - # is such that ρq_ice = ρq_liqonice + ρq_rim - # + other ice -- different than current P3 - # notation!!! And so the change of ρq_ice should - # should be the mass that's lost to rain - # (solid ice that's melted into liquid is kept)) - # - # - sink of B_rim = _F_rim * dLdt / _ρ_rim - # which is good since it should be - # = [sink of L_rim] / _ρ_rim - # if we assume d[ρ_rim]/dt = 0 - # TODO - also, when we add the parameterization of melting - # as a source for ρq_liqonice, we need to keep in mind - # that the associated dNdt is zero (ice melts to liquid - # but the entire particle doesn't melt)... - # for now we assume dN_icedt = 0.3 * dNdt since that is how much - # we're assuming goes to rain at the moment @. aux.precip_sources += to_sources( 0, 0, - 0.3 * dLdt, - -1 * (0.7 - _F_rim + (_F_rim + (0.3 - 0.7))) * dLdt, - -1 * _F_rim * dLdt, - 0.7 * dLdt, + dLdt_p3, # source of rain = sink of p3 ice + -dLdt_p3, # sink of p3 ice + -dLdt_rim, # sink of rime + dLdt_liq, # source of L_liq 0, 0, - 0.3 * dNdt, - -0.3 * dNdt, - -1 * _F_rim * dLdt / _ρ_rim, + dNdt_ice, # source of rain number + -dNdt_ice, # sink of ice number + -ddt_B_rim, # sink of B_rim ) end end diff --git a/test/experiments/KiD_driver/KiD_driver.jl b/test/experiments/KiD_driver/KiD_driver.jl index 68646985..f729c064 100644 --- a/test/experiments/KiD_driver/KiD_driver.jl +++ b/test/experiments/KiD_driver/KiD_driver.jl @@ -13,7 +13,7 @@ bound = (; ice_start = false, _magnitude = Float64(0.5), _q_flux = Float64(0.65e-4), - _N_flux = Float64(40000), + _N_flux = Float64(1e5), _F_rim = Float64(0.1), _F_liq = Float64(0.0), _ρ_r_init = Float64(900), @@ -23,9 +23,9 @@ bound = (; config = parse_commandline() config["precipitation_choice"] = "PrecipitationP3" config["moisture_choice"] = "MoistureP3" -config["n_elem"] = 24 +config["n_elem"] = 36 config["z_max"] = 3000 -config["t_end"] = 500 +config["t_end"] = 75 config["w1"] = 0 config["dt"] = Float64(0.5) config["p3_boundary_condition"] = bound diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index 9c078ae2..d1f9a203 100644 --- a/test/plotting_utils.jl +++ b/test/plotting_utils.jl @@ -97,7 +97,7 @@ function plot_final_aux_profiles(z_centers, aux, precip; output = "output") 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]") + p11 = Plots.plot(N_rai_end .* 1e-6, z_centers, xlabel = "N_rai [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]") @@ -340,7 +340,7 @@ function plot_timeheight_p3(nc_data_file, precip; output = "output") 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) + p2 = Plots.heatmap(t_plt, z_plt, (ρq_ice_plt .- ρq_liqonice_plt) .* 1e3, title = "ρq_ice (solid ice) [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) From 2c6ca60868cdd0c66923991905e81e27d08279e0 Mon Sep 17 00:00:00 2001 From: rorlija1 Date: Fri, 30 Aug 2024 15:40:58 -0700 Subject: [PATCH 22/22] updates --- src/Common/tendency.jl | 124 ++++++++++++++++++++-- test/experiments/KiD_driver/KiD_driver.jl | 6 +- 2 files changed, 121 insertions(+), 9 deletions(-) diff --git a/src/Common/tendency.jl b/src/Common/tendency.jl index 02081f21..6ed1859b 100644 --- a/src/Common/tendency.jl +++ b/src/Common/tendency.jl @@ -563,11 +563,11 @@ end (; thermo_params, air_params, common_params) = aux FT = eltype(thermo_params) (; ts, ρ, T) = aux.thermo_variables - (; ρq_liq, ρq_ice, ρq_rim, ρq_liqonice, N_ice, B_rim) = aux.microph_variables + (; ρq_liq, ρq_ice, ρq_rim, ρq_liqonice, N_ice, B_rim, ρq_rai, N_rai) = aux.microph_variables (; F_rim, ρ_rim, F_liq) = aux.p3_predicted dLdt_p3 = aux.scratch.tmp - dNdt_ice = aux.scratch.tmp2 + dNdt = aux.scratch.tmp2 dLdt_rim = aux.scratch.tmp3 dLdt_liq = aux.scratch.tmp4 ddt_B_rim = aux.scratch.tmp5 @@ -575,6 +575,8 @@ end p3 = ps.p3_params type = typeof(p3) vel = ps.Chen2022 + sb2006 = ps.sb2006 + aps = CMP.AirProperties(FT) # helper type and wrapper to populate the tuple with sources precip_sources_eltype = @NamedTuple{ @@ -604,7 +606,10 @@ end @. ρ_rim = ρ_r(ρq_rim, B_rim) @. F_liq = F_l(ρq_ice, ρq_liqonice) - if Bool(common_params.precip_sources) + # melting: turned off for testing shedding + # if Bool(common_params.precip_sources) + melting = false + if melting # melting works with corresponding CM branch! @. dLdt_p3 = limit( ρq_ice, @@ -660,7 +665,7 @@ end dt, ).dLdt_liq, ) - @. dNdt_ice = limit( + @. dNdt = limit( N_ice, dt, CMP3.ice_melt( @@ -705,11 +710,116 @@ end dLdt_liq, # source of L_liq 0, 0, - dNdt_ice, # source of rain number - -dNdt_ice, # sink of ice number + dNdt, # source of rain number + -dNdt, # sink of ice number -ddt_B_rim, # sink of B_rim ) end + shedding = false + if shedding # shedding works with corresponding CM branch! + @. dLdt_p3 = limit( + ρq_ice, + dt, + CMP3.ice_shed( + p3, + ρq_ice, + N_ice, + F_rim, + ρ_rim, + F_liq, + dt, + ).dLdt_p3_tot, + ) + @. dNdt = limit( + N_ice, + dt, + CMP3.ice_shed( + p3, + ρq_ice, + N_ice, + F_rim, + ρ_rim, + F_liq, + dt, + ).dNdt_rai, + ) + @. aux.precip_sources += to_sources( + 0, + 0, + dLdt_p3, # source of rain = sink of p3 ice + -dLdt_p3, # sink of p3 ice + 0, # sink of rime + -dLdt_p3, # sink of L_liq + 0, + 0, + dNdt, # source of rain number + 0, # sink of ice number + 0, # sink of B_rim + ) + end + collisions = true + if collisions # with rain - stable for corresponding CM branch (collision kernel) + # but haven't really been able to test much at all + # TODO - use ρ or ρ_dry? + @. dLdt_p3 = limit( + ρq_rai, + dt, + CMP3.ice_collisions( + sb2006.pdf_r, + p3, + vel, + aps, + thermo_params, + ts, + ρq_ice, + N_ice, + ρq_rai, + N_rai, + ρ, + F_rim, + ρ_rim, + F_liq, + T, + dt, + ).dLdt_p3_tot, + ) + @. dNdt = limit( + N_rai, + dt, + CMP3.ice_collisions( + sb2006.pdf_r, + p3, + vel, + aps, + thermo_params, + ts, + ρq_ice, + N_ice, + ρq_rai, + N_rai, + ρ, + F_rim, + ρ_rim, + F_liq, + T, + dt, + ).dNdt_rai, + ) + @. dLdt_p3 = ifelse(isnan(dLdt_p3), FT(0), dLdt_p3) + @. aux.precip_sources += to_sources( + 0, + 0, + -dLdt_p3, # sink of rain = source of p3 ice + dLdt_p3, # source of p3 ice + dLdt_p3, # source of rime + dLdt_p3, # source of L_liq + 0, + 0, + -dNdt, # sink of rain number + 0, # sink of ice number + dLdt_p3 / FT(900), # source of B_rim + ) + end end # helper function for precomputing aux sources for cloudy @@ -873,6 +983,8 @@ end @. dY.B_rim += aux.precip_sources.B_rim @. dY.N_ice += aux.precip_sources.N_ice @. dY.N_rai += aux.precip_sources.N_rai + @. dY.N_aer += aux.precip_sources.N_aer + @. dY.N_liq += aux.precip_sources.N_liq return dY end diff --git a/test/experiments/KiD_driver/KiD_driver.jl b/test/experiments/KiD_driver/KiD_driver.jl index f729c064..2716448f 100644 --- a/test/experiments/KiD_driver/KiD_driver.jl +++ b/test/experiments/KiD_driver/KiD_driver.jl @@ -12,10 +12,10 @@ include("parse_commandline.jl") bound = (; ice_start = false, _magnitude = Float64(0.5), - _q_flux = Float64(0.65e-4), + _q_flux = Float64(0.65e-3), _N_flux = Float64(1e5), - _F_rim = Float64(0.1), - _F_liq = Float64(0.0), + _F_rim = Float64(0.8), + _F_liq = Float64(0.5), _ρ_r_init = Float64(900), )