Merge pull request #2567 from CliMA/ln/precip-energy

Track precipitation energy in conservation check

Author: LenkaNovak

.buildkite/pipeline.yml

@@ -123,7 +123,7 @@ steps:
       - label: ":computer: orographic gravity wave parameterization unit test for 3d calculation"
         command: "julia --color=yes --project=examples test/orographic_gravity_wave/ogwd_3d.jl"
         artifact_paths: "orographic_gravity_wave_test_3d/*"
-      
+
       - label: ":computer: single column non-orographic gravity wave parameterization"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
@@ -294,6 +294,15 @@ steps:
           slurm_mem: 20GB
           slurm_constraint: icelake|cascadelake|skylake|epyc
 
+      - label: ":computer: aquaplanet (ρe_tot) slabocean equilmoist allsky radiation monin_obukhov varying insolation"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $CONFIG_PATH/aquaplanet_rhoe_equil_clearsky_tvinsol_0M_slabocean.yml
+        artifact_paths: "aquaplanet_rhoe_equil_clearsky_tvinsol_0M_slabocean/*"
+        agents:
+          slurm_mem: 20GB
+          slurm_constraint: icelake|cascadelake|skylake|epyc
+
   - group: "Sphere Examples (Topography)"
     steps:
 

config/default_configs/default_config.yml

@@ -122,7 +122,7 @@ reference_job_id:
     Identifier of job to use as the "reference" solution in the quicklook plot; the current job's results get compared to the results of the quicklook job on the main branch
   value: ~
 vert_diff:
-  help: "Vertical diffusion [`false` (default), `VerticalDiffusion`, `true` (defaults to `VerticalDiffusion`)]"
+  help: "Vertical diffusion [`false` (default), `VerticalDiffusion`, `true` (defaults to `VerticalDiffusion`), `FriersonDiffusion`]"
   value: "false"
 hyperdiff:
   help: "Hyperdiffusion [`ClimaHyperdiffusion` (or `true`; default), `none` (or `false`)]"

config/longrun_configs/longrun_aquaplanet_amip.yml

@@ -3,28 +3,28 @@ z_max: 55000.0
 z_elem: 63
 dz_bottom: 30.0
 dz_top: 3000.0
-moist: "equil" 
-precip_model: "0M" 
+moist: "equil"
+precip_model: "0M"
 override_τ_precip: false
-rad: "allskywithclear" 
+rad: "allskywithclear"
 idealized_insolation: false
 dt_rad: "1hours"
-surface_setup: "DefaultMoninObukhov" 
-turbconv: "diagnostic_edmfx" 
+surface_setup: "DefaultMoninObukhov"
+turbconv: "diagnostic_edmfx"
 implicit_diffusion: true
 approximate_linear_solve_iters: 2
 prognostic_tke: true
-edmfx_upwinding: "first_order" 
-edmfx_entr_model: "Generalized" 
-edmfx_detr_model: "Generalized" 
-edmfx_nh_pressure: true 
+edmfx_upwinding: "first_order"
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
+edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 rayleigh_sponge: true
 dt_save_state_to_disk: "1days"
-dt: "100secs" 
+dt: "100secs"
 t_end: "30days"
-job_id: "longrun_aquaplanet_amip" 
+job_id: "longrun_aquaplanet_amip"
 toml: [toml/longrun_aquaplanet_amip.toml]
 output_default_diagnostics: false
 diagnostics:

config/longrun_configs/longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_tvinsol_0M_slabocean.yml

@@ -21,3 +21,8 @@ check_conservation: true
 bubble: false
 job_id: "longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_tvinsol_0M_slabocean"
 toml: [toml/longrun_aquaplanet_rhoe_equil_55km_nz63_gray_0M.toml]
+output_default_diagnostics: false
+diagnostics:
+  - short_name: ["ts", "ta", "thetaa", "ha", "pfull", "rhoa", "ua", "va", "wa", "hfes", "hur", "hus", "clw", "cli", "evspsbl", "rsd", "rsu", "rld", "rlu"]
+    reduction_time: average
+    period: 10days

config/longrun_configs/longrun_aquaplanet_rhoe_equil_clearsky_tvinsol_0M_slabocean.yml

@@ -13,3 +13,8 @@ prognostic_surface: "PrognosticSurfaceTemperature"
 check_conservation: true
 bubble: false
 job_id: "longrun_aquaplanet_rhoe_equil_clearsky_tvinsol_0M_slabocean"
+output_default_diagnostics: false
+diagnostics:
+  - short_name: ["ts", "ta", "thetaa", "ha", "pfull", "rhoa", "ua", "va", "wa", "hfes", "hur", "hus", "clw", "cli", "evspsbl", "rsd", "rsu", "rld", "rlu"]
+    reduction_time: average
+    period: 10days

config/longrun_configs/longrun_sphere_hydrostatic_balance_rhoe.yml

@@ -2,8 +2,8 @@ h_elem: 16
 z_elem: 45
 dz_bottom: 30.0
 dt: "400secs"
-perturb_initstate: false 
-discrete_hydrostatic_balance: true 
-t_end: "800days" 
-dt_save_state_to_disk: "10days" 
+perturb_initstate: false
+discrete_hydrostatic_balance: true
+t_end: "800days"
+dt_save_state_to_disk: "10days"
 job_id: "longrun_sphere_hydrostatic_balance_rhoe"

config/longrun_configs/longrun_ssp_bw_rhoe_equil_highres.yml

@@ -1,6 +1,6 @@
 dt_save_state_to_disk: "10days"
 initial_condition: "MoistBaroclinicWave"
-max_newton_iters_ode: 3
+approximate_linear_solve_iters: 2
 z_elem: 45
 dt: "150secs"
 t_end: "100days"

config/model_configs/aquaplanet_rhoe_equil_clearsky_tvinsol_0M_slabocean.yml

@@ -0,0 +1,21 @@
+FLOAT_TYPE: "Float32"
+dt_save_state_to_disk: "10days"
+t_end: "10days"
+moist: "equil"
+precip_model: "0M"
+surface_setup: "DefaultMoninObukhov"
+vert_diff: "FriersonDiffusion"
+implicit_diffusion: true
+approximate_linear_solve_iters: 2
+rad: "clearsky"
+idealized_insolation: false
+dt_rad: "1hours"
+prognostic_surface: "PrognosticSurfaceTemperature"
+check_conservation: true
+bubble: false
+job_id: "aquaplanet_rhoe_equil_clearsky_tvinsol_0M_slabocean"
+output_default_diagnostics: false
+diagnostics:
+  - short_name: ["ts", "ta", "thetaa", "ha", "pfull", "rhoa", "ua", "va", "wa", "hfes", "hur", "hus", "clw", "cli", "evspsbl", "rsd", "rsu", "rld", "rlu"]
+    reduction_time: average
+    period: 10days

config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml

@@ -6,11 +6,10 @@ surface_setup: TRMM_LBA
 turbconv: diagnostic_edmfx
 implicit_diffusion: true
 approximate_linear_solve_iters: 2
-max_newton_iters_ode: 1
 prognostic_tke: true
 edmfx_upwinding: first_order
-edmfx_entr_model: "Generalized" 
-edmfx_detr_model: "Generalized" 
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
 edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true

config/perf_configs/flame_perf_target_frierson.yml

@@ -1,5 +1,5 @@
 vert_diff: "FriersonDiffusion"
-max_newton_iters_ode: 2
+pproximate_linear_solve_iters: 2
 use_krylov_method: true
 use_dynamic_krylov_rtol: false
 krylov_rtol: 0.99

config/perf_configs/flame_sphere_baroclinic_wave_rhoe_equilmoist_expvdiff.yml

@@ -1,4 +1,3 @@
-max_newton_iters_ode: 1
 use_krylov_method: true
 use_dynamic_krylov_rtol: false
 krylov_rtol: 0.99

examples/hybrid/driver.jl

@@ -1,5 +1,5 @@
 # When Julia 1.10+ is used interactively, stacktraces contain reduced type information to make them shorter.
-# On the other hand, the full type information is printed when julia is not run interactively. 
+# On the other hand, the full type information is printed when julia is not run interactively.
 # Given that ClimaCore objects are heavily parametrized, non-abbreviated stacktraces are hard to read,
 # so we force abbreviated stacktraces even in non-interactive runs.
 # (See also Base.type_limited_string_from_context())
@@ -58,29 +58,7 @@ if !isempty(integrator.tstops)
 end
 CA.verify_callbacks(sol.t)
 
-if ClimaComms.iamroot(config.comms_ctx)
-    @info "Plotting"
-    make_plots(Val(Symbol(reference_job_id)), simulation.output_dir)
-    @info "Plotting done"
-
-    @info "Creating tarballs"
-    Tar.create(
-        f -> endswith(f, ".nc"),
-        simulation.output_dir,
-        joinpath(simulation.output_dir, "nc_files.tar"),
-    )
-    Tar.create(
-        f -> endswith(f, r"hdf5|h5"),
-        simulation.output_dir,
-        joinpath(simulation.output_dir, "hdf5_files.tar"),
-    )
-
-    foreach(readdir(simulation.output_dir)) do f
-        endswith(f, r"nc|hdf5|h5") && rm(joinpath(simulation.output_dir, f))
-    end
-    @info "Tarballs created"
-end
-
+# Scaling check
 if CA.is_distributed(config.comms_ctx)
     nprocs = ClimaComms.nprocs(config.comms_ctx)
     comms_ctx = config.comms_ctx
@@ -105,6 +83,7 @@ if CA.is_distributed(config.comms_ctx)
     end
 end
 
+# Check if selected output has changed from the previous recorded output (bit-wise comparison)
 include(joinpath(@__DIR__, "..", "..", "regression_tests", "mse_tables.jl"))
 if config.parsed_args["regression_test"]
     # Test results against main branch
@@ -135,7 +114,9 @@ end
 @info "Callback verification, n_expected_calls: $(CA.n_expected_calls(integrator))"
 @info "Callback verification, n_measured_calls: $(CA.n_measured_calls(integrator))"
 
+# Conservation checks
 if config.parsed_args["check_conservation"]
+    @info "Checking conservation"
     FT = Spaces.undertype(axes(sol.u[end].c.ρ))
 
     # energy
@@ -154,8 +135,14 @@ if config.parsed_args["check_conservation"]
         energy_surface_change = -p.net_energy_flux_sfc[][]
     end
     energy_radiation_input = -p.net_energy_flux_toa[][]
-    @test (energy_atmos_change + energy_surface_change) / energy_total ≈
-          energy_radiation_input / energy_total atol = 5 * sqrt(eps(FT))
+
+    energy_net =
+        abs(
+            energy_atmos_change + energy_surface_change -
+            energy_radiation_input,
+        ) / energy_total
+    @info "    Net energy change: $energy_net"
+    @test (energy_net / energy_total) ≈ 0 atol = sqrt(eps(FT))
 
     if p.atmos.moisture_model isa CA.DryModel
         # density
@@ -169,12 +156,16 @@ if config.parsed_args["check_conservation"]
             water_surface_change = CA.horizontal_integral_at_boundary(
                 sol.u[end].sfc.water .- sol.u[1].sfc.water,
             )
-            @test (water_atmos_change + water_surface_change) / water_total ≈ 0 atol =
-                100 * sqrt(eps(FT))
+
+            water_net =
+                abs(water_atmos_change + water_surface_change) / water_total
+            @info "    Net water change: $water_net"
+            @test water_net ≈ 0 atol = 100 * sqrt(eps(FT))
         end
     end
 end
 
+# Precipitation characteristic checks
 if config.parsed_args["check_precipitation"]
     # run some simple tests based on the output
     FT = Spaces.undertype(axes(sol.u[end].c.ρ))
@@ -227,3 +218,27 @@ if config.parsed_args["check_precipitation"]
         @test maximum(sol.prob.p.precomputed.ᶜcloud_fraction[colidx]) <= FT(1)
     end
 end
+
+# Visualize the solution
+if ClimaComms.iamroot(config.comms_ctx)
+    @info "Plotting"
+    make_plots(Val(Symbol(reference_job_id)), simulation.output_dir)
+    @info "Plotting done"
+
+    @info "Creating tarballs"
+    Tar.create(
+        f -> endswith(f, ".nc"),
+        simulation.output_dir,
+        joinpath(simulation.output_dir, "nc_files.tar"),
+    )
+    Tar.create(
+        f -> endswith(f, r"hdf5|h5"),
+        simulation.output_dir,
+        joinpath(simulation.output_dir, "hdf5_files.tar"),
+    )
+
+    foreach(readdir(simulation.output_dir)) do f
+        endswith(f, r"nc|hdf5|h5") && rm(joinpath(simulation.output_dir, f))
+    end
+    @info "Tarballs created"
+end

post_processing/ci_plots.jl

@@ -475,6 +475,59 @@ function make_plots(
     )
 end
 
+function make_plots(
+    ::Union{
+        Val{:aquaplanet_rhoe_equil_clearsky_tvinsol_0M_slabocean},
+        Val{:longrun_aquaplanet_rhoe_equil_clearsky_tvinsol_0M_slabocean},
+        Val{
+            :longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_tvinsol_0M_slabocean,
+        },
+    },
+    simulation_path,
+)
+    simdir = SimDir(simulation_path)
+
+
+    reduction = "average"
+    period = "10d"
+    short_names_3D = [
+        "ta",
+        "thetaa",
+        "rhoa",
+        "ua",
+        "va",
+        "wa",
+        "hur",
+        "hus",
+        "clw",
+        "cli",
+        "rsd",
+        "rsu",
+        "rld",
+        "rlu",
+    ]
+    short_names_sfc = ["hfes", "evspsbl", "ts"]
+    vars_3D = [
+        get(simdir; short_name, reduction, period) |> ClimaAnalysis.average_lon for short_name in short_names_3D
+    ]
+    vars_sfc = [
+        get(simdir; short_name, reduction, period) for
+        short_name in short_names_sfc
+    ]
+    make_plots_generic(
+        simulation_path,
+        vars_3D,
+        time = LAST_SNAP,
+        more_kwargs = YLOGSCALE,
+    )
+    make_plots_generic(
+        simulation_path,
+        vars_sfc,
+        time = LAST_SNAP,
+        output_name = "summary_sfc",
+    )
+end
+
 AquaplanetPlots = Union{
     Val{:sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res},
     Val{:mpi_sphere_aquaplanet_rhoe_equilmoist_clearsky},
@@ -483,10 +536,8 @@ AquaplanetPlots = Union{
     Val{:longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_diagedmf_diffonly_0M},
     Val{:longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_diagedmf_0M},
     Val{:longrun_aquaplanet_rhoe_equil_55km_nz63_allsky_diagedmf_0M},
-    Val{:longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_tvinsol_0M_slabocean},
     Val{:longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_tvinsol_0M_earth},
     Val{:longrun_aquaplanet_rhoe_equil_highres_allsky_ft32},
-    Val{:longrun_aquaplanet_rhoe_equil_clearsky_tvinsol_0M_slabocean},
     Val{:longrun_aquaplanet_dyamond},
     Val{:longrun_aquaplanet_amip},
 }

src/cache/cache.jl

@@ -26,6 +26,7 @@ struct AtmosCache{
     RAD,
     NETFLUXTOA,
     NETFLUXSFC,
+    CONSCHECK,
     OD,
 }
     """Timestep of the simulation (in seconds). This is also used by callbacks and tendencies"""
@@ -89,6 +90,9 @@ struct AtmosCache{
     net_energy_flux_toa::NETFLUXTOA
     net_energy_flux_sfc::NETFLUXSFC
 
+    """Conservation check for prognostic surface temperature"""
+    conservation_check::CONSCHECK
+
     """Directory output."""
     output_dir::OD
 end
@@ -142,6 +146,14 @@ function build_cache(Y, atmos, params, surface_setup, sim_info)
     net_energy_flux_toa = [Geometry.WVector(FT(0))]
     net_energy_flux_sfc = [Geometry.WVector(FT(0))]
 
+    conservation_check =
+        !(atmos.precip_model isa NoPrecipitation) ?
+        (;
+            col_integrated_precip_energy_tendency = zeros(
+                axes(Fields.level(Geometry.WVector.(Y.f.u₃), half)),
+            )
+        ) : (; col_integrated_precip_energy_tendency = (;))
+
     limiter = if isnothing(atmos.numerics.limiter)
         nothing
     elseif atmos.numerics.limiter isa QuasiMonotoneLimiter
@@ -224,6 +236,7 @@ function build_cache(Y, atmos, params, surface_setup, sim_info)
         radiation,
         net_energy_flux_toa,
         net_energy_flux_sfc,
+        conservation_check,
         output_dir,
     )