From 957694a410fa00350858b15916bebadd40aa44e4 Mon Sep 17 00:00:00 2001 From: Eric Pham Date: Tue, 12 Aug 2025 15:20:15 -0700 Subject: [PATCH] add modifying config variables, add docs, add diagnostic EDMF yaml --- .buildkite/pipeline.yml | 9 ++ config/default_configs/default_config.yml | 13 ++ .../diagnostic_edmfx_dycoms_rf02_column.yml | 41 ++++++ ...tic_edmfx_perturbed_dycoms_rf02_column.yml | 47 +++++++ docs/make.jl | 5 +- docs/src/contributor_guide.md | 11 +- docs/src/perturbation_experiments.md | 112 ++++++++++++++++ post_processing/ci_plots.jl | 1 + src/initial_conditions/initial_conditions.jl | 124 ++++++++++++++---- src/solver/type_getters.jl | 9 +- 10 files changed, 341 insertions(+), 31 deletions(-) create mode 100644 config/model_configs/diagnostic_edmfx_dycoms_rf02_column.yml create mode 100644 config/model_configs/prognostic_edmfx_perturbed_dycoms_rf02_column.yml create mode 100644 docs/src/perturbation_experiments.md diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 0dd2225a23..dc55ce19e5 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -844,6 +844,15 @@ steps: artifact_paths: "prognostic_edmfx_dycoms_rf02_column/output_active/*" agents: slurm_mem: 20GB + + - label: ":umbrella: Prognostic EDMFX Perturbed Dycoms RF02 in a column" + command: > + julia --color=yes --project=.buildkite .buildkite/ci_driver.jl + --config_file $CONFIG_PATH/prognostic_edmfx_perturbed_dycoms_rf02_column.yml + --job_id prognostic_edmfx_perturbed_dycoms_rf02_column.yml + artifact_paths: "prognostic_edmfx_perturbed_dycoms_rf02_column/output_active/*" + agents: + slurm_mem: 20GB - label: ":umbrella: Prognostic EDMFX Rico in a column" command: > diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml index d49e00e50d..3a5f26b9c2 100644 --- a/config/default_configs/default_config.yml +++ b/config/default_configs/default_config.yml @@ -398,3 +398,16 @@ use_itime: 1. Surface conditions that explicitly depend on time (e.g. LifeCycleTan2018, TRMM_LBA, etc.), 2. Time dependent forcing/tendencies use time rounded to the nearest unit of time for dt" value: false +# DYCOMS-RF02 Perturbation Variables +q_tot_0_dycoms_rf02: + help: "Surface total water mixing ratio for DYCOMS RF02." + value: 9.45 +theta_0_dycoms_rf02: + help: "Surface liquid-ice potential temperature for DYCOMS RF02." + value: 288.3 +theta_i_dycoms_rf02: + help: "Initial boundary layer liquid-ice potential temperature for DYCOMS RF02." + value: 295.0 +z_i_dycoms_rf02: + help: "Initial boundary layer height." + value: 795.0 \ No newline at end of file diff --git a/config/model_configs/diagnostic_edmfx_dycoms_rf02_column.yml b/config/model_configs/diagnostic_edmfx_dycoms_rf02_column.yml new file mode 100644 index 0000000000..d95ff50885 --- /dev/null +++ b/config/model_configs/diagnostic_edmfx_dycoms_rf02_column.yml @@ -0,0 +1,41 @@ +initial_condition: DYCOMS_RF02 +subsidence: DYCOMS +scm_coriolis: DYCOMS_RF02 +rad: DYCOMS +surface_setup: DYCOMS_RF02 +turbconv: "diagnostic_edmfx" +implicit_diffusion: true +approximate_linear_solve_iters: 2 +edmfx_upwinding: first_order +edmfx_entr_model: "Generalized" +edmfx_detr_model: "Generalized" +edmfx_sgs_mass_flux: true +edmfx_sgs_diffusive_flux: true +edmfx_nh_pressure: true +edmfx_filter: true +prognostic_tke: true +moist: "nonequil" +cloud_model: "grid_scale" +precip_model: "1M" +call_cloud_diagnostics_per_stage: true +config: "column" +x_elem: 2 +y_elem: 2 +z_elem: 30 +z_max: 1500 +z_stretch: false +perturb_initstate: false +dt: 10secs +t_end: 5hours +dt_save_state_to_disk: 10mins +toml: [toml/diagnostic_edmfx_1M.toml] +netcdf_interpolation_num_points: [8, 8, 30] +diagnostics: + - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr, lwp, rwp] + period: 10mins + - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, husraup, hussnup, tke] + period: 10mins + - short_name: [entr, detr, lmix, bgrad, strain, edt, evu] + period: 10mins + - short_name: [husra, hussn] + period: 10mins diff --git a/config/model_configs/prognostic_edmfx_perturbed_dycoms_rf02_column.yml b/config/model_configs/prognostic_edmfx_perturbed_dycoms_rf02_column.yml new file mode 100644 index 0000000000..6a247e8860 --- /dev/null +++ b/config/model_configs/prognostic_edmfx_perturbed_dycoms_rf02_column.yml @@ -0,0 +1,47 @@ +initial_condition: DYCOMS_RF02 +subsidence: DYCOMS +scm_coriolis: DYCOMS_RF02 +rad: DYCOMS +surface_setup: DYCOMS_RF02 +q_tot_0_dycoms_rf02: 8.0 +theta_0_dycoms_rf02: 284.0 +theta_i_dycoms_rf02: 291.0 +z_max: 800.0 +turbconv: "prognostic_edmfx" +implicit_diffusion: false +approximate_linear_solve_iters: 2 +edmfx_upwinding: first_order +edmfx_entr_model: "Generalized" +edmfx_detr_model: "Generalized" +edmfx_sgs_mass_flux: true +edmfx_sgs_diffusive_flux: true +edmfx_nh_pressure: true +edmfx_filter: true +prognostic_tke: true +moist: "nonequil" +cloud_model: "quadrature_sgs" +precip_model: "1M" +call_cloud_diagnostics_per_stage: true +config: "column" +x_elem: 2 +y_elem: 2 +z_elem: 30 +z_max: 1500 +z_stretch: false +perturb_initstate: false +dt: 10secs +t_end: 5hours +dt_save_state_to_disk: 10mins +toml: [toml/prognostic_edmfx_1M.toml] +netcdf_interpolation_num_points: [8, 8, 30] +diagnostics: + - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr] + period: 10mins + - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, husraup, hussnup] + period: 10mins + - short_name: [waen, taen, thetaaen, haen, husen, huren, clwen, clien, husraen, hussnen, tke] + period: 10mins + - short_name: [entr, detr, lmix, bgrad, strain, edt, evu] + period: 10mins + - short_name: [husra, hussn] + period: 10mins diff --git a/docs/make.jl b/docs/make.jl index 8938d9919c..9bdd8f1e3c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,6 +10,8 @@ include(joinpath(@__DIR__, "src", "config_table.jl")) doctest(ClimaAtmos; plugins = [bib]) disable_logging(Base.CoreLogging.BelowMinLevel) # Re-enable all logging +example_pages = ["Perturbation Experiments" => "perturbation_experiments.md"] + makedocs(; plugins = [bib], modules = [ClimaAtmos], @@ -38,8 +40,9 @@ makedocs(; "Implicit Solver" => "implicit_solver.md", "Radiative Equilibrium" => "radiative_equilibrium.md", "Single Column Model" => "single_column_prospect.md", - "Restarts and checkpoints" => "restarts.md", + "Restarts and Checkpoints" => "restarts.md", "REPL scripts" => "repl_scripts.md", + "Examples" => example_pages, "Configuration" => "config.md", "Parameters" => "parameters.md", "Longruns" => "longruns.md", diff --git a/docs/src/contributor_guide.md b/docs/src/contributor_guide.md index 7038b915bc..cfd8823d9d 100644 --- a/docs/src/contributor_guide.md +++ b/docs/src/contributor_guide.md @@ -165,7 +165,16 @@ or clarification that you can take care of while you are here. Here is an example of a docstring: -TODO: add example +``` +""" + InitialCondition + +A mechanism for specifying the `LocalState` of an `AtmosModel` at every point in +the domain. Given some `initial_condition`, calling `initial_condition(params)` +returns a function of the form `local_state(local_geometry)::LocalState`. +""" +abstract type InitialCondition end +``` You can preview how the Documentation will look like after merging by building the documentation locally. From the main directory of your local repository call diff --git a/docs/src/perturbation_experiments.md b/docs/src/perturbation_experiments.md new file mode 100644 index 0000000000..b7a2b49b5d --- /dev/null +++ b/docs/src/perturbation_experiments.md @@ -0,0 +1,112 @@ + +# Creating custom arguments for perturbation experiments +In cases where it is of interest to run ensembles of the same model configuration but with perturbed initial conditions, it may be of use to add new keyword arguments to the .yml file. The idea is to allow modifications to initial conditions to be passed directly through the .yml file and then creating multiple copies of the files but with different variations and combinations of these parameters. + +As an example, we will explore modifying the total water mixing ratio (`q_tot_0`), liquid-ice potential temperature profiles (`theta_0` & `theta_i`), and initial boundary layer height (`z_i`) in the DYCOMS-RF02 simulation setup. + +### Modify initial_conditions.jl +To start, we need to go to the part of the source code responsible for defining the initial conditions. For both of the DYCOMS research flight simulation setups, the total water mixing ratio and liquid-ice potential temperature are pulled from another library called AtmosphericProfilesLibrary (APL). + +There are many ways to modify the functions to allow for perturbation to the initial conditions, but for the sake of this example we will overwrite the function signature with our own version where we can pass in our own custom profiles. We must import the functions from the APL library and then define a new function with the same signature. This looks like: + +``` +import AtmosphericProfilesLibrary: Dycoms_RF02_θ_liq_ice, Dycoms_RF02_q_tot +[...] + +# Redefine the profile functions here. Here we redefine the functions such that +# we can pass in our own values for q_tot and θ_liq_ice, as well as modifying +# the initial boundary layer height. + +""" [Ackerman2009](@cite) """ +Dycoms_RF02_θ_liq_ice(::Type{FT}, theta_0, theta_i, z_i) where {FT} = + APL.ZProfile(z -> if z <= z_i + FT(theta_0) + else + FT(theta_i) + (z - FT(z_i))^FT(1.0 / 3.0) + end) + +""" [Ackerman2009](@cite) """ +Dycoms_RF02_q_tot(::Type{FT}, q_tot_0, z_i) where {FT} = + APL.ZProfile(z -> if z <= z_i + FT(q_tot_0) / FT(1000.0) + else + (FT(5) - FT(3) * (FT(1) - exp(-(z - FT(z_i)) / FT(500)))) / FT(1000) + end) +``` + +Now that we have redefined our functions, we can pass these new functions to our intial conditions structure for the DYCOMS-RF02 setup. Here we can also begin to define our new keyword arguments that will serve as inputs for the new functions we defined. In the `DYCOMS_RF02` structure, we add the names of our new inputs: + +``` +Base.@kwdef struct DYCOMS_RF02 <: InitialCondition + prognostic_tke::Bool = false + q_tot_0_dycoms_rf02 # Define total water mixing ratio. + theta_0_dycoms_rf02 # Define liquid-ice θ at surface. + theta_i_dycoms_rf02 # Define liquid-ice θ at initial boundary layer height. + z_i_dycoms_rf02 # Define initial boundary layer height. +end + +for IC in (:Dycoms_RF01, :Dycoms_RF02) + IC_Type = Symbol(uppercase(string(IC))) + θ_func_name = Symbol(IC, :_θ_liq_ice) + q_tot_func_name = Symbol(IC, :_q_tot) + [...] + if IC == :Dycoms_RF02 + @eval function (initial_condition::$IC_Type)(params) + (; prognostic_tke, q_tot_0_dycoms_rf02, theta_0_dycoms_rf02, theta_i_dycoms_rf02, z_i_dycoms_rf02) = initial_condition #unpack the new arguments here. These arguments will be provided through the model .yml file. + FT = eltype(params) + [...] + θ = $θ_func_name(FT, FT(theta_0_dycoms_rf02), FT(theta_i_dycoms_rf02), FT(z_i_dycoms_rf02)) # Change function signature here. + q_tot = $q_tot_func_name(FT, FT(q_tot_0_dycoms_rf02), FT(z_i_dycoms_rf02)) # Change function signature here. + end + else + [...] + end +end +``` + +### Add keywords to .yml & Modify default_config.yml +Now that we have added a new keyword argument to be used in initial conditions, it is time to define the keyword argument in the YAML files responsible for configuring the model runs. In the YAML file for running a DYCOMS-RF02 single column experiment, we add and adjust the following lines: + +``` +q_tot_0_dycoms_rf02: 9.45 # Define variables here. +theta_0_dycoms_rf02: 288.3 # Define variables here. +theta_i_dycoms_rf02: 295.0 # Define variables here. +z_i_dycoms_rf02: 795.0 # Define variables here. +``` + +These are the same default values being used by the original APL function, but we can modify it in the YAML file if we want to test different initial conditions. Additionally, we need to provide a backup default value in the case that the keyword argument is not used or provided in the setup YAML file. To do this, we go to default_config.yml and add the following: + +``` +q_tot_0_dycoms_rf02: # Additional variables here. + help: "Surface total water mixing ratio for DYCOMS RF02." + value: 9.45 # Default value. +theta_0_dycoms_rf02: # Additional variables here. + help: "Surface liquid-ice potential temperature for DYCOMS RF02." + value: 288.3 # Default value. +theta_i_dycoms_rf02: # Additional variables here. + help: "Initial boundary layer liquid-ice potential temperature for DYCOMS RF02." + value: 295.0 # Default value. +z_i_dycoms_rf02: # Additional variables here. + help: "Initial boundary layer height." + value: 795.0 # Default value. +``` + +### Modify type_getters.jl +Finally, we need to update type_getters.jl with our new keyword arguments as well. We find where DYCOMS-RF02 is specified and then we replace it with a separate block that looks like: + +``` +[...] + elseif parsed_args["initial_condition"] == "DYCOMS_RF02" + return getproperty(ICs, Symbol(parsed_args["initial_condition"]))( + parsed_args["prognostic_tke"], + parsed_args["q_tot_0_dycoms_rf02"], # Add parsed args here. + parsed_args["theta_0_dycoms_rf02"], # Add parsed args here. + parsed_args["theta_i_dycoms_rf02"], # Add parsed args here. + parsed_args["z_i_dycoms_rf02"], # Add parsed args here. + ) +[...] +``` + +With this step complete, we are now ready to pass a new keyword argument that we defined ourselves to modify or change the initial conditions. + +The recommended workflow is the create multiple copies of the .yml files and label them according to the initial conditions used. For example, the default file might look like `prognostic_edmfx_dycoms_rf02_column_qtot0_6.5_theta0_284.0_thetai_290.0_zi_795.0_prescribedN_1.0e8.yml`. diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl index 3ab559139b..0615b1e1f6 100644 --- a/post_processing/ci_plots.jl +++ b/post_processing/ci_plots.jl @@ -1409,6 +1409,7 @@ DiagEDMFBoxPlotsWithPrecip = Union{ Val{:diagnostic_edmfx_dycoms_rf02_box}, Val{:diagnostic_edmfx_rico_box}, Val{:diagnostic_edmfx_trmm_box}, + Val{:diagnostic_edmfx_dycoms_rf02_column}, } """ plot_edmf_vert_profile!(grid_loc, var_group) diff --git a/src/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl index fb4f4138a7..e8200bad65 100644 --- a/src/initial_conditions/initial_conditions.jl +++ b/src/initial_conditions/initial_conditions.jl @@ -881,6 +881,7 @@ end ## # TODO: Get rid of this import AtmosphericProfilesLibrary as APL +import AtmosphericProfilesLibrary: Dycoms_RF02_θ_liq_ice, Dycoms_RF02_q_tot const FunctionOrSpline = Union{Function, APL.AbstractProfile, Intp.Extrapolation} @@ -1107,6 +1108,7 @@ hydrostatically balanced pressure profile. """ Base.@kwdef struct DYCOMS_RF01 <: InitialCondition prognostic_tke::Bool = false + end """ @@ -1117,8 +1119,29 @@ hydrostatically balanced pressure profile. """ Base.@kwdef struct DYCOMS_RF02 <: InitialCondition prognostic_tke::Bool = false + q_tot_0_dycoms_rf02::Any + theta_0_dycoms_rf02::Any + theta_i_dycoms_rf02::Any + z_i_dycoms_rf02::Any end +""" [Ackerman2009](@cite) """ +Dycoms_RF02_θ_liq_ice(::Type{FT}, theta_0, theta_i, z_i) where {FT} = + APL.ZProfile(z -> if z <= z_i + FT(theta_0) + else + FT(theta_i) + (z - FT(z_i))^FT(1.0 / 3.0) + end) + +""" [Ackerman2009](@cite) """ +Dycoms_RF02_q_tot(::Type{FT}, q_tot_0, z_i) where {FT} = APL.ZProfile( + z -> if z <= z_i + FT(q_tot_0) / FT(1000.0) + else + (FT(5) - FT(3) * (FT(1) - exp(-(z - FT(z_i)) / FT(500)))) / FT(1000) + end, +) + for IC in (:Dycoms_RF01, :Dycoms_RF02) IC_Type = Symbol(uppercase(string(IC))) θ_func_name = Symbol(IC, :_θ_liq_ice) @@ -1126,36 +1149,81 @@ for IC in (:Dycoms_RF01, :Dycoms_RF02) u_func_name = Symbol(IC, IC == :Dycoms_RF01 ? :_u0 : :_u) v_func_name = Symbol(IC, IC == :Dycoms_RF01 ? :_v0 : :_v) tke_func_name = Symbol(IC, :_tke_prescribed) - @eval function (initial_condition::$IC_Type)(params) - (; prognostic_tke) = initial_condition - FT = eltype(params) - thermo_params = CAP.thermodynamics_params(params) - p_0 = FT(101780.0) - θ = APL.$θ_func_name(FT) - q_tot = APL.$q_tot_func_name(FT) - p = hydrostatic_pressure_profile(; thermo_params, p_0, θ, q_tot) - u = APL.$u_func_name(FT) - v = APL.$v_func_name(FT) - #tke = APL.$tke_func_name(FT) - tke = APL.Dycoms_RF01_tke_prescribed(FT) #TODO - dont have the tke profile for Dycoms_RF02 - function local_state(local_geometry) - (; z) = local_geometry.coordinates - return LocalState(; - params, - geometry = local_geometry, - thermo_state = TD.PhaseEquil_pθq( - thermo_params, - p(z), - θ(z), - q_tot(z), - ), - velocity = Geometry.UVVector(u(z), v(z)), - turbconv_state = EDMFState(; - tke = prognostic_tke ? FT(0) : tke(z), - ), + if IC == :Dycoms_RF02 + @eval function (initial_condition::$IC_Type)(params) + (; + prognostic_tke, + q_tot_0_dycoms_rf02, + theta_0_dycoms_rf02, + theta_i_dycoms_rf02, + z_i_dycoms_rf02, + ) = initial_condition + FT = eltype(params) + thermo_params = CAP.thermodynamics_params(params) + p_0 = FT(101780.0) + θ = $θ_func_name( + FT, + theta_0_dycoms_rf02, + theta_i_dycoms_rf02, + z_i_dycoms_rf02, ) + q_tot = $q_tot_func_name(FT, q_tot_0_dycoms_rf02, z_i_dycoms_rf02) + p = hydrostatic_pressure_profile(; thermo_params, p_0, θ, q_tot) + u = APL.$u_func_name(FT) + v = APL.$v_func_name(FT) + #tke = APL.$tke_func_name(FT) + tke = APL.Dycoms_RF01_tke_prescribed(FT) #TODO - dont have the tke profile for Dycoms_RF02 + function local_state(local_geometry) + (; z) = local_geometry.coordinates + return LocalState(; + params, + geometry = local_geometry, + thermo_state = TD.PhaseEquil_pθq( + thermo_params, + p(z), + θ(z), + q_tot(z), + ), + velocity = Geometry.UVVector(u(z), v(z)), + turbconv_state = EDMFState(; + tke = prognostic_tke ? FT(0) : tke(z), + ), + ) + end + return local_state + end + else + @eval function (initial_condition::$IC_Type)(params) + (; prognostic_tke) = initial_condition + FT = eltype(params) + thermo_params = CAP.thermodynamics_params(params) + p_0 = FT(101780.0) + θ = APL.$θ_func_name(FT) + q_tot = APL.$q_tot_func_name(FT) + p = hydrostatic_pressure_profile(; thermo_params, p_0, θ, q_tot) + u = APL.$u_func_name(FT) + v = APL.$v_func_name(FT) + #tke = APL.$tke_func_name(FT) + tke = APL.Dycoms_RF01_tke_prescribed(FT) #TODO - dont have the tke profile for Dycoms_RF02 + function local_state(local_geometry) + (; z) = local_geometry.coordinates + return LocalState(; + params, + geometry = local_geometry, + thermo_state = TD.PhaseEquil_pθq( + thermo_params, + p(z), + θ(z), + q_tot(z), + ), + velocity = Geometry.UVVector(u(z), v(z)), + turbconv_state = EDMFState(; + tke = prognostic_tke ? FT(0) : tke(z), + ), + ) + end + return local_state end - return local_state end end diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index bbea414285..82311c5be5 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -381,7 +381,6 @@ function get_initial_condition(parsed_args, atmos) "LifeCycleTan2018", "ARM_SGP", "DYCOMS_RF01", - "DYCOMS_RF02", "Rico", "TRMM_LBA", "SimplePlume", @@ -389,6 +388,14 @@ function get_initial_condition(parsed_args, atmos) return getproperty(ICs, Symbol(parsed_args["initial_condition"]))( parsed_args["prognostic_tke"], ) + elseif parsed_args["initial_condition"] == "DYCOMS_RF02" + return getproperty(ICs, Symbol(parsed_args["initial_condition"]))( + parsed_args["prognostic_tke"], + parsed_args["q_tot_0_dycoms_rf02"], + parsed_args["theta_0_dycoms_rf02"], + parsed_args["theta_i_dycoms_rf02"], + parsed_args["z_i_dycoms_rf02"], + ) elseif parsed_args["initial_condition"] == "ISDAC" ICs.ISDAC( parsed_args["prognostic_tke"],