LES exp ClimaAtmos (updated CC to use cartesian gradient evaluations)

	modified:   src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl
	modified:   config/model_configs/les_rico_box.yml
	modified:   post_processing/ci_plots.jl
	modified:   config/model_configs/les_bomex_box.yml
Tested BOMEX (1M) + RICO (1M)
	modified:   src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl

Clean up before review

config/default_configs/default_config.yml

@@ -142,6 +142,10 @@ vert_diff:
 hyperdiff:
   help: "Hyperdiffusion [`ClimaHyperdiffusion` (or `true`; default), `none` (or `false`)]"
   value: "CAM_SE"
+smagorinsky_lilly:
+  value: false
+c_smag: 
+  value: 0.2
 bubble:
   help: "Enable bubble correction for more accurate surface areas"
   value: true

config/model_configs/box_density_current_test.yml

@@ -1,16 +1,22 @@
+job_id: "box_density_current_test"
 dt_save_state_to_disk: "10secs"
 initial_condition: "DryDensityCurrentProfile"
 x_max: 51200.0
 z_elem: 45
-dt: "0.1secs"
-t_end: "10.0secs"
+dt: "0.3secs"
+t_end: "1000secs"
 discrete_hydrostatic_balance: true
-y_max: 6400.0
-y_elem: 15
+hyperdiff: "false"
+y_max: 51200.0
+y_elem: 45
 z_stretch: false
 x_elem: 45
 config: "box"
+hyperdiff: "true"
 z_max: 6400.0
+smagorinsky_lilly: true
+c_smag: 0.25
+netcdf_interpolation_num_points: [40, 40, 80]
 diagnostics:
-  - short_name: thetaa
-    period: 10secs
+  - short_name: [thetaa, wa, ua, ta, va, hfes, rhoa, rv, ha]
+    period: 30secs

config/model_configs/les_bomex_box.yml

@@ -0,0 +1,36 @@
+initial_condition: "Bomex" 
+edmf_coriolis: "Bomex"
+subsidence: "Bomex" 
+ls_adv: "Bomex" 
+surface_setup: "Bomex"
+moist: "equil"
+config: "box" 
+hyperdiff: "false" 
+implicit_diffusion: false
+precip_model: "1M"
+dt_cloud_fraction: "10mins"
+ode_algo: "SSPKnoth"
+smagorinsky_lilly: true
+call_cloud_diagnostics_per_stage: true
+c_smag: 0.17
+x_max: 3.2e3
+y_max: 3.2e3 
+z_max: 4.5e3
+x_elem: 8
+y_elem: 8
+z_elem: 100
+z_stretch: false
+dt: "0.2secs" 
+output_default_diagnostics: false
+t_end: "21600secs" 
+dt_save_state_to_disk: "30mins"
+rayleigh_sponge: true
+toml: [toml/bomex_les.toml]
+restart_file: "/home/asridhar/ClimaAtmos.jl/restart/bomex/Micro1M_day0.0_Prod.hdf5"
+netcdf_interpolation_num_points: [40, 40, 80]
+diagnostics:
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, hfes]
+    period: 10mins
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, hfes]
+    reduction: average
+    period: 60mins

config/model_configs/les_dycoms_box.yml

@@ -0,0 +1,36 @@
+job_id: "les_dycoms_box"
+initial_condition: DYCOMS_RF01 
+edmf_coriolis: DYCOMS_RF01
+subsidence: DYCOMS 
+rad: DYCOMS
+surface_setup: DYCOMS_RF01
+moist: "equil"
+config: box 
+hyperdiff: "false"
+ode_algo: "SSPKnoth"
+precip_model: "1M"
+x_max: 500.0
+y_max: 500.0
+z_max: 1.5e3
+x_elem: 4
+y_elem: 4
+z_elem: 300
+z_stretch: false
+smagorinsky_lilly: true
+tracer_upwinding: third_order
+c_smag: 0.15
+dt: "0.05secs"
+t_end: "14400secs"
+dt_save_state_to_disk: "30mins"
+restart_file: "/home/asridhar/Codes/ClimaAtmos.jl/restart/dycoms/day0_1m.hdf5"
+#restart_file: "/home/asridhar/Codes/ClimaAtmos.jl/restart/dycoms/resolution_10m_day0.0.hdf5"
+dt_cloud_fraction: "10mins"
+rayleigh_sponge: true
+toml: [toml/dycoms_box_rhoe.toml]
+netcdf_interpolation_num_points: [30,30,75]
+diagnostics:
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, hfes]
+    period: 10mins
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, hfes]
+    reduction_time: average
+    period: 60mins

config/model_configs/les_gabls_box.yml

@@ -0,0 +1,26 @@
+job_id: "les_gabls_box"
+edmf_coriolis: GABLS
+initial_condition: GABLS
+surface_setup: GABLS
+implicit_diffusion: false
+moist: "equil"
+precip_model: "1M"
+config: "box"
+hyperdiff: "false"
+x_max: 1e3
+y_max: 1e3
+z_max: 500
+x_elem: 5
+y_elem: 5
+z_elem: 10
+z_stretch: false
+dt: "0.10secs"
+t_end: "100secs"
+dt_save_state_to_disk: "30mins"
+perturb_initstate: true
+netcdf_output_at_levels: true
+netcdf_interpolation_num_points: [20, 20, 20]
+#restart_file: "restart/gabls/day0.0.hdf5"
+diagnostics:
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+    period: 10secs

config/model_configs/les_rico_box.yml

@@ -0,0 +1,37 @@
+job_id: les_rico_box
+initial_condition: Rico
+subsidence: Rico
+edmf_coriolis: Rico
+ls_adv: Rico
+surface_setup: Rico
+implicit_diffusion: false
+approximate_linear_solve_iters: 2
+moist: equil
+hyperdiff: "false"
+apply_limiter: false
+precip_model: "1M"
+config: box
+dt_cloud_fraction: "10mins"
+ode_algo: "SSPKnoth"
+smagorinsky_lilly: true
+c_smag: 0.17
+x_max: 4e3
+y_max: 4e3
+x_elem: 10
+y_elem: 10
+z_elem: 150
+z_max: 6e3
+z_stretch: false
+rayleigh_sponge: true
+dt: "0.20secs"
+t_end: 12hours
+restart_file: "/home/asridhar/Codes/ClimaAtmos.jl/restart/rico/day0.0.hdf5"
+toml: [toml/rico_les.toml]
+dt_save_state_to_disk: 30mins
+netcdf_interpolation_num_points: [30, 30, 150]
+diagnostics:
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, husra, cl, clw, cli, hussfc, evspsbl, hfes, pr]
+    period: 10mins
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, husra, cl, clw, cli, hussfc, evspsbl, hfes, pr]
+    reduction_time: average
+    period: 60mins

post_processing/ci_plots.jl

@@ -42,7 +42,7 @@ import CairoMakie
 import CairoMakie.Makie
 import ClimaAnalysis
 import ClimaAnalysis: Visualize as viz
-import ClimaAnalysis: SimDir, slice, read_var, average_xy
+import ClimaAnalysis: SimDir, slice, read_var, average_xy, window, average_time
 import ClimaAnalysis.Utils: kwargs as ca_kwargs
 
 import ClimaCoreSpectra: power_spectrum_2d
@@ -577,6 +577,7 @@ end
 
 function make_plots(
     ::Val{:box_density_current_test},
+    ::Val{:box_rising_thermal_test},
     output_paths::Vector{<:AbstractString},
 )
     simdirs = SimDir.(output_paths)
@@ -1065,6 +1066,117 @@ function make_plots(::Aquaplanet1MPlots, output_paths::Vector{<:AbstractString})
     )
 end
 
+LESBoxPlots= Union{
+    Val{:les_rico_box},
+    Val{:les_dycoms_box},
+    Val{:les_bomex_box},
+    Val{:les_gabls_box},
+}
+
+"""
+    plot_les_vert_profile!(grid_loc, var_group)
+
+Helper function for `make_plots_generic`. Takes a list of variables and plots
+them on the same axis.
+"""
+function plot_les_vert_profile!(grid_loc, var_group)
+    z = var_group[1].dims["z"]
+    units = var_group[1].attributes["units"]
+    ax = CairoMakie.Axis(
+        grid_loc[1, 1],
+        ylabel = "z [$(var_group[1].dim_attributes["z"]["units"])]",
+        xlabel = "$(short_name(var_group[1])) [$units]",
+        title = parse_var_attributes(var_group[1]),
+    )
+
+    for var in var_group
+        CairoMakie.lines!(ax, var.data, z, label = short_name(var))
+    end
+    length(var_group) > 1 && Makie.axislegend(ax)
+end
+
+function make_plots(
+    sim_type::Union{LESBoxPlots},
+    output_paths::Vector{<:AbstractString},
+)
+    simdirs = SimDir.(output_paths)
+
+    reduction = "inst"
+    short_names = [
+        "wa",
+        "ua",
+        "va",
+        "ta",
+        "thetaa",
+        "ha",
+        "hus",
+        "hur",
+        "cl",
+        "clw",
+        "cli",
+    ]
+    available_periods = ClimaAnalysis.available_periods(
+        simdirs[1];
+        short_name = short_names[1],
+        reduction,
+    )
+    if "5m" in available_periods
+        period = "5m"
+    elseif "10m" in available_periods
+        period = "10m"
+    elseif "30m" in available_periods
+        period = "30m"
+    elseif "1h" in available_periods
+        period = "1h"
+    end
+
+    # Window average from instantaneous snapshots?
+    function horizontal_average(var)
+        return average_xy(var)
+    end
+    function windowed_reduction(var)
+        hours = 3600.0
+        window_end = last(var.dims["time"])
+        window_start = window_end - 2hours
+        var_window = ClimaAnalysis.window(var, "time"; left=window_start, right=window_end)
+        var_reduced = horizontal_average(average_time(var_window))
+        return var_reduced
+    end
+
+    var_groups_xyt_reduced =
+        map_comparison(simdirs, short_names) do simdir, short_name
+            return [
+                get(simdir; short_name, reduction, period) |> windowed_reduction
+            ]
+        end
+
+    var_groups_xy_reduced =
+        map_comparison(simdirs, short_names) do simdir, short_name
+            return [
+                get(simdir; short_name, reduction, period) |> horizontal_average
+            ]
+        end
+
+    tmp_file = make_plots_generic(
+        output_paths,
+        var_groups_xyt_reduced,
+        output_name = "tmp";
+        plot_fn = plot_les_vert_profile!,
+        MAX_NUM_COLS = 2,
+        MAX_NUM_ROWS = 4,
+    )
+
+    make_plots_generic(
+        output_paths,
+        vcat(var_groups_xy_reduced...),
+        plot_fn = plot_parsed_attribute_title!,
+        summary_files = [tmp_file],
+        MAX_NUM_COLS = 2,
+        MAX_NUM_ROWS = 4,
+    )
+end
+
+
 EDMFBoxPlots = Union{
     Val{:diagnostic_edmfx_gabls_box},
     Val{:diagnostic_edmfx_bomex_box},
@@ -1097,7 +1209,6 @@ EDMFBoxPlotsWithPrecip = Union{
     Val{:diagnostic_edmfx_trmm_stretched_box},
 }
 
-
 """
     plot_edmf_vert_profile!(grid_loc, var_group)
 

src/ClimaAtmos.jl

@@ -110,6 +110,7 @@ include(
 )
 include(joinpath("parameterized_tendencies", "sponge", "rayleigh_sponge.jl"))
 include(joinpath("parameterized_tendencies", "sponge", "viscous_sponge.jl"))
+include(joinpath("parameterized_tendencies", "les_sgs_models", "smagorinsky_lilly.jl"))
 include(joinpath("prognostic_equations", "advection.jl"))
 
 include(joinpath("cache", "temporary_quantities.jl"))

src/cache/cache.jl

@@ -16,6 +16,7 @@ struct AtmosCache{
     DSS,
     RS,
     VS,
+    SL, 
     PR,
     SUB,
     LSAD,
@@ -79,6 +80,7 @@ struct AtmosCache{
     """Additional parameters used by the various tendencies"""
     rayleigh_sponge::RS
     viscous_sponge::VS
+    smagorinsky_lilly::SL
     precipitation::PR
     subsidence::SUB
     large_scale_advection::LSAD
@@ -187,6 +189,7 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names)
     hyperdiff = hyperdiffusion_cache(Y, atmos)
     rayleigh_sponge = rayleigh_sponge_cache(Y, atmos)
     viscous_sponge = viscous_sponge_cache(Y, atmos)
+    smagorinsky_lilly = smagorinsky_lilly_cache(Y, atmos)
     precipitation = precipitation_cache(Y, atmos)
     subsidence = subsidence_cache(Y, atmos)
     large_scale_advection = large_scale_advection_cache(Y, atmos)
@@ -216,6 +219,7 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names)
         do_dss,
         rayleigh_sponge,
         viscous_sponge,
+        smagorinsky_lilly,
         precipitation,
         subsidence,
         large_scale_advection,

src/cache/temporary_quantities.jl

@@ -53,6 +53,14 @@ function temporary_quantities(Y, atmos)
             typeof(UVW(FT(0), FT(0), FT(0)) * UVW(FT(0), FT(0), FT(0))'),
             face_space,
         ), # ᶠstrain_rate
+        ᶜtemp_strain = Fields.Field(
+            typeof(UVW(FT(0), FT(0), FT(0)) * UVW(FT(0), FT(0), FT(0))'),
+            center_space,
+        ), # ᶜstrain_rate
+        ᶠtemp_strain = Fields.Field(
+            typeof(UVW(FT(0), FT(0), FT(0)) * UVW(FT(0), FT(0), FT(0))'),
+            face_space,
+        ), # ᶠstrain_rate
         # TODO: Remove this hack
         sfc_temp_C3 = Fields.Field(C3{FT}, Spaces.level(face_space, half)), # ρ_flux_χ
         # Implicit solver cache: