diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 8d938e44b4..dccc69b96d 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -281,6 +281,15 @@ steps: artifact_paths: "deep_sphere_baroclinic_wave_rhoe_equilmoist/output_active/*" agents: slurm_constraint: icelake|cascadelake|skylake|epyc + + - label: ":computer: test DYAMOND interpolated initial conditions" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $GPU_CONFIG_PATH/gpu_aquaplanet_dyamond_summer.yml + --job_id gpu_aquaplanet_dyamond_summer + artifact_paths: "gpu_aquaplanet_dyamond_summer/output_active/*" + agents: + slurm_constraint: icelake|cascadelake|skylake|epyc # Add this back when we figure out what to do with SSP and zalesak # - label: ":computer: SSP zalesak tracer & energy upwind baroclinic wave (ρe_tot) equilmoist" diff --git a/Artifacts.toml b/Artifacts.toml index 70062a55ca..1bbe36cfbc 100644 --- a/Artifacts.toml +++ b/Artifacts.toml @@ -43,3 +43,11 @@ git-tree-sha1 = "10742e0a2e343d13bb04df379e300a83402d4106" [[era5_cloud.download]] sha256 = "bb51e2f2d315b487e05a8d38944d4ad937ee4a40c43b68541220c5d54425e24a" url = "https://caltech.box.com/shared/static/b6ur4ap4vo04j09vdulem96z9fxqlgyn.gz" + +[DYAMOND_SUMMER_ICS_p98deg] +git-tree-sha1 = "46aa28a9177dd7e4cbb813eeba3e02597c3ba071" +lazy = true + + [[DYAMOND_SUMMER_ICS_p98deg.download]] + sha256 = "9a5efd20d68d9b954af2fd7803bccdda3fc7ef4ff60421ed13d18f768312d2e3" + url = "https://caltech.box.com/shared/static/whn01xr83v0s4p8tc59ds3nvfsfe6ebs.gz" diff --git a/NEWS.md b/NEWS.md index 97edf8545a..51f86e288e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,10 +5,29 @@ Main ------- ### Features +### Read initial conditions from NetCDF files + +Added functionality to allow initial conditions to be overwritten by +interpolated NetCDF datasets. + +To use this feature from the YAML interface, just pass the path of the file. +We expect the file to contain the following variables: +- `p`, for pressure, +- `t`, for temperature, +- `q`, for humidity, +- `u, v, w`, for velocity, +- `cswc, crwc` for snow and rain water content (for 1 moment microphysics). + +For example, to use the DYAMONDSummer initial condition, set +``` +initial_condition: "artifact\"DYAMONDSummer\"/DYAMOND_SUMMER_ICS_p98deg.nc" +``` +in your configuration file. + ### Write diagnostics to text files Added functionality to write diagnostics in DictWriter to text files. -This is useful for outputing scalar diagnostics, such as total mass of +This is useful for outputting scalar diagnostics, such as total mass of the atmosphere. PR [3476](https://github.com/CliMA/ClimaAtmos.jl/pull/3476) v0.28.0 diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml index 7acb4764bc..69d4a30d3a 100644 --- a/config/default_configs/default_config.yml +++ b/config/default_configs/default_config.yml @@ -182,7 +182,7 @@ surface_temperature: help: "Prescribed surface temperature functional form ['ZonallySymmetric' (default), 'ZonallyAsymmetric', 'RCEMIPII']" value: "ZonallySymmetric" initial_condition: - help: "Initial condition [`DryBaroclinicWave`, `MoistBaroclinicWave`, `DecayingProfile`, `IsothermalProfile`, `Bomex`, `DryDensityCurrentProfile`, `AgnesiHProfile`, `ScharProfile`, `RisingThermalBubbleProfile`, `ISDAC`]" + help: "Initial condition [`DryBaroclinicWave`, `MoistBaroclinicWave`, `DecayingProfile`, `IsothermalProfile`, `Bomex`, `DryDensityCurrentProfile`, `AgnesiHProfile`, `ScharProfile`, `RisingThermalBubbleProfile`, `ISDAC`], or a file path for a NetCDF file (read documentation about requirements)." value: "DecayingProfile" perturb_initstate: help: "Add a perturbation to the initial condition [`false`, `true` (default)]" diff --git a/config/gpu_configs/gpu_aquaplanet_dyamond_summer.yml b/config/gpu_configs/gpu_aquaplanet_dyamond_summer.yml new file mode 100644 index 0000000000..95cb7013de --- /dev/null +++ b/config/gpu_configs/gpu_aquaplanet_dyamond_summer.yml @@ -0,0 +1,31 @@ +dt_save_state_to_disk: "Inf" +dt_save_to_sol: "Inf" +output_default_diagnostics: false +h_elem: 16 +z_max: 60000.0 +z_elem: 63 +dz_bottom: 30.0 +rayleigh_sponge: true +viscous_sponge: true +moist: "equil" +precip_model: "1M" +rad: "allskywithclear" +insolation: "timevarying" +dt_rad: "1hours" +dt_cloud_fraction: "1hours" +vert_diff: "FriersonDiffusion" +implicit_diffusion: true +approximate_linear_solve_iters: 2 +surface_setup: "DefaultMoninObukhov" +dt: "5secs" +t_end: "30secs" +toml: [toml/longrun_aquaplanet.toml] +prescribe_ozone: true +aerosol_radiation: true +prescribed_aerosols: ["CB1", "CB2", "DST01", "OC1", "OC2", "SO4", "SSLT01"] +start_date: "20160801" +initial_condition: "artifact\"DYAMOND_SUMMER_ICS_p98deg\"/DYAMOND_SUMMER_ICS_p98deg.nc" +topography: "Earth" +diagnostics: + - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr, rv] + period: "30secs" diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl index da9252eff3..beefe9d664 100644 --- a/post_processing/ci_plots.jl +++ b/post_processing/ci_plots.jl @@ -679,6 +679,7 @@ end SphereOrographyPlots = Union{ Val{:sphere_baroclinic_wave_rhoe_topography_dcmip_rs}, + Val{:gpu_aquaplanet_dyamond_summer}, Val{:sphere_baroclinic_wave_rhoe_hughes2023}, } diff --git a/src/cache/tracer_cache.jl b/src/cache/tracer_cache.jl index 21e09eff36..3c4ef831df 100644 --- a/src/cache/tracer_cache.jl +++ b/src/cache/tracer_cache.jl @@ -1,4 +1,5 @@ import Dates: Year +import ClimaUtilities.TimeVaryingInputs import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput, LinearPeriodFillingInterpolation import Interpolations as Intp @@ -14,7 +15,10 @@ function ozone_cache(::PrescribedOzone, Y, start_date) reference_date = start_date, regridder_type = :InterpolationsRegridder, regridder_kwargs = (; extrapolation_bc), - method = LinearPeriodFillingInterpolation(Year(1)), + method = LinearPeriodFillingInterpolation( + Year(1), + TimeVaryingInputs.Flat(), + ), ) return (; o3, prescribed_o3_timevaryinginput) end diff --git a/src/initial_conditions/InitialConditions.jl b/src/initial_conditions/InitialConditions.jl index 73e0e58151..4d283036b4 100644 --- a/src/initial_conditions/InitialConditions.jl +++ b/src/initial_conditions/InitialConditions.jl @@ -9,8 +9,11 @@ import ..Microphysics0Moment import ..Microphysics1Moment import ..PrescribedSurfaceTemperature import ..PrognosticSurfaceTemperature +import ..ᶜinterp +import ..ᶠinterp import ..C3 import ..C12 +import ..compute_kinetic! import ..PrognosticEDMFX import ..DiagnosticEDMFX import ..n_mass_flux_subdomains @@ -30,6 +33,7 @@ import SciMLBase import Interpolations as Intp import NCDatasets as NC import Statistics: mean +import ClimaUtilities.SpaceVaryingInputs include("local_state.jl") include("atmos_state.jl") diff --git a/src/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl index c3f9edef99..8e61c1f81c 100644 --- a/src/initial_conditions/initial_conditions.jl +++ b/src/initial_conditions/initial_conditions.jl @@ -160,6 +160,33 @@ function (initial_condition::DecayingProfile)(params) return local_state end +""" + MoistFromFile(file_path) + +This function assigns an empty initial condition for , populating the `LocalState` with +`NaN`, and later overwriting it with the content of the given file +""" +struct MoistFromFile <: InitialCondition + file_path::String +end + +function (initial_condition::MoistFromFile)(params) + function local_state(local_geometry) + FT = eltype(params) + grav = CAP.grav(params) + thermo_params = CAP.thermodynamics_params(params) + + T, p = FT(NaN), FT(NaN) # placeholder values + + return LocalState(; + params, + geometry = local_geometry, + thermo_state = TD.PhaseDry_pT(thermo_params, p, T), + ) + end + return local_state +end + """ AgnesiHProfile(; perturb = false) @@ -351,10 +378,117 @@ function (initial_condition::RisingThermalBubbleProfile)(params) return local_state end +""" + overwrite_initial_conditions!(initial_condition, args...) + +Do-nothing fallback method for the operation overwriting initial conditions +(this functionality required in instances where we interpolate initial conditions from NetCDF files). +Future work may revisit this design choice. +""" +function overwrite_initial_conditions!( + initial_condition::InitialCondition, + args..., +) + return nothing +end + +""" + overwrite_initial_conditions!(initial_condition::MoistFromFile, Y, thermo_params, config) + +Given a prognostic state `Y`, an `initial condition` (specifically, where initial values are +assigned from interpolations of existing datasets), a `thermo_state`, this function +overwrites the default initial condition and populates prognostic variables with +interpolated values using the `SpaceVaryingInputs` tool. To mitigate issues related to +unbalanced states following the interpolation operation, we recompute vertical pressure +levels assuming hydrostatic balance, given the surface pressure. + +We expect the file to contain the following variables: +- `p`, for pressure, +- `t`, for temperature, +- `q`, for humidity, +- `u, v, w`, for velocity, +- `cswc, crwc` for snow and rain water content (for 1 moment microphysics). +""" +function overwrite_initial_conditions!( + initial_conditions::MoistFromFile, + Y, + thermo_params, +) + file_path = initial_conditions.file_path + isfile(file_path) || error("$(file_path) is not a file") + @info "Overwriting initial conditions with data from file $(file_path)" + center_space = Fields.axes(Y.c) + face_space = Fields.axes(Y.f) + # Using surface pressure, air temperature and specific humidity + # from the dataset, compute air pressure. + p_sfc = Fields.level( + SpaceVaryingInputs.SpaceVaryingInput(file_path, "p", face_space), + Fields.half, + ) + ᶜT = SpaceVaryingInputs.SpaceVaryingInput(file_path, "t", center_space) + ᶜq_tot = SpaceVaryingInputs.SpaceVaryingInput(file_path, "q", center_space) + + # With the known temperature (ᶜT) and moisture (ᶜq_tot) profile, + # recompute the pressure levels assuming hydrostatic balance is maintained. + # Uses the ClimaCore `column_integral_indefinite!` function to solve + # ∂(ln𝑝)/∂z = -g/(Rₘ(q)T), where + # p is the local pressure + # g is the gravitational constant + # q is the specific humidity + # Rₘ is the gas constant for moist air + # T is the air temperature + # p is then updated with the integral result, given p_sfc, + # following which the thermodynamic state is constructed. + ᶜ∂lnp∂z = @. -thermo_params.grav / + (TD.gas_constant_air(thermo_params, TD.PhasePartition(ᶜq_tot)) * ᶜT) + ᶠlnp_over_psfc = zeros(face_space) + Operators.column_integral_indefinite!(ᶠlnp_over_psfc, ᶜ∂lnp∂z) + ᶠp = p_sfc .* exp.(ᶠlnp_over_psfc) + ᶜts = TD.PhaseEquil_pTq.(thermo_params, ᶜinterp.(ᶠp), ᶜT, ᶜq_tot) + + # Assign prognostic variables from equilibrium moisture models + Y.c.ρ .= TD.air_density.(thermo_params, ᶜts) + # Velocity is first assigned on cell-centers and then interpolated onto + # cell faces. + vel = + Geometry.UVWVector.( + SpaceVaryingInputs.SpaceVaryingInput(file_path, "u", center_space), + SpaceVaryingInputs.SpaceVaryingInput(file_path, "v", center_space), + SpaceVaryingInputs.SpaceVaryingInput(file_path, "w", center_space), + ) + Y.c.uₕ .= C12.(Geometry.UVVector.(vel)) + Y.f.u₃ .= ᶠinterp.(C3.(Geometry.WVector.(vel))) + e_kin = similar(ᶜT) + compute_kinetic!(e_kin, Y.c.uₕ, Y.f.u₃) + e_pot = Fields.coordinate_field(Y.c).z .* thermo_params.grav + Y.c.ρe_tot .= TD.total_energy.(thermo_params, ᶜts, e_kin, e_pot) .* Y.c.ρ + if hasproperty(Y.c, :ρq_tot) + Y.c.ρq_tot .= ᶜq_tot .* Y.c.ρ + else + error( + "`dry` configurations are incompatible with the interpolated initial conditions.", + ) + end + if hasproperty(Y.c, :ρq_sno) && hasproperty(Y.c, :ρq_rai) + Y.c.ρq_sno .= + SpaceVaryingInputs.SpaceVaryingInput( + file_path, + "cswc", + center_space, + ) .* Y.c.ρ + Y.c.ρq_rai .= + SpaceVaryingInputs.SpaceVaryingInput( + file_path, + "crwc", + center_space, + ) .* Y.c.ρ + end + return nothing +end + ## ## Baroclinic Wave ## - function shallow_atmos_baroclinic_wave_values(z, ϕ, λ, params, perturb) FT = eltype(params) R_d = CAP.R_d(params) diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index 51cde77ff1..fde49a911c 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -296,6 +296,8 @@ function get_initial_condition(parsed_args) "PrecipitatingColumn", ] return getproperty(ICs, Symbol(parsed_args["initial_condition"]))() + elseif isfile(parsed_args["initial_condition"]) + return ICs.MoistFromFile(parsed_args["initial_condition"]) elseif parsed_args["initial_condition"] == "GCM" @assert parsed_args["prognostic_tke"] == true return ICs.GCMDriven( @@ -623,7 +625,6 @@ end function get_simulation(config::AtmosConfig) params = create_parameter_set(config) atmos = get_atmos(config, params) - sim_info = get_sim_info(config) job_id = sim_info.job_id output_dir = sim_info.output_dir @@ -652,7 +653,6 @@ function get_simulation(config::AtmosConfig) initial_condition = get_initial_condition(config.parsed_args) surface_setup = get_surface_setup(config.parsed_args) - if !sim_info.restart s = @timed_str begin Y = ICs.atmos_state( @@ -664,6 +664,17 @@ function get_simulation(config::AtmosConfig) t_start = Spaces.undertype(axes(Y.c))(0) end @info "Allocating Y: $s" + + # In instances where we wish to interpolate existing datasets, e.g. + # NetCDF files containing spatially varying thermodynamic properties, + # this call to `overwrite_initial_conditions` overwrites the variables + # in `Y` (specific to `initial_condition`) with those computed using the + # `SpaceVaryingInputs` tool. + CA.InitialConditions.overwrite_initial_conditions!( + initial_condition, + Y, + params.thermodynamics_params, + ) end tracers = get_tracers(config.parsed_args) diff --git a/test/runtests.jl b/test/runtests.jl index 5181583738..2c46493029 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,6 +24,7 @@ using Test @safetestset "Topography tests" begin @time include("topography.jl") end @safetestset "Restarts" begin @time include("restart.jl") end @safetestset "Reproducibility infra" begin @time include("unit_reproducibility_infra.jl") end +@safetestset "Init with file" begin @time include("test_init_with_file.jl") end #! format: on diff --git a/test/test_init_with_file.jl b/test/test_init_with_file.jl new file mode 100644 index 0000000000..f62434fa99 --- /dev/null +++ b/test/test_init_with_file.jl @@ -0,0 +1,14 @@ +import ClimaAtmos as CA + +simulation = CA.get_simulation( + CA.AtmosConfig( + Dict( + "initial_condition" => "artifact\"DYAMOND_SUMMER_ICS_p98deg\"/DYAMOND_SUMMER_ICS_p98deg.nc", + "moist" => "equil", + ), + job_id = "test_init_with_file_dyamond", + ), +) + +# Just a small test to see that we got here +@test maximum(simulation.integrator.u.c.ρ) > 0