diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index dc2925531ad..61e64b8b7f1 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -194,6 +194,33 @@ steps: --job_id plane_density_current_test artifact_paths: "plane_density_current_test/output_active/*" + - label: ":computer: Analytic No Topography Test (2D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_no_topography_test.yml + --job_id plane_analytic_no_topography_test + artifact_paths: "plane_analytic_no_topography_test/output_active/*" + + - label: ":computer: Analytic No Topography Test (2D, Float32)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_no_topography_float32_test.yml + --job_id plane_analytic_no_topography_float32_test + artifact_paths: "plane_analytic_no_topography_float32_test/output_active/*" + + - label: ":computer: Analytic Cosine Hills Test (2D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_cosine_hills_test.yml + --job_id plane_analytic_cosine_hills_test + artifact_paths: "plane_analytic_cosine_hills_test/output_active/*" + + - label: ":computer: Analytic Cosine Hills Test (2D, Float32)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_cosine_hills_float32_test.yml + --job_id plane_analytic_cosine_hills_float32_test + artifact_paths: "plane_analytic_cosine_hills_float32_test/output_active/*" - group: "Conservation check" steps: @@ -868,6 +895,391 @@ steps: slurm_gpus: 1 slurm_mem: 20G + - label: "GPU: Analytic No Topography Test (2D, Long Duration)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_no_topography_long_test.yml + --job_id gpu_plane_analytic_no_topography_long_test + artifact_paths: "gpu_plane_analytic_no_topography_long_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic No Topography Test (2D, Long Duration, Float32)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_no_topography_long_float32_test.yml + --job_id gpu_plane_analytic_no_topography_long_float32_test + artifact_paths: "gpu_plane_analytic_no_topography_long_float32_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic No Topography Test (2D, Long Duration, No Hyperdiff)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_no_topography_no_hyperdiff_long_test.yml + --job_id gpu_plane_analytic_no_topography_no_hyperdiff_long_test + artifact_paths: "gpu_plane_analytic_no_topography_no_hyperdiff_long_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic No Topography Test (2D, Long Duration, No Hyperdiff, Float32)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_no_topography_no_hyperdiff_long_float32_test.yml + --job_id gpu_plane_analytic_no_topography_no_hyperdiff_long_float32_test + artifact_paths: "gpu_plane_analytic_no_topography_no_hyperdiff_long_float32_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic No Topography Test (2D, Long Duration, Discrete Balance)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_no_topography_discrete_balance_long_test.yml + --job_id gpu_plane_analytic_no_topography_discrete_balance_long_test + artifact_paths: "gpu_plane_analytic_no_topography_discrete_balance_long_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic No Topography Test (2D, Long Duration, Discrete Balance, Float32)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_no_topography_discrete_balance_long_float32_test.yml + --job_id gpu_plane_analytic_no_topography_discrete_balance_long_float32_test + artifact_paths: "gpu_plane_analytic_no_topography_discrete_balance_long_float32_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic No Topography Test (Extruded 2D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/extruded_plane_analytic_no_topography_test.yml + --job_id gpu_extruded_plane_analytic_no_topography_test + artifact_paths: "gpu_extruded_plane_analytic_no_topography_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic No Topography Test (3D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/box_analytic_no_topography_test.yml + --job_id gpu_box_analytic_no_topography_test + artifact_paths: "gpu_box_analytic_no_topography_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic No Topography Test (3D, Float32)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/box_analytic_no_topography_float32_test.yml + --job_id gpu_box_analytic_no_topography_float32_test + artifact_paths: "gpu_box_analytic_no_topography_float32_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Cosine Hills Test (2D, Long Duration)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_cosine_hills_long_test.yml + --job_id gpu_plane_analytic_cosine_hills_long_test + artifact_paths: "gpu_plane_analytic_cosine_hills_long_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Cosine Hills Test (2D, Long Duration, Float32)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_cosine_hills_long_float32_test.yml + --job_id gpu_plane_analytic_cosine_hills_long_float32_test + artifact_paths: "gpu_plane_analytic_cosine_hills_long_float32_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Cosine Hills Test (2D, Long Duration, Stronger Sponge)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_cosine_hills_strong_sponge_long_test.yml + --job_id gpu_plane_analytic_cosine_hills_strong_sponge_long_test + artifact_paths: "gpu_plane_analytic_cosine_hills_strong_sponge_long_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Cosine Hills Test (2D, Long Duration, Weaker Sponge)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_cosine_hills_weak_sponge_long_test.yml + --job_id gpu_plane_analytic_cosine_hills_weak_sponge_long_test + artifact_paths: "gpu_plane_analytic_cosine_hills_weak_sponge_long_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Cosine Hills Test (2D, Long Duration, Higher Domain Top)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_cosine_hills_high_top_long_test.yml + --job_id gpu_plane_analytic_cosine_hills_high_top_long_test + artifact_paths: "gpu_plane_analytic_cosine_hills_high_top_long_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Cosine Hills Test (2D, Long Duration, Higher Domain Top, Higher Sponge)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_cosine_hills_high_top_high_sponge_long_test.yml + --job_id gpu_plane_analytic_cosine_hills_high_top_high_sponge_long_test + artifact_paths: "gpu_plane_analytic_cosine_hills_high_top_high_sponge_long_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Cosine Hills Test (Extruded 2D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/extruded_plane_analytic_cosine_hills_test.yml + --job_id gpu_extruded_plane_analytic_cosine_hills_test + artifact_paths: "gpu_extruded_plane_analytic_cosine_hills_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Cosine Hills Test (3D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/box_analytic_cosine_hills_test.yml + --job_id gpu_box_analytic_cosine_hills_test + artifact_paths: "gpu_box_analytic_cosine_hills_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Cosine Hills Test (3D, Float32)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/box_analytic_cosine_hills_float32_test.yml + --job_id gpu_box_analytic_cosine_hills_float32_test + artifact_paths: "gpu_box_analytic_cosine_hills_float32_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Agnesi Mountain Test (2D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_agnesi_mountain_test.yml + --job_id gpu_plane_analytic_agnesi_mountain_test + artifact_paths: "gpu_plane_analytic_agnesi_mountain_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Huge Schar Mountain Test (2D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_huge_schar_mountain_test.yml + --job_id gpu_plane_analytic_huge_schar_mountain_test + artifact_paths: "gpu_plane_analytic_huge_schar_mountain_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Big Schar Mountain Test (2D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_big_schar_mountain_test.yml + --job_id gpu_plane_analytic_big_schar_mountain_test + artifact_paths: "gpu_plane_analytic_big_schar_mountain_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Schar Mountain Test (2D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_schar_mountain_test.yml + --job_id gpu_plane_analytic_schar_mountain_test + artifact_paths: "gpu_plane_analytic_schar_mountain_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Small Schar Mountain Test (2D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_small_schar_mountain_test.yml + --job_id gpu_plane_analytic_small_schar_mountain_test + artifact_paths: "gpu_plane_analytic_small_schar_mountain_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Tiny Schar Mountain Test (2D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_tiny_schar_mountain_test.yml + --job_id gpu_plane_analytic_tiny_schar_mountain_test + artifact_paths: "gpu_plane_analytic_tiny_schar_mountain_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Teeny Tiny Schar Mountain Test (2D)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_teeny_tiny_schar_mountain_test.yml + --job_id gpu_plane_analytic_teeny_tiny_schar_mountain_test + artifact_paths: "gpu_plane_analytic_teeny_tiny_schar_mountain_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Schar Mountain Test (2D, Float32)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_schar_mountain_float32_test.yml + --job_id gpu_plane_analytic_schar_mountain_float32_test + artifact_paths: "gpu_plane_analytic_schar_mountain_float32_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Tiny Schar Mountain Test (2D, Float32)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_tiny_schar_mountain_float32_test.yml + --job_id gpu_plane_analytic_tiny_schar_mountain_float32_test + artifact_paths: "gpu_plane_analytic_tiny_schar_mountain_float32_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Schar Mountain Test (2D, No Hyperdiff)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_schar_mountain_no_hyperdiff_test.yml + --job_id gpu_plane_analytic_schar_mountain_no_hyperdiff_test + artifact_paths: "gpu_plane_analytic_schar_mountain_no_hyperdiff_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Schar Mountain Test (2D, No Hyperdiff, Float32)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_schar_mountain_no_hyperdiff_float32_test.yml + --job_id gpu_plane_analytic_schar_mountain_no_hyperdiff_float32_test + artifact_paths: "gpu_plane_analytic_schar_mountain_no_hyperdiff_float32_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Schar Mountain Test (2D, Stretched Grid)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_schar_mountain_stretched_grid_test.yml + --job_id gpu_plane_analytic_schar_mountain_stretched_grid_test + artifact_paths: "gpu_plane_analytic_schar_mountain_stretched_grid_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Schar Mountain Test (2D, Higher Horizontal Resolution)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_schar_mountain_high_horz_res_test.yml + --job_id gpu_plane_analytic_schar_mountain_high_horz_res_test + artifact_paths: "gpu_plane_analytic_schar_mountain_high_horz_res_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Schar Mountain Test (2D, Higher Vertical Resolution)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_schar_mountain_high_vert_res_test.yml + --job_id gpu_plane_analytic_schar_mountain_high_vert_res_test + artifact_paths: "gpu_plane_analytic_schar_mountain_high_vert_res_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Schar Mountain Test (2D, Higher Domain Top)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_schar_mountain_high_top_test.yml + --job_id gpu_plane_analytic_schar_mountain_high_top_test + artifact_paths: "gpu_plane_analytic_schar_mountain_high_top_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Schar Mountain Test (2D, Higher Domain Top, Higher Sponge)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_schar_mountain_high_top_high_sponge_test.yml + --job_id gpu_plane_analytic_schar_mountain_high_top_high_sponge_test + artifact_paths: "gpu_plane_analytic_schar_mountain_high_top_high_sponge_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + + - label: "GPU: Analytic Schar Mountain Test (2D, Higher Domain Top, Weaker and Higher Sponge)" + command: > + julia --color=yes --project=examples examples/hybrid/driver.jl + --config_file $CONFIG_PATH/plane_analytic_schar_mountain_high_top_weak_high_sponge_test.yml + --job_id gpu_plane_analytic_schar_mountain_high_top_weak_high_sponge_test + artifact_paths: "gpu_plane_analytic_schar_mountain_high_top_weak_high_sponge_test/output_active/*" + env: + CLIMACOMMS_DEVICE: "CUDA" + agents: + slurm_gpus: 1 + - group: "Benchmarks" steps: diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml index 34d4f1684c7..82372187696 100644 --- a/config/default_configs/default_config.yml +++ b/config/default_configs/default_config.yml @@ -55,6 +55,9 @@ topo_smoothing: smoothing_order: help: "Define the surface smoothing kernel factor (integer) [`3 (default)`]" value: 3 +max_topography_height: + help: "The maximum height in meters for certain choices of `topography` [`25.0 (default)`]" + value: 25.0 # ODE use_newton_rtol: help: "Whether to check if the current iteration of Newton's method has an error within a relative tolerance, instead of always taking the maximum number of iterations (only for ClimaTimeSteppers.jl)" @@ -154,6 +157,9 @@ start_date: forcing: help: "Forcing [`nothing` (default), `held_suarez`]" value: ~ +analytic_check: + help: "Test final state against analytic reference, if one is available [`false` (default), `true`]" + value: false test_dycore_consistency: help: "Test dycore consistency [`false` (default), `true`]" value: false diff --git a/config/model_configs/box_analytic_cosine_hills_float32_test.yml b/config/model_configs/box_analytic_cosine_hills_float32_test.yml new file mode 100644 index 00000000000..0e1674ad6c6 --- /dev/null +++ b/config/model_configs/box_analytic_cosine_hills_float32_test.yml @@ -0,0 +1,20 @@ +config: "box" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine3D" +x_max: 300e3 +y_max: 300e3 +z_max: 21e3 +x_elem: 15 +y_elem: 15 +z_elem: 50 +z_stretch: false +dt: "12secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/box_analytic_cosine_hills_test.yml b/config/model_configs/box_analytic_cosine_hills_test.yml new file mode 100644 index 00000000000..b4b70a1f255 --- /dev/null +++ b/config/model_configs/box_analytic_cosine_hills_test.yml @@ -0,0 +1,20 @@ +config: "box" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine3D" +x_max: 300e3 +y_max: 300e3 +z_max: 21e3 +x_elem: 15 +y_elem: 15 +z_elem: 50 +z_stretch: false +dt: "12secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/box_analytic_no_topography_float32_test.yml b/config/model_configs/box_analytic_no_topography_float32_test.yml new file mode 100644 index 00000000000..a14330d2045 --- /dev/null +++ b/config/model_configs/box_analytic_no_topography_float32_test.yml @@ -0,0 +1,21 @@ +config: "box" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine3D" +x_max: 300e3 +y_max: 300e3 +z_max: 21e3 +x_elem: 15 +y_elem: 15 +z_elem: 50 +z_stretch: false +max_topography_height: 0 +dt: "12secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/box_analytic_no_topography_test.yml b/config/model_configs/box_analytic_no_topography_test.yml new file mode 100644 index 00000000000..ee50c33696b --- /dev/null +++ b/config/model_configs/box_analytic_no_topography_test.yml @@ -0,0 +1,21 @@ +config: "box" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine3D" +x_max: 300e3 +y_max: 300e3 +z_max: 21e3 +x_elem: 15 +y_elem: 15 +z_elem: 50 +z_stretch: false +max_topography_height: 0 +dt: "12secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/extruded_plane_analytic_cosine_hills_test.yml b/config/model_configs/extruded_plane_analytic_cosine_hills_test.yml new file mode 100644 index 00000000000..abd281236b6 --- /dev/null +++ b/config/model_configs/extruded_plane_analytic_cosine_hills_test.yml @@ -0,0 +1,20 @@ +config: "box" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +y_max: 40e3 +z_max: 21e3 +x_elem: 15 +y_elem: 2 +z_elem: 50 +z_stretch: false +dt: "12secs" +t_end: "5days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/extruded_plane_analytic_no_topography_test.yml b/config/model_configs/extruded_plane_analytic_no_topography_test.yml new file mode 100644 index 00000000000..54d65515adc --- /dev/null +++ b/config/model_configs/extruded_plane_analytic_no_topography_test.yml @@ -0,0 +1,21 @@ +config: "box" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +y_max: 40e3 +z_max: 21e3 +x_elem: 15 +y_elem: 2 +z_elem: 50 +z_stretch: false +max_topography_height: 0 +dt: "12secs" +t_end: "5days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_agnesi_mountain_test.yml b/config/model_configs/plane_analytic_agnesi_mountain_test.yml new file mode 100644 index 00000000000..8f8510f91ea --- /dev/null +++ b/config/model_configs/plane_analytic_agnesi_mountain_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Agnesi" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +z_stretch: false +dt: "0.5secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_big_schar_mountain_test.yml b/config/model_configs/plane_analytic_big_schar_mountain_test.yml new file mode 100644 index 00000000000..9422d3e932a --- /dev/null +++ b/config/model_configs/plane_analytic_big_schar_mountain_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +z_stretch: false +max_topography_height: 250 +dt: "0.5secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_cosine_hills_float32_test.yml b/config/model_configs/plane_analytic_cosine_hills_float32_test.yml new file mode 100644 index 00000000000..8599300b6e1 --- /dev/null +++ b/config/model_configs/plane_analytic_cosine_hills_float32_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +dt: "12secs" +t_end: "10days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_cosine_hills_high_top_high_sponge_long_test.yml b/config/model_configs/plane_analytic_cosine_hills_high_top_high_sponge_long_test.yml new file mode 100644 index 00000000000..fb3980528f6 --- /dev/null +++ b/config/model_configs/plane_analytic_cosine_hills_high_top_high_sponge_long_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 84e3 +x_elem: 15 +z_elem: 200 +z_stretch: false +dt: "12secs" +t_end: "100days" +rayleigh_sponge: true +toml: [toml/analytic_topography_high_sponge_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_cosine_hills_high_top_long_test.yml b/config/model_configs/plane_analytic_cosine_hills_high_top_long_test.yml new file mode 100644 index 00000000000..cb6eea6bf56 --- /dev/null +++ b/config/model_configs/plane_analytic_cosine_hills_high_top_long_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 84e3 +x_elem: 15 +z_elem: 200 +z_stretch: false +dt: "12secs" +t_end: "100days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_cosine_hills_long_float32_test.yml b/config/model_configs/plane_analytic_cosine_hills_long_float32_test.yml new file mode 100644 index 00000000000..9af9b6aa0ac --- /dev/null +++ b/config/model_configs/plane_analytic_cosine_hills_long_float32_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +dt: "12secs" +t_end: "100days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_cosine_hills_long_test.yml b/config/model_configs/plane_analytic_cosine_hills_long_test.yml new file mode 100644 index 00000000000..32a9187e204 --- /dev/null +++ b/config/model_configs/plane_analytic_cosine_hills_long_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +dt: "12secs" +t_end: "100days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_cosine_hills_strong_sponge_long_test.yml b/config/model_configs/plane_analytic_cosine_hills_strong_sponge_long_test.yml new file mode 100644 index 00000000000..d02a180941a --- /dev/null +++ b/config/model_configs/plane_analytic_cosine_hills_strong_sponge_long_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +dt: "12secs" +t_end: "100days" +rayleigh_sponge: true +toml: [toml/analytic_topography_strong_sponge_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_cosine_hills_test.yml b/config/model_configs/plane_analytic_cosine_hills_test.yml new file mode 100644 index 00000000000..0586af9c056 --- /dev/null +++ b/config/model_configs/plane_analytic_cosine_hills_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +dt: "12secs" +t_end: "10days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_cosine_hills_weak_sponge_long_test.yml b/config/model_configs/plane_analytic_cosine_hills_weak_sponge_long_test.yml new file mode 100644 index 00000000000..ea9dea67be4 --- /dev/null +++ b/config/model_configs/plane_analytic_cosine_hills_weak_sponge_long_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +dt: "12secs" +t_end: "100days" +rayleigh_sponge: true +toml: [toml/analytic_topography_weak_sponge_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_huge_schar_mountain_test.yml b/config/model_configs/plane_analytic_huge_schar_mountain_test.yml new file mode 100644 index 00000000000..8eddab2c68e --- /dev/null +++ b/config/model_configs/plane_analytic_huge_schar_mountain_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +z_stretch: false +max_topography_height: 2500 +dt: "0.05secs" +t_end: "10hours" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_no_topography_discrete_balance_long_float32_test.yml b/config/model_configs/plane_analytic_no_topography_discrete_balance_long_float32_test.yml new file mode 100644 index 00000000000..df23e353a3f --- /dev/null +++ b/config/model_configs/plane_analytic_no_topography_discrete_balance_long_float32_test.yml @@ -0,0 +1,20 @@ +config: "plane" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +discrete_hydrostatic_balance: true +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +max_topography_height: 0 +dt: "12secs" +t_end: "100days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_no_topography_discrete_balance_long_test.yml b/config/model_configs/plane_analytic_no_topography_discrete_balance_long_test.yml new file mode 100644 index 00000000000..df23e353a3f --- /dev/null +++ b/config/model_configs/plane_analytic_no_topography_discrete_balance_long_test.yml @@ -0,0 +1,20 @@ +config: "plane" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +discrete_hydrostatic_balance: true +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +max_topography_height: 0 +dt: "12secs" +t_end: "100days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_no_topography_float32_test.yml b/config/model_configs/plane_analytic_no_topography_float32_test.yml new file mode 100644 index 00000000000..e1c1950dbff --- /dev/null +++ b/config/model_configs/plane_analytic_no_topography_float32_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +max_topography_height: 0 +dt: "12secs" +t_end: "10days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_no_topography_long_float32_test.yml b/config/model_configs/plane_analytic_no_topography_long_float32_test.yml new file mode 100644 index 00000000000..e7375d09d9c --- /dev/null +++ b/config/model_configs/plane_analytic_no_topography_long_float32_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +max_topography_height: 0 +dt: "12secs" +t_end: "100days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_no_topography_long_test.yml b/config/model_configs/plane_analytic_no_topography_long_test.yml new file mode 100644 index 00000000000..12a2767b380 --- /dev/null +++ b/config/model_configs/plane_analytic_no_topography_long_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +max_topography_height: 0 +dt: "12secs" +t_end: "100days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_no_topography_no_hyperdiff_long_float32_test.yml b/config/model_configs/plane_analytic_no_topography_no_hyperdiff_long_float32_test.yml new file mode 100644 index 00000000000..73c153e3970 --- /dev/null +++ b/config/model_configs/plane_analytic_no_topography_no_hyperdiff_long_float32_test.yml @@ -0,0 +1,20 @@ +config: "plane" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +max_topography_height: 0 +dt: "12secs" +t_end: "100days" +hyperdiff: false +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_no_topography_no_hyperdiff_long_test.yml b/config/model_configs/plane_analytic_no_topography_no_hyperdiff_long_test.yml new file mode 100644 index 00000000000..7b07f63ebf7 --- /dev/null +++ b/config/model_configs/plane_analytic_no_topography_no_hyperdiff_long_test.yml @@ -0,0 +1,20 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +max_topography_height: 0 +dt: "12secs" +t_end: "100days" +hyperdiff: false +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_no_topography_test.yml b/config/model_configs/plane_analytic_no_topography_test.yml new file mode 100644 index 00000000000..cdb4e41960b --- /dev/null +++ b/config/model_configs/plane_analytic_no_topography_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Cosine2D" +x_max: 300e3 +z_max: 21e3 +x_elem: 15 +z_elem: 50 +z_stretch: false +max_topography_height: 0 +dt: "12secs" +t_end: "10days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_schar_mountain_float32_test.yml b/config/model_configs/plane_analytic_schar_mountain_float32_test.yml new file mode 100644 index 00000000000..0a87fc1c045 --- /dev/null +++ b/config/model_configs/plane_analytic_schar_mountain_float32_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +z_stretch: false +dt: "0.5secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_schar_mountain_high_horz_res_test.yml b/config/model_configs/plane_analytic_schar_mountain_high_horz_res_test.yml new file mode 100644 index 00000000000..6e870079117 --- /dev/null +++ b/config/model_configs/plane_analytic_schar_mountain_high_horz_res_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 200 +z_elem: 100 +z_stretch: false +dt: "0.25secs" +t_end: "1.5days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_schar_mountain_high_top_high_sponge_test.yml b/config/model_configs/plane_analytic_schar_mountain_high_top_high_sponge_test.yml new file mode 100644 index 00000000000..5e4742ad785 --- /dev/null +++ b/config/model_configs/plane_analytic_schar_mountain_high_top_high_sponge_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 84e3 +x_elem: 100 +z_elem: 400 +z_stretch: false +dt: "0.5secs" +t_end: "1days" +rayleigh_sponge: true +toml: [toml/analytic_topography_high_sponge_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_schar_mountain_high_top_test.yml b/config/model_configs/plane_analytic_schar_mountain_high_top_test.yml new file mode 100644 index 00000000000..6ef2826f9c5 --- /dev/null +++ b/config/model_configs/plane_analytic_schar_mountain_high_top_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 84e3 +x_elem: 100 +z_elem: 400 +z_stretch: false +dt: "0.5secs" +t_end: "1days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_schar_mountain_high_top_weak_high_sponge_test.yml b/config/model_configs/plane_analytic_schar_mountain_high_top_weak_high_sponge_test.yml new file mode 100644 index 00000000000..9e5b53c8580 --- /dev/null +++ b/config/model_configs/plane_analytic_schar_mountain_high_top_weak_high_sponge_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 84e3 +x_elem: 100 +z_elem: 400 +z_stretch: false +dt: "0.5secs" +t_end: "1days" +rayleigh_sponge: true +toml: [toml/analytic_topography_weak_high_sponge_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_schar_mountain_high_vert_res_test.yml b/config/model_configs/plane_analytic_schar_mountain_high_vert_res_test.yml new file mode 100644 index 00000000000..6a34e6b3fb3 --- /dev/null +++ b/config/model_configs/plane_analytic_schar_mountain_high_vert_res_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 400 +z_stretch: false +dt: "0.5secs" +t_end: "1.5days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_schar_mountain_no_hyperdiff_float32_test.yml b/config/model_configs/plane_analytic_schar_mountain_no_hyperdiff_float32_test.yml new file mode 100644 index 00000000000..f8918cf2c6d --- /dev/null +++ b/config/model_configs/plane_analytic_schar_mountain_no_hyperdiff_float32_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +z_stretch: false +dt: "0.2secs" +t_end: "2days" +hyperdiff: false +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_schar_mountain_no_hyperdiff_test.yml b/config/model_configs/plane_analytic_schar_mountain_no_hyperdiff_test.yml new file mode 100644 index 00000000000..9ffdad17f84 --- /dev/null +++ b/config/model_configs/plane_analytic_schar_mountain_no_hyperdiff_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +z_stretch: false +dt: "0.2secs" +t_end: "2days" +hyperdiff: false +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_schar_mountain_stretched_grid_test.yml b/config/model_configs/plane_analytic_schar_mountain_stretched_grid_test.yml new file mode 100644 index 00000000000..287bac31b60 --- /dev/null +++ b/config/model_configs/plane_analytic_schar_mountain_stretched_grid_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +dz_bottom: 50 +dz_top: 1000 +dt: "0.5secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_schar_mountain_test.yml b/config/model_configs/plane_analytic_schar_mountain_test.yml new file mode 100644 index 00000000000..c3b5f2138fc --- /dev/null +++ b/config/model_configs/plane_analytic_schar_mountain_test.yml @@ -0,0 +1,18 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +z_stretch: false +dt: "0.5secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_small_schar_mountain_test.yml b/config/model_configs/plane_analytic_small_schar_mountain_test.yml new file mode 100644 index 00000000000..1f8e7ea5611 --- /dev/null +++ b/config/model_configs/plane_analytic_small_schar_mountain_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +z_stretch: false +max_topography_height: 2.5 +dt: "0.5secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_teeny_tiny_schar_mountain_test.yml b/config/model_configs/plane_analytic_teeny_tiny_schar_mountain_test.yml new file mode 100644 index 00000000000..b8bf2432ed8 --- /dev/null +++ b/config/model_configs/plane_analytic_teeny_tiny_schar_mountain_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +z_stretch: false +max_topography_height: 0.00000025 +dt: "0.5secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_tiny_schar_mountain_float32_test.yml b/config/model_configs/plane_analytic_tiny_schar_mountain_float32_test.yml new file mode 100644 index 00000000000..172b18849df --- /dev/null +++ b/config/model_configs/plane_analytic_tiny_schar_mountain_float32_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float32" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +z_stretch: false +max_topography_height: 0.25 +dt: "0.5secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/config/model_configs/plane_analytic_tiny_schar_mountain_test.yml b/config/model_configs/plane_analytic_tiny_schar_mountain_test.yml new file mode 100644 index 00000000000..21f94015829 --- /dev/null +++ b/config/model_configs/plane_analytic_tiny_schar_mountain_test.yml @@ -0,0 +1,19 @@ +config: "plane" +FLOAT_TYPE: "Float64" +initial_condition: "ConstantBuoyancyFrequencyProfile" +topography: "Schar" +x_max: 100e3 +z_max: 21e3 +x_elem: 100 +z_elem: 100 +z_stretch: false +max_topography_height: 0.25 +dt: "0.5secs" +t_end: "4days" +rayleigh_sponge: true +toml: [toml/analytic_topography_test.toml] +analytic_check: true +output_default_diagnostics: false +diagnostics: + - short_name: [orog, ua, wa, uapredicted, wapredicted, uaerror, waerror, c1error, c3error] + period: 1hours diff --git a/examples/hybrid/disable_topography_dss.jl b/examples/hybrid/disable_topography_dss.jl new file mode 100644 index 00000000000..2457c2df863 --- /dev/null +++ b/examples/hybrid/disable_topography_dss.jl @@ -0,0 +1,166 @@ +import ClimaCore: + DataLayouts, + Geometry, + Topologies, + Grids, + Hypsography, + Spaces, + Fields, + Operators +import ClimaTimeSteppers +import ClimaComms +import LinearAlgebra: adjoint, ldiv!, DenseMatrix, lu, norm + +ClimaTimeSteppers.NVTX.@annotate function ClimaTimeSteppers.solve_newton!( + alg::ClimaTimeSteppers.NewtonsMethod, + cache, + x, + f!, + j! = nothing, + pre_iteration! = nothing, + post_solve! = nothing, +) + (; max_iters, update_j, krylov_method, convergence_checker, verbose) = alg + (; krylov_method_cache, convergence_checker_cache) = cache + (; Δx, f, j) = cache + if (!isnothing(j)) && ClimaTimeSteppers.needs_update!( + update_j, + ClimaTimeSteppers.NewNewtonSolve(), + ) + j!(j, x) + end + for n in 1:max_iters + # Compute Δx[n]. + if (!isnothing(j)) && ClimaTimeSteppers.needs_update!( + update_j, + ClimaTimeSteppers.NewNewtonIteration(), + ) + j!(j, x) + end + f!(f, x) + if isnothing(krylov_method) + if j isa DenseMatrix + ldiv!(Δx, lu(j), f) # Highly inefficient! Only used for testing. + else + ldiv!(Δx, j, f) + end + else + ClimaTimeSteppers.solve_krylov!( + krylov_method, + krylov_method_cache, + Δx, + x, + f!, + f, + n, + pre_iteration!, + j, + ) + end + ClimaTimeSteppers.is_verbose(verbose) && + @info "Newton iteration $n: ‖x‖ = $(norm(x)), ‖Δx‖ = $(norm(Δx))" + + c_dss_buffer = Spaces.create_dss_buffer(Δx.c) + f_dss_buffer = Spaces.create_dss_buffer(Δx.f) + Spaces.weighted_dss!(Δx.c => c_dss_buffer, Δx.f => f_dss_buffer) + + x .-= Δx + # Update x[n] with Δx[n - 1], and exit the loop if Δx[n] is not needed. + # Check for convergence if necessary. + if ClimaTimeSteppers.is_converged!( + convergence_checker, + convergence_checker_cache, + x, + Δx, + n, + ) + isnothing(post_solve!) || post_solve!(x) + break + elseif n == max_iters + isnothing(post_solve!) || post_solve!(x) + else + isnothing(pre_iteration!) || pre_iteration!(x) + end + if ClimaTimeSteppers.is_verbose(verbose) && n == max_iters + @warn "Newton's method did not converge within $n iterations: ‖x‖ = $(norm(x)), ‖Δx‖ = $(norm(Δx))" + end + end +end + +function Grids._ExtrudedFiniteDifferenceGrid( + horizontal_grid::Grids.AbstractGrid, + vertical_grid::Grids.FiniteDifferenceGrid, + adaption::Grids.HypsographyAdaption, + global_geometry::Geometry.AbstractGlobalGeometry, + face_z::DataLayouts.AbstractData{Geometry.ZPoint{FT}}, +) where {FT} + # construct the "flat" grid + # avoid cached constructor so that it gets cleaned up automatically + flat_grid = Grids._ExtrudedFiniteDifferenceGrid( + horizontal_grid, + vertical_grid, + Grids.Flat(), + global_geometry, + ) + center_flat_space = Spaces.space(flat_grid, Grids.CellCenter()) + face_flat_space = Spaces.space(flat_grid, Grids.CellFace()) + + # compute the "z-only local geometry" based on face z coords + ArrayType = ClimaComms.array_type(horizontal_grid.topology) + # currently only works on Arrays + (center_z_local_geometry, face_z_local_geometry) = Grids.fd_geometry_data( + Adapt.adapt(Array, face_z), + Val(Topologies.isperiodic(vertical_grid.topology)), + ) + + center_z_local_geometry = Adapt.adapt(ArrayType, center_z_local_geometry) + face_z_local_geometry = Adapt.adapt(ArrayType, face_z_local_geometry) + + # compute ∇Z at face and centers + grad = Operators.Gradient() + + center_∇Z_field = + grad.( + Fields.Field( + center_z_local_geometry, + center_flat_space, + ).coordinates.z + ) + face_∇Z_field = + grad.( + Fields.Field(face_z_local_geometry, face_flat_space).coordinates.z + ) + # buffer = (; + # c = Spaces.create_dss_buffer(center_∇Z_field), + # f = Spaces.create_dss_buffer(face_∇Z_field), + # ) + + # Spaces.weighted_dss!(center_∇Z_field => buffer.c, face_∇Z_field => buffer.f) + + # construct full local geometry + center_local_geometry = + Geometry.product_geometry.( + horizontal_grid.local_geometry, + center_z_local_geometry, + Ref(global_geometry), + Ref(Geometry.WVector(1)) .* + adjoint.(Fields.field_values(center_∇Z_field)), + ) + face_local_geometry = + Geometry.product_geometry.( + horizontal_grid.local_geometry, + face_z_local_geometry, + Ref(global_geometry), + Ref(Geometry.WVector(1)) .* + adjoint.(Fields.field_values(face_∇Z_field)), + ) + + return Grids.ExtrudedFiniteDifferenceGrid( + horizontal_grid, + vertical_grid, + adaption, + global_geometry, + center_local_geometry, + face_local_geometry, + ) +end diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl index 75ba3ffe2a5..81ba7a913d4 100644 --- a/examples/hybrid/driver.jl +++ b/examples/hybrid/driver.jl @@ -10,6 +10,10 @@ import ClimaAtmos as CA import Random Random.seed!(1234) +include("fix_asymmetric_metric_tensor.jl") +include("fix_1d_spectral_space_on_gpu.jl") +include("disable_topography_dss.jl") + if !(@isdefined config) (; config_file, job_id) = CA.commandline_kwargs() config = CA.AtmosConfig(config_file; job_id) @@ -22,9 +26,10 @@ sol_res = CA.solve_atmos!(simulation) (; p) = integrator import ClimaCore -import ClimaCore: Topologies, Quadratures, Spaces +import ClimaCore: Topologies, Quadratures, Spaces, Fields import ClimaAtmos.InitialConditions as ICs using Statistics: mean +import LinearAlgebra: norm_sqr import ClimaAtmos.Parameters as CAP import Thermodynamics as TD import ClimaComms @@ -114,6 +119,39 @@ end @info "Callback verification, n_expected_calls: $(CA.n_expected_calls(integrator))" @info "Callback verification, n_measured_calls: $(CA.n_measured_calls(integrator))" +if config.parsed_args["analytic_check"] + @info "Comparing final state against predicted steady-state solution" + + Y_end = integrator.sol.u[end] + (; predicted_steady_state, params) = integrator.p + @assert !isnothing(predicted_steady_state) + + FT = eltype(Y_end) + (; zd_rayleigh) = params + + ᶜuₕ_error_squared = + norm_sqr.(Y_end.c.uₕ .- CA.C12.(predicted_steady_state.ᶜu)) + ᶠu₃_error_squared = + norm_sqr.(Y_end.f.u₃ .- CA.C3.(predicted_steady_state.ᶠu)) + + ᶜsponge_mask = FT.(Fields.coordinate_field(Y_end.c).z .< zd_rayleigh) + ᶠsponge_mask = FT.(Fields.coordinate_field(Y_end.f).z .< zd_rayleigh) + uₕ_rmse = sqrt(sum(ᶜuₕ_error_squared .* ᶜsponge_mask) / sum(ᶜsponge_mask)) + u₃_rmse = sqrt(sum(ᶠu₃_error_squared .* ᶠsponge_mask) / sum(ᶠsponge_mask)) + + uₕ_rmse_by_level = map(1:5) do level + sqrt(mean(Fields.level(ᶜuₕ_error_squared, level))) + end + u₃_rmse_by_level = map(1:5) do level + sqrt(mean(Fields.level(ᶠu₃_error_squared, level - Fields.half))) + end + + @info " RMSE of uₕ below sponge layer: $uₕ_rmse" + @info " RMSE of u₃ below sponge layer: $u₃_rmse" + @info " RMSE of uₕ on first 5 levels: $uₕ_rmse_by_level" + @info " RMSE of u₃ on first 5 levels: $u₃_rmse_by_level" +end + # Conservation checks if config.parsed_args["check_conservation"] FT = Spaces.undertype(axes(sol.u[end].c.ρ)) diff --git a/examples/hybrid/fix_1d_spectral_space_on_gpu.jl b/examples/hybrid/fix_1d_spectral_space_on_gpu.jl new file mode 100644 index 00000000000..c1364bfc4b6 --- /dev/null +++ b/examples/hybrid/fix_1d_spectral_space_on_gpu.jl @@ -0,0 +1,1090 @@ +import ClimaCore: + ClimaCore, + Utilities, + RecursiveApply, + DataLayouts, + Geometry, + Meshes, + Topologies, + Quadratures, + Grids, + Spaces, + Fields, + Operators, + Hypsography, + Remapping, + slab +import ClimaCore.RecursiveApply: ⊞, ⊟, ⊠ +import ClimaComms +import CUDA +import Adapt +import StaticArrays: SMatrix + +ClimaCoreCUDAExt = Base.get_extension(ClimaCore, :ClimaCoreCUDAExt) + +########## +## ClimaCore.jl/ext/cuda/data_layouts.jl +########## + +function Base.copyto!( + dest::DataLayouts.IFH{S, Ni}, + bc::DataLayouts.BroadcastedUnionIFH{S, Ni}, + ::ClimaCoreCUDAExt.ToCUDA, +) where {S, Ni} + _, _, _, _, Nh = size(bc) + if Nh > 0 + ClimaCoreCUDAExt.auto_launch!( + ClimaCoreCUDAExt.knl_copyto!, + (dest, bc), + dest; + threads_s = (Ni, 1), + blocks_s = (Nh, 1), + ) + end + return dest +end + +function Base.copyto!( + dest::DataLayouts.VIFH{S, Nv, Ni}, + bc::DataLayouts.BroadcastedUnionVIFH{S, Nv, Ni}, + ::ClimaCoreCUDAExt.ToCUDA, +) where {S, Nv, Ni} + _, _, _, _, Nh = size(bc) + if Nv > 0 && Nh > 0 + Nv_per_block = min(Nv, fld(256, Ni)) + Nv_blocks = cld(Nv, Nv_per_block) + ClimaCoreCUDAExt.auto_launch!( + ClimaCoreCUDAExt.knl_copyto!, + (dest, bc), + dest; + threads_s = (Ni, 1, Nv_per_block), + blocks_s = (Nh, Nv_blocks), + ) + end + return dest +end + +########## +## ClimaCore.jl/ext/cuda/fields.jl +########## + +function ClimaCoreCUDAExt.mapreduce_cuda( + f, + op, + field::Fields.Field{V}; + weighting = false, + opargs..., +) where { + S, + V <: Union{ + DataLayouts.VF{S}, + DataLayouts.IFH{S}, + DataLayouts.IJFH{S}, + DataLayouts.VIFH{S}, + DataLayouts.VIJFH{S}, + }, +} + data = Fields.field_values(field) + pdata = parent(data) + T = eltype(pdata) + (Ni, Nj, Nk, Nv, Nh) = size(data) + Nf = div(length(pdata), prod(size(data))) # length of field dimension + wt = Spaces.weighted_jacobian(axes(field)) + pwt = parent(wt) + + nitems = Nv * Ni * Nj * Nk * Nh + max_threads = 256# 512 1024 + nthreads = min(max_threads, nitems) + # perform n ops during loading to shmem (this is a tunable parameter) + n_ops_on_load = cld(nitems, nthreads) == 1 ? 0 : 7 + effective_blksize = nthreads * (n_ops_on_load + 1) + nblocks = cld(nitems, effective_blksize) + + reduce_cuda = CUDA.CuArray{T}(undef, nblocks, Nf) + shmemsize = nthreads + # place each field on a different block + CUDA.@cuda always_inline = true threads = (nthreads) blocks = (nblocks, Nf) ClimaCoreCUDAExt.mapreduce_cuda_kernel!( + reduce_cuda, + f, + op, + pdata, + pwt, + weighting, + n_ops_on_load, + Val(shmemsize), + nitems * Nf, + ClimaCoreCUDAExt._dataview, + ClimaCoreCUDAExt._get_gidx, + ClimaCoreCUDAExt._cuda_intrablock_reduce!, + ) + # reduce block data + if nblocks > 1 + nthreads = min(32, nblocks) + shmemsize = nthreads + CUDA.@cuda always_inline = true threads = (nthreads) blocks = (Nf) ClimaCoreCUDAExt.reduce_cuda_blocks_kernel!( + reduce_cuda, + op, + Val(shmemsize), + ) + end + return DataLayouts.DataF{S}(Array(Array(reduce_cuda)[1, :])) +end + +function ClimaCoreCUDAExt.mapreduce_cuda_kernel!( + reduce_cuda::AbstractArray{T, 2}, + f, + op, + pdata::AbstractArray{T, N}, + pwt::AbstractArray{T, N}, + weighting::Bool, + n_ops_on_load::Int, + ::Val{shmemsize}, + nitems, + _dataview, + _get_gidx, + _cuda_intrablock_reduce!, +) where {T, N, shmemsize} + blksize = CUDA.blockDim().x + nblk = CUDA.gridDim().x + tidx = CUDA.threadIdx().x + bidx = CUDA.blockIdx().x + fidx = CUDA.blockIdx().y + dataview = _dataview(pdata, fidx) + effective_blksize = blksize * (n_ops_on_load + 1) + gidx = _get_gidx(tidx, bidx, effective_blksize) + reduction = CUDA.CuStaticSharedArray(T, shmemsize) + reduction[tidx] = 0 + + # load shmem + if gidx ≤ nitems + if weighting + reduction[tidx] = f(dataview[gidx]) * pwt[gidx] + for n_ops in 1:n_ops_on_load + gidx2 = + _get_gidx(tidx + blksize * n_ops, bidx, effective_blksize) + if gidx2 ≤ nitems + reduction[tidx] = + op(reduction[tidx], f(dataview[gidx2]) * pwt[gidx2]) + end + end + else + reduction[tidx] = f(dataview[gidx]) + for n_ops in 1:n_ops_on_load + gidx2 = + _get_gidx(tidx + blksize * n_ops, bidx, effective_blksize) + if gidx2 ≤ nitems + reduction[tidx] = op(reduction[tidx], f(dataview[gidx2])) + end + end + end + end + CUDA.sync_threads() + _cuda_intrablock_reduce!(op, reduction, tidx, blksize) + + tidx == 1 && (reduce_cuda[bidx, fidx] = reduction[1]) + return nothing +end + +@inline ClimaCoreCUDAExt._dataview( + pdata::AbstractArray{FT, 3}, + fidx, +) where {FT} = view(pdata, :, fidx:fidx, :) + +########## +## ClimaCore.jl/ext/cuda/remapping_distributed.jl +########## + +function ClimaCoreCUDAExt.set_interpolated_values_kernel!( + out::AbstractArray, + (I,)::NTuple{1}, + local_horiz_indices, + vert_interpolation_weights, + vert_bounding_indices, + field_values, +) + # TODO: Check the memory access pattern. This was not optimized and likely inefficient! + num_horiz = length(local_horiz_indices) + num_vert = length(vert_bounding_indices) + num_fields = length(field_values) + + hindex = + (CUDA.blockIdx().x - Int32(1)) * CUDA.blockDim().x + CUDA.threadIdx().x + vindex = + (CUDA.blockIdx().y - Int32(1)) * CUDA.blockDim().y + CUDA.threadIdx().y + findex = + (CUDA.blockIdx().z - Int32(1)) * CUDA.blockDim().z + CUDA.threadIdx().z + + totalThreadsX = CUDA.gridDim().x * CUDA.blockDim().x + totalThreadsY = CUDA.gridDim().y * CUDA.blockDim().y + totalThreadsZ = CUDA.gridDim().z * CUDA.blockDim().z + + _, Nq = size(I) + + for i in hindex:totalThreadsX:num_horiz + h = local_horiz_indices[i] + for j in vindex:totalThreadsY:num_vert + v_lo, v_hi = vert_bounding_indices[j] + A, B = vert_interpolation_weights[j] + for k in findex:totalThreadsZ:num_fields + if i ≤ num_horiz && j ≤ num_vert && k ≤ num_fields + out[i, j, k] = 0 + for t in 1:Nq + out[i, j, k] += + I[i, t] * ( + A * + field_values[k][t, nothing, nothing, v_lo, h] + + B * + field_values[k][t, nothing, nothing, v_hi, h] + ) + end + end + end + end + end + return nothing +end + +########## +## ClimaCore.jl/src/Remapping/distributed_remapping.jl +########## + +function Remapping._set_interpolated_values!( + out::AbstractArray, + fields::AbstractArray{<:Fields.Field}, + _scratch_field_values, + local_horiz_indices, + local_horiz_interpolation_weights, + ::Nothing, + ::Nothing, +) + CUDA.@allowscalar Remapping._set_interpolated_values_device!( + out, + fields, + _scratch_field_values, + local_horiz_indices, + local_horiz_interpolation_weights, + nothing, + nothing, + ClimaComms.CPUSingleThreaded(), + ) # TODO: Fix the kernel function above and remove the @allowscalar. +end + +########## +## ClimaCore.jl/ext/cuda/operators_sem_shmem.jl +########## + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_shmem( + space, + ::Val{Nvt}, + op::Operators.Divergence{(1,)}, + arg, +) where {Nvt} + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + # allocate temp output + RT = Operators.operator_return_eltype(op, eltype(arg)) + Jv¹ = CUDA.CuStaticSharedArray(RT, (Nq, Nvt)) + return (Jv¹,) +end +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_fill_shmem!( + op::Operators.Divergence{(1,)}, + (Jv¹,), + space, + ij, + slabidx, + arg, +) + vt = CUDA.threadIdx().z + local_geometry = Operators.get_local_geometry(space, ij, slabidx) + i, _ = ij.I + Jv¹[i, vt] = + local_geometry.J ⊠ RecursiveApply.rmap( + v -> Geometry.contravariant1(v, local_geometry), + arg, + ) +end + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_shmem( + space, + ::Val{Nvt}, + op::Operators.WeakDivergence{(1,)}, + arg, +) where {Nvt} + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + # allocate temp output + RT = Operators.operator_return_eltype(op, eltype(arg)) + Nf = DataLayouts.typesize(FT, RT) + WJv¹ = CUDA.CuStaticSharedArray(RT, (Nq, Nvt)) + return (WJv¹,) +end +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_fill_shmem!( + op::Operators.WeakDivergence{(1,)}, + (WJv¹,), + space, + ij, + slabidx, + arg, +) + vt = CUDA.threadIdx().z + local_geometry = Operators.get_local_geometry(space, ij, slabidx) + i, _ = ij.I + WJv¹[i, vt] = + local_geometry.WJ ⊠ RecursiveApply.rmap( + v -> Geometry.contravariant1(v, local_geometry), + arg, + ) +end + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_shmem( + space, + ::Val{Nvt}, + op::Operators.Gradient{(1,)}, + arg, +) where {Nvt} + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + # allocate temp output + IT = eltype(arg) + input = CUDA.CuStaticSharedArray(IT, (Nq, Nvt)) + return (input,) +end +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_fill_shmem!( + op::Operators.Gradient{(1,)}, + (input,), + space, + ij, + slabidx, + arg, +) + vt = CUDA.threadIdx().z + i, _ = ij.I + input[i, vt] = arg +end + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_shmem( + space, + ::Val{Nvt}, + op::Operators.WeakGradient{(1,)}, + arg, +) where {Nvt} + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + # allocate temp output + IT = eltype(arg) + Wf = CUDA.CuStaticSharedArray(IT, (Nq, Nvt)) + return (Wf,) +end +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_fill_shmem!( + op::Operators.WeakGradient{(1,)}, + (Wf,), + space, + ij, + slabidx, + arg, +) + vt = CUDA.threadIdx().z + local_geometry = Operators.get_local_geometry(space, ij, slabidx) + W = local_geometry.WJ * local_geometry.invJ + i, _ = ij.I + Wf[i, vt] = W ⊠ arg +end + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_shmem( + space, + ::Val{Nvt}, + op::Operators.Curl{(1,)}, + arg, +) where {Nvt} + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + IT = eltype(arg) + ET = eltype(IT) + RT = Operators.operator_return_eltype(op, IT) + # allocate temp output + if RT <: Geometry.Contravariant3Vector + v₂ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt)) + return (nothing, v₂) + elseif RT <: Geometry.Contravariant2Vector + v₃ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt)) + return (v₃,) + elseif RT <: Geometry.Contravariant23Vector + v₂ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt)) + v₃ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt)) + return (nothing, v₂, v₃) + else + error("invalid return type") + end +end +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_fill_shmem!( + op::Operators.Curl{(1,)}, + work, + space, + ij, + slabidx, + arg, +) + vt = CUDA.threadIdx().z + i, _ = ij.I + local_geometry = Operators.get_local_geometry(space, ij, slabidx) + RT = Operators.operator_return_eltype(op, typeof(arg)) + if RT <: Geometry.Contravariant3Vector + _, v₂ = work + v₂[i, vt] = Geometry.covariant2(arg, local_geometry) + elseif RT <: Geometry.Contravariant2Vector + (v₃,) = work + v₃[i, vt] = Geometry.covariant3(arg, local_geometry) + else + _, v₂, v₃ = work + v₂[i, vt] = Geometry.covariant2(arg, local_geometry) + v₃[i, vt] = Geometry.covariant3(arg, local_geometry) + end +end + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_shmem( + space, + ::Val{Nvt}, + op::Operators.WeakCurl{(1,)}, + arg, +) where {Nvt} + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + IT = eltype(arg) + ET = eltype(IT) + RT = Operators.operator_return_eltype(op, IT) + # allocate temp output + if RT <: Geometry.Contravariant3Vector + Wv₂ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt)) + return (nothing, Wv₂) + elseif RT <: Geometry.Contravariant2Vector + Wv₃ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt)) + return (Wv₃,) + elseif RT <: Geometry.Contravariant23Vector + Wv₂ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt)) + Wv₃ = CUDA.CuStaticSharedArray(ET, (Nq, Nvt)) + return (nothing, Wv₂, Wv₃) + else + error("invalid return type") + end +end +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_fill_shmem!( + op::Operators.WeakCurl{(1,)}, + work, + space, + ij, + slabidx, + arg, +) + vt = CUDA.threadIdx().z + i, _ = ij.I + local_geometry = Operators.get_local_geometry(space, ij, slabidx) + W = local_geometry.WJ * local_geometry.invJ + RT = Operators.operator_return_eltype(op, typeof(arg)) + if RT <: Geometry.Contravariant3Vector + _, Wv₂ = work + Wv₂[i, vt] = W ⊠ Geometry.covariant2(arg, local_geometry) + elseif RT <: Geometry.Contravariant2Vector + (Wv₃,) = work + Wv₃[i, vt] = W ⊠ Geometry.covariant3(arg, local_geometry) + else + _, Wv₂, Wv₃ = work + Wv₂[i, vt] = W ⊠ Geometry.covariant2(arg, local_geometry) + Wv₃[i, vt] = W ⊠ Geometry.covariant3(arg, local_geometry) + end +end + +########## +## ClimaCore.jl/ext/cuda/operators_spectral_element.jl +########## + +function Base.copyto!( + out::Fields.Field, + sbc::Union{ + Operators.SpectralBroadcasted{ClimaCoreCUDAExt.CUDASpectralStyle}, + Base.Broadcast.Broadcasted{ClimaCoreCUDAExt.CUDASpectralStyle}, + }, +) + space = axes(out) + (Ni, Nj, _, Nv, Nh) = size(Fields.field_values(out)) + max_threads = 256 + @assert Ni * Nj ≤ max_threads + Nvthreads = fld(max_threads, Ni * Nj) + Nvblocks = cld(Nv, Nvthreads) + # executed + args = ( + Operators.strip_space(out, space), + Operators.strip_space(sbc, space), + space, + Val(Nvthreads), + ) + ClimaCoreCUDAExt.auto_launch!( + ClimaCoreCUDAExt.copyto_spectral_kernel!, + args, + out; + threads_s = (Ni, Nj, Nvthreads), + blocks_s = (Nh, Nvblocks), + ) + return out +end + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_evaluate( + op::Operators.Divergence{(1,)}, + (Jv¹,), + space, + ij, + slabidx, +) + vt = CUDA.threadIdx().z + i, _ = ij.I + + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + D = Quadratures.differentiation_matrix(FT, QS) + + local_geometry = Operators.get_local_geometry(space, ij, slabidx) + + DJv = D[i, 1] ⊠ Jv¹[1, vt] + for k in 2:Nq + DJv = DJv ⊞ D[i, k] ⊠ Jv¹[k, vt] + end + return RecursiveApply.rmul(DJv, local_geometry.invJ) +end + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_evaluate( + op::Operators.WeakDivergence{(1,)}, + (WJv¹,), + space, + ij, + slabidx, +) + vt = CUDA.threadIdx().z + i, _ = ij.I + + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + D = Quadratures.differentiation_matrix(FT, QS) + + local_geometry = Operators.get_local_geometry(space, ij, slabidx) + + Dᵀ₁WJv¹ = D[1, i] ⊠ WJv¹[1, vt] + for k in 2:Nq + Dᵀ₁WJv¹ = Dᵀ₁WJv¹ ⊞ D[k, i] ⊠ WJv¹[k, vt] + end + return ⊟(RecursiveApply.rdiv(Dᵀ₁WJv¹, local_geometry.WJ)) +end + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_evaluate( + op::Operators.Gradient{(1,)}, + (input,), + space, + ij, + slabidx, +) + vt = CUDA.threadIdx().z + i, _ = ij.I + + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + D = Quadratures.differentiation_matrix(FT, QS) + + ∂f∂ξ₁ = D[i, 1] * input[1, vt] + for k in 2:Nq + ∂f∂ξ₁ += D[i, k] * input[k, vt] + end + return Geometry.Covariant1Vector(∂f∂ξ₁) +end + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_evaluate( + op::Operators.WeakGradient{(1,)}, + (Wf,), + space, + ij, + slabidx, +) + vt = CUDA.threadIdx().z + i, _ = ij.I + + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + D = Quadratures.differentiation_matrix(FT, QS) + + local_geometry = Operators.get_local_geometry(space, ij, slabidx) + W = local_geometry.WJ * local_geometry.invJ + + Dᵀ₁Wf = D[1, i] ⊠ Wf[1, vt] + for k in 2:Nq + Dᵀ₁Wf = Dᵀ₁Wf ⊞ D[k, i] ⊠ Wf[k, vt] + end + return Geometry.Covariant1Vector(⊟(RecursiveApply.rdiv(Dᵀ₁Wf, W))) +end + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_evaluate( + op::Operators.Curl{(1,)}, + work, + space, + ij, + slabidx, +) + vt = CUDA.threadIdx().z + i, _ = ij.I + + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + D = Quadratures.differentiation_matrix(FT, QS) + local_geometry = Operators.get_local_geometry(space, ij, slabidx) + + if length(work) == 2 + _, v₂ = work + D₁v₂ = D[i, 1] ⊠ v₂[1, vt] + for k in 2:Nq + D₁v₂ = D₁v₂ ⊞ D[i, k] ⊠ v₂[k, vt] + end + return Geometry.Contravariant3Vector( + RecursiveApply.rmul(D₁v₂, local_geometry.invJ), + ) + elseif length(work) == 1 + (v₃,) = work + D₁v₃ = D[i, 1] ⊠ v₃[1, vt] + for k in 2:Nq + D₁v₃ = D₁v₃ ⊞ D[i, k] ⊠ v₃[k, vt] + end + return Geometry.Contravariant2Vector( + ⊟(RecursiveApply.rmul(D₁v₃, local_geometry.invJ)), + ) + else + _, v₂, v₃ = work + D₁v₂ = D[i, 1] ⊠ v₂[1, vt] + D₁v₃ = D[i, 1] ⊠ v₃[1, vt] + @simd for k in 2:Nq + D₁v₂ = D₁v₂ ⊞ D[i, k] ⊠ v₂[k, vt] + D₁v₃ = D₁v₃ ⊞ D[i, k] ⊠ v₃[k, vt] + end + return Geometry.Contravariant23Vector( + ⊟(RecursiveApply.rmul(D₁v₃, local_geometry.invJ)), + RecursiveApply.rmul(D₁v₂, local_geometry.invJ), + ) + end +end + +Base.@propagate_inbounds function ClimaCoreCUDAExt.operator_evaluate( + op::Operators.WeakCurl{(1,)}, + work, + space, + ij, + slabidx, +) + vt = CUDA.threadIdx().z + i, _ = ij.I + + FT = Spaces.undertype(space) + QS = Spaces.quadrature_style(space) + Nq = Quadratures.degrees_of_freedom(QS) + D = Quadratures.differentiation_matrix(FT, QS) + local_geometry = Operators.get_local_geometry(space, ij, slabidx) + + if length(work) == 2 + _, Wv₂ = work + Dᵀ₁Wv₂ = D[1, i] ⊠ Wv₂[1, vt] + for k in 2:Nq + Dᵀ₁Wv₂ = Dᵀ₁Wv₂ ⊞ D[k, i] ⊠ Wv₂[k, vt] + end + return Geometry.Contravariant3Vector( + RecursiveApply.rdiv(⊟(Dᵀ₁Wv₂), local_geometry.WJ), + ) + elseif length(work) == 1 + (Wv₃,) = work + Dᵀ₁Wv₃ = D[1, i] ⊠ Wv₃[1, vt] + for k in 2:Nq + Dᵀ₁Wv₃ = Dᵀ₁Wv₃ ⊞ D[k, i] ⊠ Wv₃[k, vt] + end + return Geometry.Contravariant2Vector( + RecursiveApply.rdiv(Dᵀ₁Wv₃, local_geometry.WJ), + ) + else + _, Wv₂, Wv₃ = work + Dᵀ₁Wv₂ = D[1, i] ⊠ Wv₂[1, vt] + Dᵀ₁Wv₃ = D[1, i] ⊠ Wv₃[1, vt] + @simd for k in 2:Nq + Dᵀ₁Wv₂ = Dᵀ₁Wv₂ ⊞ D[k, i] ⊠ Wv₂[k, vt] + Dᵀ₁Wv₃ = Dᵀ₁Wv₃ ⊞ D[k, i] ⊠ Wv₃[k, vt] + end + return Geometry.Contravariant23Vector( + RecursiveApply.rdiv(Dᵀ₁Wv₃, local_geometry.WJ), + RecursiveApply.rdiv(⊟(Dᵀ₁Wv₂), local_geometry.WJ), + ) + end +end + +########## +## ClimaCore.jl/ext/cuda/operators_finite_difference.jl +########## + +function Base.copyto!( + out::Fields.Field, + bc::Union{ + Operators.StencilBroadcasted{ClimaCoreCUDAExt.CUDAColumnStencilStyle}, + Base.Broadcast.Broadcasted{ClimaCoreCUDAExt.CUDAColumnStencilStyle}, + }, +) + space = axes(out) + (Ni, Nj, _, Nv, Nh) = size(Fields.field_values(out)) + (li, lw, rw, ri) = bounds = Operators.window_bounds(space, bc) + @assert Nv == ri - li + 1 + max_threads = 256 + nitems = Nv * Ni * Nj * Nh # # of independent items + (nthreads, nblocks) = + ClimaCoreCUDAExt._configure_threadblock(max_threads, nitems) + args = ( + Operators.strip_space(out, space), + Operators.strip_space(bc, space), + axes(out), + bounds, + Ni, + Nj, + Nh, + Nv, + ) + ClimaCoreCUDAExt.auto_launch!( + ClimaCoreCUDAExt.copyto_stencil_kernel!, + args, + out; + threads_s = (nthreads,), + blocks_s = (nblocks,), + ) + return out +end + +function ClimaCoreCUDAExt.copyto_stencil_kernel!( + out, + bc, + space, + bds, + Ni, + Nj, + Nh, + Nv, +) + gid = CUDA.threadIdx().x + (CUDA.blockIdx().x - 1) * CUDA.blockDim().x + if gid ≤ Nv * Ni * Nj * Nh + (li, lw, rw, ri) = bds + (v, i, j, h) = Utilities.cart_ind((Nv, Ni, Nj, Nh), gid).I + hidx = (i, j, h) + idx = v - 1 + li + window = + idx < lw ? + Operators.LeftBoundaryWindow{Spaces.left_boundary_name(space)}() : + ( + idx > rw ? + Operators.RightBoundaryWindow{ + Spaces.right_boundary_name(space), + }() : Operators.Interior() + ) + Operators.setidx!( + space, + out, + idx, + hidx, + Operators.getidx(space, bc, window, idx, hidx), + ) + end + return nothing +end + +########## +## ClimaCore.jl/ext/cuda/topologies_dss.jl +########## + +function Topologies.dss_1d!( + ::ClimaComms.CUDADevice, + data::Union{DataLayouts.VIFH, DataLayouts.IFH}, + topology::Topologies.IntervalTopology, + lg = nothing, + weight = nothing, +) + (_, _, _, Nv, Nh) = size(data) + Nfaces = Topologies.isperiodic(topology) ? Nh : Nh - 1 + if Nfaces > 0 + nthreads, nblocks = ClimaCoreCUDAExt._configure_threadblock(Nv * Nfaces) + ClimaCoreCUDAExt.auto_launch!( + dss_1d_kernel!, + (data, lg, weight, Nfaces), + data; + threads_s = nthreads, + blocks_s = nblocks, + ) + end + return nothing +end + +function dss_1d_kernel!(data, lg, weight, Nfaces) + T = eltype(data) + (Ni, _, _, Nv, Nh) = size(data) + gidx = CUDA.threadIdx().x + (CUDA.blockIdx().x - 1) * CUDA.blockDim().x + if gidx ≤ Nv * Nfaces + left_face_elem = cld(gidx, Nv) + level = gidx - (left_face_elem - 1) * Nv + right_face_elem = (left_face_elem % Nh) + 1 + left_idx = CartesianIndex(Ni, 1, 1, level, left_face_elem) + right_idx = CartesianIndex(1, 1, 1, level, right_face_elem) + val = + Topologies.dss_transform(data, lg, weight, left_idx) ⊞ + Topologies.dss_transform(data, lg, weight, right_idx) + data[left_idx] = Topologies.dss_untransform(T, val, lg, left_idx) + data[right_idx] = Topologies.dss_untransform(T, val, lg, right_idx) + end + return nothing +end + +########## +## ClimaCore.jl/src/Topologies/dss.jl +########## + +function Topologies.dss_1d!( + ::ClimaComms.AbstractCPUDevice, + data::Union{DataLayouts.VIFH, DataLayouts.IFH}, + topology::Topologies.IntervalTopology, + lg = nothing, + weight = nothing, +) + T = eltype(data) + (Ni, _, _, Nv, Nh) = size(data) + Nfaces = Topologies.isperiodic(topology) ? Nh : Nh - 1 + @inbounds for left_face_elem in 1:Nfaces, level in 1:Nv + right_face_elem = (left_face_elem % Nh) + 1 + left_idx = CartesianIndex(Ni, 1, 1, level, left_face_elem) + right_idx = CartesianIndex(1, 1, 1, level, right_face_elem) + val = + Topologies.dss_transform(data, lg, weight, left_idx) ⊞ + Topologies.dss_transform(data, lg, weight, right_idx) + data[left_idx] = Topologies.dss_untransform(T, val, lg, left_idx) + data[right_idx] = Topologies.dss_untransform(T, val, lg, right_idx) + end +end + +########## +## ClimaCore.jl/src/Spaces/dss.jl +########## + +function Spaces.weighted_dss_internal!( + data::Union{ + DataLayouts.IFH, + DataLayouts.VIFH, + DataLayouts.IJFH, + DataLayouts.VIJFH, + }, + space::Union{ + Spaces.AbstractSpectralElementSpace, + Spaces.ExtrudedFiniteDifferenceSpace, + }, + hspace::Spaces.AbstractSpectralElementSpace, + dss_buffer::Union{Topologies.DSSBuffer, Nothing}, +) + Spaces.assert_same_eltype(data, dss_buffer) + length(parent(data)) == 0 && return nothing + device = ClimaComms.device(Spaces.topology(hspace)) + if hspace isa Spaces.SpectralElementSpace1D + Topologies.dss_1d!( + device, + data, + Spaces.topology(hspace), + Spaces.local_geometry_data(space), + Spaces.local_dss_weights(space), + ) + else + Topologies.dss_transform!( + device, + dss_buffer, + data, + Spaces.local_geometry_data(space), + Spaces.local_dss_weights(space), + Spaces.perimeter(hspace), + dss_buffer.internal_elems, + ) + Topologies.dss_local!( + device, + dss_buffer.perimeter_data, + Spaces.perimeter(hspace), + Spaces.topology(hspace), + ) + Topologies.dss_untransform!( + device, + dss_buffer, + data, + Spaces.local_geometry_data(space), + Spaces.perimeter(hspace), + dss_buffer.internal_elems, + ) + end + return nothing +end + +########## +## ClimaCore.jl/src/Grids/spectralelement.jl +########## + +function Grids._SpectralElementGrid1D( + topology::Topologies.IntervalTopology, + quadrature_style::Quadratures.QuadratureStyle, +) + DA = ClimaComms.array_type(topology) + global_geometry = Geometry.CartesianGlobalGeometry() + CoordType = Topologies.coordinate_type(topology) + AIdx = Geometry.coordinate_axis(CoordType) + FT = eltype(CoordType) + nelements = Topologies.nlocalelems(topology) + Nq = Quadratures.degrees_of_freedom(quadrature_style) + LG = Geometry.LocalGeometry{AIdx, CoordType, FT, SMatrix{1, 1, FT, 1}} + quad_points, quad_weights = + Quadratures.quadrature_points(FT, quadrature_style) + + # Initialize the local_geometry with Array{FT} as the backing array type. + local_geometry_cpu = DataLayouts.IFH{LG, Nq}(Array{FT}, nelements) + + for elem in 1:nelements + local_geometry_slab = slab(local_geometry_cpu, elem) + for i in 1:Nq + ξ = quad_points[i] + vcoords = Topologies.vertex_coordinates(topology, elem) + x = Geometry.linear_interpolate(vcoords, ξ) + ∂x∂ξ = + ( + Geometry.component(vcoords[2], 1) - + Geometry.component(vcoords[1], 1) + ) / 2 + J = abs(∂x∂ξ) + WJ = J * quad_weights[i] + local_geometry_slab[i] = Geometry.LocalGeometry( + x, + J, + WJ, + Geometry.AxisTensor( + ( + Geometry.LocalAxis{AIdx}(), + Geometry.CovariantAxis{AIdx}(), + ), + ∂x∂ξ, + ), + ) + end + end + + # If needed, move the local_geometry onto the GPU. + local_geometry = DataLayouts.rebuild(local_geometry_cpu, DA) + + # TODO: Why was this code setting dss_weights = 1 / DSS(1)? That's just 1?? + # dss_weights = copy(local_geometry.J) + # dss_weights .= one(FT) + # Topologies.dss_1d!(topology, dss_weights) + # dss_weights = one(FT) ./ dss_weights + + dss_weights = copy(local_geometry.J) + Topologies.dss_1d!(ClimaComms.device(topology), dss_weights, topology) + dss_weights .= local_geometry.J ./ dss_weights + + return Grids.SpectralElementGrid1D( + topology, + quadrature_style, + global_geometry, + local_geometry, + dss_weights, + ) +end + +########## +## ClimaCore.jl/src/Hypsography/Hypsography.jl +########## + +function Grids._ExtrudedFiniteDifferenceGrid( + horizontal_grid::Grids.AbstractGrid, + vertical_grid::Grids.FiniteDifferenceGrid, + adaption::Grids.HypsographyAdaption, + global_geometry::Geometry.AbstractGlobalGeometry, +) + @assert Spaces.grid(axes(adaption.surface)) == horizontal_grid + z_surface = Fields.field_values(adaption.surface) + + face_z_ref = + Grids.local_geometry_data(vertical_grid, Grids.CellFace()).coordinates + vertical_domain = Topologies.domain(vertical_grid) + z_top = vertical_domain.coord_max + + face_z = + modified_ref_z_to_physical_z.( + Ref(typeof(adaption)), + face_z_ref, + z_surface, + Ref(z_top), + Ref(adaption_parameters(adaption)), + ) + face_z_cpu = + modified_ref_z_to_physical_z.( + Ref(typeof(adaption)), + Adapt.adapt(Array, face_z_ref), + Adapt.adapt(Array, z_surface), + Ref(z_top), + Ref(adaption_parameters(adaption)), + ) + + return Grids._ExtrudedFiniteDifferenceGrid( + horizontal_grid, + vertical_grid, + adaption, + global_geometry, + face_z, + ) +end + +adaption_parameters(::Grids.Flat) = (;) +adaption_parameters(::Hypsography.LinearAdaption) = (;) +adaption_parameters((; ηₕ, s)::Hypsography.SLEVEAdaption) = (; ηₕ, s) + +function modified_ref_z_to_physical_z( + ::Type{<:Grids.Flat}, + z_ref::Geometry.ZPoint, + z_surface::Geometry.ZPoint, + z_top::Geometry.ZPoint, + _, +) + return z_ref +end +function modified_ref_z_to_physical_z( + ::Type{<:Hypsography.LinearAdaption}, + z_ref::Geometry.ZPoint, + z_surface::Geometry.ZPoint, + z_top::Geometry.ZPoint, + _, +) + Geometry.ZPoint(z_ref.z + (1 - z_ref.z / z_top.z) * z_surface.z) +end +function modified_ref_z_to_physical_z( + ::Type{<:Hypsography.SLEVEAdaption}, + z_ref::Geometry.ZPoint, + z_surface::Geometry.ZPoint, + z_top::Geometry.ZPoint, + (; ηₕ, s), +) + if s * z_top.z <= z_surface.z + error("Decay scale (s*z_top) must be higher than max surface elevation") + end + + η = z_ref.z / z_top.z + if η <= ηₕ + return Geometry.ZPoint( + η * z_top.z + + z_surface.z * (sinh((ηₕ - η) / s / ηₕ)) / (sinh(1 / s)), + ) + else + return Geometry.ZPoint(η * z_top.z) + end +end diff --git a/examples/hybrid/fix_asymmetric_metric_tensor.jl b/examples/hybrid/fix_asymmetric_metric_tensor.jl new file mode 100644 index 00000000000..7e13c95ea69 --- /dev/null +++ b/examples/hybrid/fix_asymmetric_metric_tensor.jl @@ -0,0 +1,25 @@ +import ClimaCore: Geometry, DataLayouts +import LinearAlgebra: issymmetric + +@inline function Geometry.LocalGeometry( + coordinates, + J, + WJ, + ∂x∂ξ::Geometry.Axis2Tensor{ + FT, + Tuple{Geometry.LocalAxis{I}, Geometry.CovariantAxis{I}}, + S, + }, +) where {FT, I, S} + ∂ξ∂x = inv(∂x∂ξ) + C = typeof(coordinates) + Jinv = inv(J) + gⁱʲ = (∂ξ∂x * ∂ξ∂x' + (∂ξ∂x * ∂ξ∂x')') / 2 + gᵢⱼ = (∂x∂ξ' * ∂x∂ξ + (∂x∂ξ' * ∂x∂ξ)') / 2 + issymmetric(Geometry.components(gⁱʲ)) || error("gⁱʲ is not symmetric.") + issymmetric(Geometry.components(gᵢⱼ)) || error("gᵢⱼ is not symmetric.") + return DataLayouts.bypass_constructor( + Geometry.LocalGeometry{I, C, FT, S}, + (coordinates, J, WJ, Jinv, ∂x∂ξ, ∂ξ∂x, gⁱʲ, gᵢⱼ), + ) +end diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl index 8b52d5f74cf..d9418ef7b75 100644 --- a/post_processing/ci_plots.jl +++ b/post_processing/ci_plots.jl @@ -42,13 +42,15 @@ import CairoMakie import CairoMakie.Makie import ClimaAnalysis import ClimaAnalysis: Visualize as viz -import ClimaAnalysis: SimDir, slice, read_var +import ClimaAnalysis: SimDir, slice, window, read_var import ClimaAnalysis.Utils: kwargs as ca_kwargs import ClimaCoreSpectra: power_spectrum_2d using Poppler_jll: pdfunite import Base.Filesystem +import LinearAlgebra: norm +import StaticArrays: SVector const days = 86400 @@ -279,6 +281,90 @@ function make_plots_generic( return output_file end +# TODO: Move the following to ClimaAnalysis. +import Interpolations +linear_interpolation(ranges, values) = Interpolations.linear_interpolation( + ranges, + values; + extrapolation_bc = Interpolations.Line(), +) +function plot_streamlines!(place, var) + length(var.dims) == 2 || error("Can only plot 2D variables") + + dim1_name, dim2_name = var.index2dim + dim1 = var.dims[dim1_name] + dim2 = var.dims[dim2_name] + dim1_units = var.dim_attributes[dim1_name]["units"] + dim2_units = var.dim_attributes[dim2_name]["units"] + + var_units = var.attributes["units"] + var_name = var.attributes["short_name"] + + CairoMakie.Axis( + place[1, 1]; + title = var.attributes["long_name"], + xlabel = "$dim1_name [$dim1_units]", + ylabel = "$dim2_name [$dim2_units]", + xgridvisible = false, + ygridvisible = false, + limits = (extrema(dim1), extrema(dim2)), + ) + + u_spline = linear_interpolation((dim1, dim2), first.(var.data)) + w_spline = linear_interpolation((dim1, dim2), last.(var.data)) + spline = point -> Makie.Point2f(u_spline(point...), w_spline(point...)) + plot = CairoMakie.streamplot!( + spline, + dim1, + dim2; + stepsize = 0.4, + maxsteps = 60000, + gridsize = (100, 20), + arrow_size = 7, + rasterize = 10, + ) + + CairoMakie.Colorbar(place[1, 2], plot; label = "$var_name [$var_units]") +end +function plot_contours!(place, var) + length(var.dims) == 2 || error("Can only plot 2D variables") + + dim1_name, dim2_name = var.index2dim + dim1 = var.dims[dim1_name] + dim2 = var.dims[dim2_name] + dim1_units = var.dim_attributes[dim1_name]["units"] + dim2_units = var.dim_attributes[dim2_name]["units"] + + var_name = var.attributes["short_name"] + var_units = var.attributes["units"] + + CairoMakie.Axis( + place[1, 1]; + title = var.attributes["long_name"], + xlabel = "$dim1_name [$dim1_units]", + ylabel = "$dim2_name [$dim2_units]", + limits = (extrema(dim1), extrema(dim2)), + ) + + n_levels = 21 + var_mean = round(mean(var.data)) + var_delta = maximum(value -> abs(value - var_mean), var.data) + levels = range(var_mean - var_delta, var_mean + var_delta, n_levels + 1) + + spectral_colors = Makie.to_colormap(:Spectral) + plot_kwargs = (; + levels, + colormap = setindex!(spectral_colors, Makie.RGBA(1, 1, 1, 0), 6), + extendhigh = spectral_colors[11], + extendlow = spectral_colors[1], + ) + plot = CairoMakie.contourf!(dim1, dim2, var.data; plot_kwargs...) + + if var_delta / abs(var_mean) > 1e-7 + CairoMakie.Colorbar(place[1, 2], plot; label = "$var_name [$var_units]") + end # Colorbar's LineAxis throws an error when var_delta ≪ abs(var_mean). +end + """ compute_spectrum(var::ClimaAnalysis.OutputVar; mass_weight = nothing) @@ -589,6 +675,219 @@ function make_plots( make_plots_generic(output_paths, vars, y = 0.0, time = LAST_SNAP) end +rms(var; dims) = sqrt.(mean(x -> x^2, var; dims)) +function level_rms(var) + reduced_var = ClimaAnalysis.Var._reduce_over(rms, "x", var) + if haskey(var.dims, "y") + reduced_var = ClimaAnalysis.Var._reduce_over(rms, "y", reduced_var) + end + if haskey(var.attributes, "long_name") + long_name = reduced_var.attributes["long_name"] + reduced_var.attributes["long_name"] = + "RMS " * long_name * " averaged over levels" + end + return reduced_var +end +function column_rms(var) + reduced_var = ClimaAnalysis.Var._reduce_over(rms, "z_reference", var) + if haskey(var.attributes, "long_name") + long_name = reduced_var.attributes["long_name"] + reduced_var.attributes["long_name"] = + "RMS " * long_name * " averaged over columns" + end + return reduced_var +end + +const AnalyticMountainTest = Union{ + Val{:gpu_plane_analytic_agnesi_mountain_test}, + Val{:gpu_plane_analytic_huge_schar_mountain_test}, + Val{:gpu_plane_analytic_big_schar_mountain_test}, + Val{:gpu_plane_analytic_schar_mountain_test}, + Val{:gpu_plane_analytic_small_schar_mountain_test}, + Val{:gpu_plane_analytic_tiny_schar_mountain_test}, + Val{:gpu_plane_analytic_teeny_tiny_schar_mountain_test}, + Val{:gpu_plane_analytic_schar_mountain_float32_test}, + Val{:gpu_plane_analytic_tiny_schar_mountain_float32_test}, + Val{:gpu_plane_analytic_schar_mountain_no_hyperdiff_test}, + Val{:gpu_plane_analytic_schar_mountain_no_hyperdiff_float32_test}, + Val{:gpu_plane_analytic_schar_mountain_stretched_grid_test}, + Val{:gpu_plane_analytic_schar_mountain_high_horz_res_test}, + Val{:gpu_plane_analytic_schar_mountain_high_vert_res_test}, + Val{:gpu_plane_analytic_schar_mountain_high_top_test}, + Val{:gpu_plane_analytic_schar_mountain_high_top_high_sponge_test}, + Val{:gpu_plane_analytic_schar_mountain_high_top_weak_high_sponge_test}, +} +const AnalyticTopographyTest = Union{ + Val{:plane_analytic_no_topography_test}, + Val{:plane_analytic_no_topography_float32_test}, + Val{:gpu_plane_analytic_no_topography_long_test}, + Val{:gpu_plane_analytic_no_topography_long_float32_test}, + Val{:gpu_plane_analytic_no_topography_no_hyperdiff_long_test}, + Val{:gpu_plane_analytic_no_topography_no_hyperdiff_long_float32_test}, + Val{:gpu_plane_analytic_no_topography_discrete_balance_long_test}, + Val{:gpu_plane_analytic_no_topography_discrete_balance_long_float32_test}, + Val{:gpu_extruded_plane_analytic_no_topography_test}, + Val{:gpu_box_analytic_no_topography_test}, + Val{:gpu_box_analytic_no_topography_float32_test}, + Val{:plane_analytic_cosine_hills_test}, + Val{:plane_analytic_cosine_hills_float32_test}, + Val{:gpu_plane_analytic_cosine_hills_long_test}, + Val{:gpu_plane_analytic_cosine_hills_long_float32_test}, + Val{:gpu_plane_analytic_cosine_hills_strong_sponge_long_test}, + Val{:gpu_plane_analytic_cosine_hills_weak_sponge_long_test}, + Val{:gpu_plane_analytic_cosine_hills_high_top_long_test}, + Val{:gpu_plane_analytic_cosine_hills_high_top_high_sponge_long_test}, + Val{:gpu_extruded_plane_analytic_cosine_hills_test}, + Val{:gpu_box_analytic_cosine_hills_test}, + Val{:gpu_box_analytic_cosine_hills_float32_test}, + AnalyticMountainTest, +} + +function make_plots( + val::AnalyticTopographyTest, + output_paths::Vector{<:AbstractString}, +) + simdirs = SimDir.(output_paths) + + is_mountain_test = val isa AnalyticMountainTest + is_3d = val in ( + Val(:gpu_extruded_plane_analytic_no_topography_test), + Val(:gpu_box_analytic_no_topography_test), + Val(:gpu_box_analytic_no_topography_float32_test), + Val(:gpu_extruded_plane_analytic_cosine_hills_test), + Val(:gpu_box_analytic_cosine_hills_test), + Val(:gpu_box_analytic_cosine_hills_float32_test), + ) + + z_bounds = (; left = 0, right = 13e3) # The Rayleigh sponge starts at 13 km. + + error_short_name_pairs = [["uaerror", "waerror"], ["c1error", "c3error"]] + rms_error_vars = Iterators.flatmap(error_short_name_pairs) do short_names + Iterators.flatmap((level_rms, column_rms)) do rms_func + Iterators.flatmap(short_names) do short_name + Iterators.map(simdirs) do simdir + data = get(simdir; short_name) + data = window(data, "z_reference"; z_bounds...) + data = slice(data; time = Inf) + rms_func(data) + end + end + end + end + orog_vars = map(simdirs) do simdir + slice(get(simdir; short_name = "orog"); time = Inf) + end + make_plots_generic( + output_paths, + [rms_error_vars..., orog_vars...], + output_name = "final_rms_errors", + ) + + all_short_names = [ + "ua", + "wa", + "uaerror", + "waerror", + "c1error", + "c3error", + "uapredicted", + "wapredicted", + ] + slice_summary_vars = Iterators.flatmap(all_short_names) do short_name + hours_to_plot = + endswith(short_name, "predicted") ? (Inf,) : (1, 2, 24, Inf) + Iterators.flatmap(hours_to_plot) do n_hours + Iterators.map(simdirs) do simdir + data = get(simdir; short_name) + data = window(data, "z_reference"; z_bounds...) + data = slice(data; time = n_hours * 60 * 60) + is_3d ? slice(data; y = 0) : data + end + end + end + make_plots_generic( + output_paths, + collect(slice_summary_vars); + plot_fn = plot_contours!, + output_name = "slice_summary", + ) + + if is_mountain_test + vertical_data_short_names = ["wa", "wapredicted", "waerror", "c3error"] + mountain_closeup_vars = + Iterators.flatmap(vertical_data_short_names) do short_name + Iterators.map(simdirs) do simdir + data = get(simdir; short_name) + data = window(data, "z_reference"; z_bounds...) + data = window(data, "x"; left = 35e3, right = 65e3) + data = slice(data; time = Inf) + is_3d ? slice(data; y = 0) : data + end + end + make_plots_generic( + output_paths, + collect(mountain_closeup_vars); + plot_fn = plot_contours!, + output_name = "final_mountain_closeup", + ) + end + + # collected_vars = collect(slice_summary_vars) + # (; dims, dim_attributes, dim2index, index2dim) = collected_vars[1] + # u_data = collected_vars[4].data + # w_data = collected_vars[8].data + # predicted_u_data = collected_vars[9].data + # predicted_w_data = collected_vars[10].data + + # FT = eltype(predicted_u_data) + # background_velocity = SVector(CA.background_u(FT), FT(0)) + # perturbation_data = @. SVector(u_data, w_data) - (background_velocity,) + # predicted_perturbation_data = + # @. SVector(predicted_u_data, predicted_w_data) - (background_velocity,) + + # # Rescale the predicted perturbations so that their range no less than the + # # background air speed. Apply the same scaling to the actual perturbations. + # background_speed = norm(background_velocity) + # predicted_perturbation_range = + # -(reverse(extrema(norm, predicted_perturbation_data))...) + # s = max(1, round(Int, background_speed / predicted_perturbation_range)) + # uw_data = @. s * perturbation_data + (background_velocity,) + # predicted_uw_data = + # @. s * predicted_perturbation_data + (background_velocity,) + + # uw_var = ClimaAnalysis.OutputVar( + # Dict{String, String}( + # "units" => "m s^-1", + # "short_name" => "ū + $s u′", + # "long_name" => "Air Speed with $(s)x Perturbations, Instantaneous", + # ), + # dims, + # dim_attributes, + # uw_data, + # dim2index, + # index2dim, + # ) + # predicted_uw_var = ClimaAnalysis.OutputVar( + # Dict{String, String}( + # "units" => "m s^-1", + # "short_name" => "predicted ū + $s u′", + # "long_name" => "Predicted Air Speed with $(s)x Perturbations, Instantaneous", + # ), + # dims, + # dim_attributes, + # predicted_uw_data, + # dim2index, + # index2dim, + # ) + + # make_plots_generic( + # output_paths, + # [uw_var, predicted_uw_var]; + # plot_fn = plot_streamlines!, + # output_name = "streamlines_24h", + # ) +end + MountainPlots = Union{ Val{:plane_agnesi_mountain_test_uniform}, Val{:plane_agnesi_mountain_test_stretched}, diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl index 85d5545cec1..d67f87a1e0e 100644 --- a/src/ClimaAtmos.jl +++ b/src/ClimaAtmos.jl @@ -16,6 +16,7 @@ include(joinpath("utils", "debug_utils.jl")) include(joinpath("topography", "topography.jl")) include(joinpath("utils", "variable_manipulations.jl")) include(joinpath("utils", "read_gcm_driven_scm_data.jl")) +include(joinpath("analytic_solutions", "analytic_solutions.jl")) include( joinpath("parameterized_tendencies", "radiation", "radiation_utilities.jl"), diff --git a/src/analytic_solutions/analytic_solutions.jl b/src/analytic_solutions/analytic_solutions.jl new file mode 100644 index 00000000000..c6656f1c169 --- /dev/null +++ b/src/analytic_solutions/analytic_solutions.jl @@ -0,0 +1,309 @@ +#= +References: +https://www.cosmo-model.org/content/model/documentation/newsLetters/newsLetter09/cnl9-04.pdf +http://dx.doi.org/10.1175/1520-0493(2003)131%3C1229:NCOMTI%3E2.0.CO;2 +https://atmos.uw.edu/academics/classes/2010Q1/536/1503AP_lee_waves.pdf +=# +import StaticArrays: @SVector + +background_N(::Type{FT}) where {FT} = FT(0.01) +background_T_sfc(::Type{FT}) where {FT} = FT(288) +background_u(::Type{FT}) where {FT} = FT(10) + +function initial_thermo_state(params, coord) + FT = eltype(params) + g = CAP.grav(params) + R_d = CAP.R_d(params) + cp_d = CAP.cp_d(params) + p_sfc = CAP.MSLP(params) + N = background_N(FT) + T_sfc = background_T_sfc(FT) + + β = N^2 / g + a = g / (cp_d * T_sfc * β) + + # Replace ζ with z in the constant-N profile so that the initial condition + # is hydrostatically balanced. + z = coord.z + + # Replace the constant-N profile with an isothermal profile above T = 100 to + # avoid unreasonably small or negative values of T. + T_min = FT(100) + z_iso = log((T_min / T_sfc - a) / (1 - a)) / β + p_iso = p_sfc * ((1 - a) / (1 - a * T_sfc / T_min))^(cp_d / R_d) + + # The perturbation in approximate_FΔuvw is computed around a background + # state that is not hydrostatically balanced and is unphysical above 30 km, + # but the initial condition does not need to exactly match this state. + p = + z > z_iso ? p_iso * exp(-g / (R_d * T_min) * (z - z_iso)) : + p_sfc * (1 - a + a * exp(-β * z))^(cp_d / R_d) + T = z > z_iso ? T_min : T_sfc * ((1 - a) * exp(β * z) + a) + + return TD.PhaseDry_pT(CAP.thermodynamics_params(params), p, T) +end + +function approximate_FΔuvw(params, k_x, k_y, ζ, Fh) + FT = eltype(params) + g = CAP.grav(params) + R_d = CAP.R_d(params) + cp_d = CAP.cp_d(params) + p_sfc = CAP.MSLP(params) + N = background_N(FT) + T_sfc = background_T_sfc(FT) + u = background_u(FT) + (; z_top) = CAP.topography_params(params) + + β = N^2 / g + a = g / (cp_d * T_sfc * β) + γ = cp_d / (cp_d - R_d) + α = g / (γ * R_d * T_sfc * β) + r = k_y / k_x + k_h² = k_x^2 + k_y^2 + + p = p_sfc * (1 - a + a * exp(-β * ζ))^(cp_d / R_d) + T = T_sfc * ((1 - a) * exp(β * ζ) + a) + + FΔζ_xh = -im * k_x * Fh # ≈ F(-∂h/∂x / (1 - h/z_top)) to first order in h + FΔζ_yh = -im * k_y * Fh # ≈ F(-∂h/∂y / (1 - h/z_top)) to first order in h + FΔζ_z = Fh / z_top # ≈ F(h/z_top / (1 - h/z_top)) to first order in h + + ρ_sfc = p_sfc / (R_d * T_sfc) + m_sfc = u^2 / (γ * R_d * T_sfc) + μ_sfc = 1 - k_x^2 / k_h² * m_sfc + ν_sfc = k_x^2 / k_h² * m_sfc / μ_sfc + d_sfc = ρ_sfc / μ_sfc + FΔζ_x_sfc = FΔζ_xh + FΔζ_y_sfc = FΔζ_yh + + ρ = p / (R_d * T) + m = u^2 / (γ * R_d * T) + μ = 1 - k_x^2 / k_h² * m + ν = k_x^2 / k_h² * m / μ + d = ρ / μ + FΔζ_x = FΔζ_xh * (1 - ζ / z_top) + FΔζ_y = FΔζ_yh * (1 - ζ / z_top) + + k_ζ_sfc² = + k_h² * (-μ_sfc + (N / (k_x * u))^2 * (1 + ν_sfc * (1 - a))) - + β^2 / 4 * ( + (1 + α)^2 + + 2 / μ_sfc * α * (1 - a) + + ν_sfc * (1 + 3 / μ_sfc) * (1 - a)^2 + ) + Δf_sfc = + g / u * ( + im * k_h² / k_x * μ_sfc * FΔζ_z + + ((1 - m_sfc) * FΔζ_xh + r * FΔζ_yh) / z_top + + β * r * (-r * FΔζ_x_sfc + FΔζ_y_sfc) * (1 + ν_sfc * (1 - a)) + ) + + k_ζ² = + k_h² * (-μ + (N / (k_x * u))^2 * (1 + ν * (1 - a))) - + β^2 / 4 * + ((1 + α)^2 + 2 / μ * α * (1 - a) + ν * (1 + 3 / μ) * (1 - a)^2) + Δf = + g / u * ( + im * k_h² / k_x * μ * FΔζ_z + + ((1 - m) * FΔζ_xh + r * FΔζ_yh) / z_top + + β * r * (-r * FΔζ_x + FΔζ_y) * (1 + ν * (1 - a)) + ) + + plus_or_minus = k_ζ² > 0 ? sign(k_x * u) : 1 + + # Approximate FΔw and ∂FΔw_∂ζ to zeroth order in β. + FΔw = + Δf / k_ζ² + + (sqrt(d_sfc / d) * (im * k_x * u * Fh - Δf_sfc / k_ζ_sfc²)) * + exp(plus_or_minus * im * sqrt(Complex(k_ζ²)) * ζ) + ∂FΔw_∂ζ = plus_or_minus * im * sqrt(Complex(k_ζ²)) * (FΔw - Δf / k_ζ²) + + FΔp = + -(im * k_x / k_h² * g * d) * + ((1 - m) * FΔζ_x + r * FΔζ_y - m / u * FΔw + u / g * ∂FΔw_∂ζ) + FΔu = -(im / k_x * g * ρ * FΔζ_x + FΔp) / (ρ * u) + FΔv = -(im / k_x * g * ρ * FΔζ_y + r * FΔp) / (ρ * u) + + return @SVector([FΔu, FΔv, FΔw]) +end + +# Integrate f(x) from x1 to x2 with n_calls function calls. Note that QuadGK is +# not GPU-compatible, so it should not be used here. +function definite_integral(f::F, x1, x2, n_calls) where {F} + dx = (x2 - x1) / n_calls + return sum(index -> f(x1 + index * dx), 2:n_calls; init = f(x1 + dx)) * dx +end + +function predicted_velocity_mountain_2d( + params, + coord, + x_center, + k_x_max, + mountain_topography::F1, + centered_mountain_topography_transform::F2, +) where {F1, F2} + FT = eltype(params) + u = background_u(FT) + (; z_top) = CAP.topography_params(params) + (; x, z) = coord + + h = mountain_topography(params, coord) + ζ = (z - h) / (1 - h / z_top) + + # Compute the inverse Fourier transform of FΔuvw. Since the integral is + # symmetric, we can just evaluate it over positive values of k_x and + # multiply the result by 2. The integral should go to k_x = Inf, but we can + # assume that it is negligibly small when k_x > k_x_max. + (Δu, Δv, Δw) = 2 * definite_integral(FT(0), k_x_max, 50000) do k_x + # Use the transform of the centered topography, h(x - x_center), + # instead of the actual topography, h(x), to avoid propagating the + # quantity exp(-im * k_x * x_center) through approximate_FΔuvw, as + # this can introduce numerical errors. + Fh = centered_mountain_topography_transform(k_x) + FΔuvw = approximate_FΔuvw(params, k_x, 0, ζ, Fh) + real.(FΔuvw * exp(im * k_x * (x - x_center))) + end + + return UVW(u + Δu, Δv, Δw) +end + +################################################################################ + +agnesi_topography_params(::Type{FT}) where {FT} = + (; a = FT(5e3), x_center = FT(50e3)) +schar_topography_params(::Type{FT}) where {FT} = + (; λ = FT(4e3), a = FT(5e3), x_center = FT(50e3)) +cosine_topography_params(::Type{FT}) where {FT} = (; λ = FT(100e3)) + +################################################################################ + +function topography_agnesi(params, coord) + FT = eltype(params) + (; h_max) = CAP.topography_params(params) + (; a, x_center) = agnesi_topography_params(FT) + (; x) = coord + + return h_max / (1 + ((x - x_center) / a)^2) +end + +function predicted_velocity_agnesi(params, coord) + FT = eltype(params) + (; h_max) = CAP.topography_params(params) + (; a, x_center) = agnesi_topography_params(FT) + + centered_topography_transform_agnesi(k_x) = + h_max * a / 2 * exp(-a * abs(k_x)) + + n_efolding_intervals = -log(eps(FT)) + k_x_max = n_efolding_intervals / a + # centered_topography_transform_agnesi(k_x) < eps(FT) when k_x > k_x_max + + return predicted_velocity_mountain_2d( + params, + coord, + x_center, + k_x_max, + topography_agnesi, + centered_topography_transform_agnesi, + ) +end + +################################################################################ + +function topography_schar(params, coord) + FT = eltype(params) + (; h_max) = CAP.topography_params(params) + (; λ, a, x_center) = schar_topography_params(FT) + (; x) = coord + + return h_max * exp(-(x - x_center)^2 / a^2) * cospi((x - x_center) / λ)^2 +end + +function predicted_velocity_schar(params, coord) + FT = eltype(params) + (; h_max) = CAP.topography_params(params) + (; λ, a, x_center) = schar_topography_params(FT) + + k_peak = 2 * FT(π) / λ + centered_topography_transform_schar(k_x::FT) where {FT} = + h_max * a / (8 * sqrt(FT(π))) * ( + exp(-a^2 / 4 * (k_x + k_peak)^2) + + 2 * exp(-a^2 / 4 * k_x^2) + + exp(-a^2 / 4 * (k_x - k_peak)^2) + ) # FT needs to be redefined for type stability + + n_efolding_intervals = -log(eps(FT)) + k_x_max = k_peak + 2 * sqrt(n_efolding_intervals) / a + # centered_topography_transform_schar(k_x) < eps(FT) when k_x > k_x_max + + return predicted_velocity_mountain_2d( + params, + coord, + x_center, + k_x_max, + topography_schar, + centered_topography_transform_schar, + ) +end + +################################################################################ + +function topography_cosine_2d(params, coord) + FT = eltype(params) + (; h_max) = CAP.topography_params(params) + (; λ) = cosine_topography_params(FT) + (; x) = coord + + return h_max * cospi(2 * x / λ) +end + +function predicted_velocity_cosine_2d(params, coord) + FT = eltype(params) + u = background_u(FT) + (; h_max, z_top) = CAP.topography_params(params) + (; λ) = cosine_topography_params(FT) + (; x, z) = coord + + h = topography_cosine_2d(params, coord) + ζ = (z - h) / (1 - h / z_top) + + # Instead of using Fh = h_max / 2 * (δ(k_x + k_peak) + δ(k_x - k_peak)) and + # integrating FΔuvw * exp(im * k_x * x) over all values of k_x, we can drop + # the delta functions and just evaluate the integrand at k_x = ±k_peak. + k_peak = 2 * FT(π) / λ + Fh_δ_coef = h_max / 2 + FΔuvw_δ_coef = approximate_FΔuvw(params, k_peak, 0, ζ, Fh_δ_coef) + (Δu, Δv, Δw) = 2 * real.(FΔuvw_δ_coef * exp(im * k_peak * x)) + + return UVW(u + Δu, Δv, Δw) +end + +################################################################################ + +function topography_cosine_3d(params, coord) + FT = eltype(params) + (; h_max) = CAP.topography_params(params) + (; λ) = cosine_topography_params(FT) + (; x, y) = coord + + return h_max * cospi(2 * x / λ) * cospi(2 * y / λ) +end + +function predicted_velocity_cosine_3d(params, coord) + FT = eltype(params) + u = background_u(FT) + (; h_max, z_top) = CAP.topography_params(params) + (; λ) = cosine_topography_params(FT) + (; x, y, z) = coord + + h = topography_cosine_3d(params, coord) + ζ = (z - h) / (1 - h / z_top) + + # Drop all delta functions and evaluate the integrand at k_x, k_y = ±k_peak. + k_peak = 2 * FT(π) / λ + Fh_δ_coef = h_max / 4 + FΔuvw_δ_coef = approximate_FΔuvw(params, k_peak, k_peak, ζ, Fh_δ_coef) + (Δu, Δv, Δw) = 4 * real.(FΔuvw_δ_coef * exp(im * k_peak * (x + y))) + + return UVW(u + Δu, Δv, Δw) +end diff --git a/src/cache/cache.jl b/src/cache/cache.jl index 6c1feb4e016..ab5f2cef99f 100644 --- a/src/cache/cache.jl +++ b/src/cache/cache.jl @@ -28,6 +28,7 @@ struct AtmosCache{ TRAC, NETFLUXTOA, NETFLUXSFC, + PSS, CONSCHECK, OD, } @@ -94,6 +95,9 @@ struct AtmosCache{ net_energy_flux_toa::NETFLUXTOA net_energy_flux_sfc::NETFLUXSFC + """Components of the predicted steady-state solution, if one is available""" + predicted_steady_state::PSS + """Conservation check for prognostic surface temperature""" conservation_check::CONSCHECK @@ -113,7 +117,15 @@ end # The model also depends on f_plane_coriolis_frequency(params) # This is a constant Coriolis frequency that is only used if space is flat -function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names) +function build_cache( + Y, + atmos, + params, + surface_setup, + sim_info, + aerosol_names, + predicted_steady_state, +) (; dt, t_end, start_date, output_dir) = sim_info FT = eltype(params) @@ -252,6 +264,7 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names) tracers, net_energy_flux_toa, net_energy_flux_sfc, + predicted_steady_state, conservation_check, output_dir, ) diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl index 520b6ac9424..e0c560c0293 100644 --- a/src/diagnostics/core_diagnostics.jl +++ b/src/diagnostics/core_diagnostics.jl @@ -733,3 +733,160 @@ add_diagnostic_variable!( comments = "Elevation of the horizontal coordinates", compute! = compute_orog!, ) + +### +# Analytic Steady-State Approximations +### + +add_diagnostic_variable!( + short_name = "uapredicted", + long_name = "Predicted Eastward Wind", + standard_name = "predicted_eastward_wind", + units = "m s^-1", + comments = "Predicted value of eastward (zonal) wind component", + compute! = (out, state, cache, time) -> begin + if isnothing(out) + u_component.(cache.predicted_steady_state.ᶜu) + else + out .= u_component.(cache.predicted_steady_state.ᶜu) + end + end, +) + +add_diagnostic_variable!( + short_name = "vapredicted", + long_name = "Predicted Northward Wind", + standard_name = "predicted_northward_wind", + units = "m s^-1", + comments = "Predicted value of northward (meridional) wind component", + compute! = (out, state, cache, time) -> begin + if isnothing(out) + v_component.(cache.predicted_steady_state.ᶜu) + else + out .= v_component.(cache.predicted_steady_state.ᶜu) + end + end, +) + +add_diagnostic_variable!( + short_name = "wapredicted", + long_name = "Predicted Upward Air Velocity", + standard_name = "predicted_upward_air_velocity", + units = "m s^-1", + comments = "Predicted value of vertical wind component", + compute! = (out, state, cache, time) -> begin + if isnothing(out) + w_component.(cache.predicted_steady_state.ᶜu) + else + out .= w_component.(cache.predicted_steady_state.ᶜu) + end + end, +) + +add_diagnostic_variable!( + short_name = "uaerror", + long_name = "Error of Eastward Wind", + standard_name = "error_eastward_wind", + units = "m s^-1", + comments = "Error of eastward (zonal) wind component against predicted value", + compute! = (out, state, cache, time) -> begin + if isnothing(out) + u_component.(Geometry.UVWVector.(cache.precomputed.ᶜu)) .- + u_component.(cache.predicted_steady_state.ᶜu) + else + out .= + u_component.(Geometry.UVWVector.(cache.precomputed.ᶜu)) .- + u_component.(cache.predicted_steady_state.ᶜu) + end + end, +) + +add_diagnostic_variable!( + short_name = "vaerror", + long_name = "Error of Northward Wind", + standard_name = "error_northward_wind", + units = "m s^-1", + comments = "Error of northward (meridional) wind component against predicted value", + compute! = (out, state, cache, time) -> begin + if isnothing(out) + v_component.(Geometry.UVWVector.(cache.precomputed.ᶜu)) .- + v_component.(cache.predicted_steady_state.ᶜu) + else + out .= + v_component.(Geometry.UVWVector.(cache.precomputed.ᶜu)) .- + v_component.(cache.predicted_steady_state.ᶜu) + end + end, +) + +add_diagnostic_variable!( + short_name = "waerror", + long_name = "Error of Upward Air Velocity", + standard_name = "error_upward_air_velocity", + units = "m s^-1", + comments = "Error of vertical wind component against predicted value", + compute! = (out, state, cache, time) -> begin + if isnothing(out) + w_component.(Geometry.UVWVector.(cache.precomputed.ᶜu)) .- + w_component.(cache.predicted_steady_state.ᶜu) + else + out .= + w_component.(Geometry.UVWVector.(cache.precomputed.ᶜu)) .- + w_component.(cache.predicted_steady_state.ᶜu) + end + end, +) + +add_diagnostic_variable!( + short_name = "c1error", + long_name = "Error of First Covariant Wind Component", + standard_name = "error_covariant1_wind", + units = "m^2 s^-1", + comments = "Error of first covariant wind component against predicted value", + compute! = (out, state, cache, time) -> begin + if isnothing(out) + state.c.uₕ.components.data.:1 .- + C12.(cache.predicted_steady_state.ᶜu).components.data.:1 + else + out .= + state.c.uₕ.components.data.:1 .- + C12.(cache.predicted_steady_state.ᶜu).components.data.:1 + end + end, +) + +add_diagnostic_variable!( + short_name = "c2error", + long_name = "Error of Second Covariant Wind Component", + standard_name = "error_covariant2_wind", + units = "m^2 s^-1", + comments = "Error of second covariant wind component against predicted value", + compute! = (out, state, cache, time) -> begin + if isnothing(out) + state.c.uₕ.components.data.:2 .- + C12.(cache.predicted_steady_state.ᶜu).components.data.:2 + else + out .= + state.c.uₕ.components.data.:2 .- + C12.(cache.predicted_steady_state.ᶜu).components.data.:2 + end + end, +) + +add_diagnostic_variable!( + short_name = "c3error", + long_name = "Error of Third Covariant Wind Component", + standard_name = "error_covariant3_wind", + units = "m^2 s^-1", + comments = "Error of third covariant wind component against predicted value", + compute! = (out, state, cache, time) -> begin + if isnothing(out) + state.f.u₃.components.data.:1 .- + C3.(cache.predicted_steady_state.ᶠu).components.data.:1 + else + out .= + state.f.u₃.components.data.:1 .- + C3.(cache.predicted_steady_state.ᶠu).components.data.:1 + end + end, +) diff --git a/src/initial_conditions/InitialConditions.jl b/src/initial_conditions/InitialConditions.jl index 73e0e581518..04380ebd0ea 100644 --- a/src/initial_conditions/InitialConditions.jl +++ b/src/initial_conditions/InitialConditions.jl @@ -17,6 +17,8 @@ import ..n_mass_flux_subdomains import ..gcm_driven_profile import ..gcm_height import ..gcm_driven_profile_tmean +import ..background_u +import ..initial_thermo_state import Thermodynamics.TemperatureProfiles: DecayingTemperatureProfile, DryAdiabaticProfile diff --git a/src/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl index e31f614ebd4..21412ed5041 100644 --- a/src/initial_conditions/initial_conditions.jl +++ b/src/initial_conditions/initial_conditions.jl @@ -121,6 +121,28 @@ function (initial_condition::IsothermalProfile)(params) return local_state end +""" + ConstantBuoyancyFrequencyProfile() + +An `InitialCondition` with a constant Brunt-Vaisala frequency of 0.01 Hz, a +surface temperature of 288 K, an eastward wind velocity of 10 m/s, and a +hydrostatically balanced pressure profile. +""" +struct ConstantBuoyancyFrequencyProfile <: InitialCondition end +function (::ConstantBuoyancyFrequencyProfile)(params) + function local_state(local_geometry) + FT = eltype(params) + coord = local_geometry.coordinates + return LocalState(; + params, + geometry = local_geometry, + thermo_state = initial_thermo_state(params, coord), + velocity = Geometry.UVector(background_u(FT)), + ) + end + return local_state +end + """ DecayingProfile(; perturb = true) diff --git a/src/parameterized_tendencies/sponge/viscous_sponge.jl b/src/parameterized_tendencies/sponge/viscous_sponge.jl index 33abde389de..fc5aef6157c 100644 --- a/src/parameterized_tendencies/sponge/viscous_sponge.jl +++ b/src/parameterized_tendencies/sponge/viscous_sponge.jl @@ -28,14 +28,14 @@ end function viscous_sponge_tendency!(Yₜ, Y, p, t, ::ViscousSponge) (; ᶜβ_viscous, ᶠβ_viscous) = p.viscous_sponge (; ᶜh_tot, ᶜspecific) = p.precomputed - ᶜuₕ = Y.c.uₕ - @. Yₜ.c.uₕ += - ᶜβ_viscous * ( - wgradₕ(divₕ(ᶜuₕ)) - Geometry.project( - Geometry.Covariant12Axis(), - wcurlₕ(Geometry.project(Geometry.Covariant3Axis(), curlₕ(ᶜuₕ))), - ) - ) + point_type = eltype(Fields.coordinate_field(Y.c)) + + if point_type <: Geometry.Abstract3DPoint + @. Yₜ.c.uₕ -= ᶜβ_viscous * C12(wcurlₕ(C3(curlₕ(Y.c.uₕ)))) + end + @. Yₜ.c.uₕ += ᶜβ_viscous * C12(wgradₕ(divₕ(Y.c.uₕ))) + # Without the C12(), the right-hand side would be a C1 or C2 in 2D space. + @. Yₜ.f.u₃.components.data.:1 += ᶠβ_viscous * wdivₕ(gradₕ(Y.f.u₃.components.data.:1)) diff --git a/src/parameters/Parameters.jl b/src/parameters/Parameters.jl index 30b0b5e47e7..8d19b212953 100644 --- a/src/parameters/Parameters.jl +++ b/src/parameters/Parameters.jl @@ -57,7 +57,7 @@ end Base.@kwdef struct ClimaAtmosParameters{ FT, - TP, + TDP, RP, IP, MPC, @@ -66,8 +66,9 @@ Base.@kwdef struct ClimaAtmosParameters{ SFP, TCP, STP, + TP, } <: ACAP - thermodynamics_params::TP + thermodynamics_params::TDP rrtmgp_params::RP insolation_params::IP microphysics_cloud_params::MPC @@ -76,6 +77,7 @@ Base.@kwdef struct ClimaAtmosParameters{ surface_fluxes_params::SFP turbconv_params::TCP surface_temp_params::STP + topography_params::TP Omega::FT f_plane_coriolis_frequency::FT planet_radius::FT diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl index f745791f44c..6fd5804e840 100644 --- a/src/parameters/create_parameters.jl +++ b/src/parameters/create_parameters.jl @@ -61,7 +61,7 @@ function create_parameter_set(config::AtmosConfig) TCP = typeof(turbconv_params) thermodynamics_params = ThermodynamicsParameters(toml_dict) - TP = typeof(thermodynamics_params) + TDP = typeof(thermodynamics_params) rrtmgp_params = RRTMGPParameters(toml_dict) RP = typeof(rrtmgp_params) @@ -79,6 +79,16 @@ function create_parameter_set(config::AtmosConfig) surface_temp_params = SurfaceTemperatureParameters(toml_dict) STP = typeof(surface_temp_params) + topography_params = if parsed_args["topography"] != "NoWarp" + (; + h_max = FT(parsed_args["max_topography_height"]), + z_top = FT(parsed_args["z_max"]), + ) + else + nothing + end + TP = typeof(topography_params) + moisture_model = parsed_args["moist"] microphysics_cloud_params = if moisture_model == "nonequil" (; @@ -149,7 +159,19 @@ function create_parameter_set(config::AtmosConfig) :optics_lookup_temperature_max => :optics_lookup_temperature_max, ) parameters = CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos") - return CAP.ClimaAtmosParameters{FT, TP, RP, IP, MPC, MPP, WP, SFP, TCP, STP}(; + return CAP.ClimaAtmosParameters{ + FT, + TDP, + RP, + IP, + MPC, + MPP, + WP, + SFP, + TCP, + STP, + TP, + }(; parameters..., thermodynamics_params, rrtmgp_params, @@ -160,5 +182,6 @@ function create_parameter_set(config::AtmosConfig) surface_fluxes_params, turbconv_params, surface_temp_params, + topography_params, ) end diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index 44dcbf6ba57..d76d4360579 100644 --- a/src/prognostic_equations/advection.jl +++ b/src/prognostic_equations/advection.jl @@ -39,7 +39,8 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t) @. Yₜ.c.sgs⁰.ρatke -= wdivₕ(Y.c.sgs⁰.ρatke * ᶜu⁰) end - @. Yₜ.c.uₕ -= C12(gradₕ(ᶜp) / Y.c.ρ + gradₕ(ᶜK + ᶜΦ)) + # @. Yₜ.c.uₕ -= C12(gradₕ(ᶜp) / Y.c.ρ + gradₕ(ᶜK + ᶜΦ)) + @. Yₜ.c.uₕ -= C12(ᶜp / Y.c.ρ * gradₕ(log(ᶜp)) + gradₕ(ᶜK + ᶜΦ)) # Without the C12(), the right-hand side would be a C1 or C2 in 2D space. return nothing end diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl index 645ec5d079f..6548508bcdb 100644 --- a/src/prognostic_equations/implicit/implicit_tendency.jl +++ b/src/prognostic_equations/implicit/implicit_tendency.jl @@ -162,7 +162,8 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t) ) end - @. Yₜ.f.u₃ -= ᶠgradᵥ(ᶜp) / ᶠinterp(Y.c.ρ) + ᶠgradᵥ_ᶜΦ + # @. Yₜ.f.u₃ -= ᶠgradᵥ(ᶜp) / ᶠinterp(Y.c.ρ) + ᶠgradᵥ_ᶜΦ + @. Yₜ.f.u₃ -= ᶠinterp(ᶜp / Y.c.ρ) * ᶠgradᵥ(log(ᶜp)) + ᶠgradᵥ_ᶜΦ if rayleigh_sponge isa RayleighSponge (; ᶠβ_rayleigh_w) = p.rayleigh_sponge diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index 3c19e7fb6ac..fde1a45fb2e 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -139,13 +139,25 @@ function get_spaces(parsed_args, params, comms_ctx) bubble = parsed_args["bubble"] deep = parsed_args["deep_atmosphere"] - @assert topography in ("NoWarp", "DCMIP200", "Earth", "Agnesi", "Schar") + @assert topography in ( + "NoWarp", + "DCMIP200", + "Earth", + "Agnesi", + "Schar", + "Cosine2D", + "Cosine3D", + ) if topography == "DCMIP200" warp_function = topography_dcmip200 elseif topography == "Agnesi" warp_function = topography_agnesi elseif topography == "Schar" warp_function = topography_schar + elseif topography == "Cosine2D" + warp_function = topography_cosine_2d + elseif topography == "Cosine3D" + warp_function = topography_cosine_3d elseif topography == "NoWarp" warp_function = nothing elseif topography == "Earth" @@ -194,6 +206,7 @@ function get_spaces(parsed_args, params, comms_ctx) z_max, z_elem, z_stretch; + params, parsed_args = parsed_args, surface_warp = warp_function, deep, @@ -248,6 +261,7 @@ function get_spaces(parsed_args, params, comms_ctx) z_max, z_elem, z_stretch; + params, parsed_args, surface_warp = warp_function, deep, @@ -272,6 +286,7 @@ function get_spaces(parsed_args, params, comms_ctx) z_max, z_elem, z_stretch; + params, parsed_args, surface_warp = warp_function, deep, @@ -345,6 +360,7 @@ function get_initial_condition(parsed_args) ) elseif parsed_args["initial_condition"] in [ "IsothermalProfile", + "ConstantBuoyancyFrequencyProfile", "AgnesiHProfile", "DryDensityCurrentProfile", "RisingThermalBubbleProfile", @@ -369,6 +385,31 @@ function get_surface_setup(parsed_args) return getproperty(SurfaceConditions, Symbol(parsed_args["surface_setup"]))() end +function get_predicted_steady_state(params, Y, parsed_args) + parsed_args["analytic_check"] || return nothing + @assert parsed_args["initial_condition"] == + "ConstantBuoyancyFrequencyProfile" + topography = parsed_args["topography"] + predicted_velocity = if topography == "Agnesi" + predicted_velocity_agnesi + elseif topography == "Schar" + predicted_velocity_schar + elseif topography == "Cosine2D" + predicted_velocity_cosine_2d + elseif topography == "Cosine3D" + predicted_velocity_cosine_3d + else + error("Analytic solution for `topography = $topography` is not defined") + end + @info "Computing linearized approximation of steady-state solution" + s = @timed_str begin + ᶜu = predicted_velocity.(params, Fields.coordinate_field(Y.c)) + ᶠu = predicted_velocity.(params, Fields.coordinate_field(Y.f)) + end + @info "Steady-state approximation completed: $s" + return (; ᶜu, ᶠu) +end + is_explicit_CTS_algo_type(alg_or_tableau) = alg_or_tableau <: CTS.ERKAlgorithmName @@ -643,6 +684,8 @@ function get_simulation(config::AtmosConfig) @info "Allocating Y: $s" end + predicted_steady_state = + get_predicted_steady_state(params, Y, config.parsed_args) tracers = get_tracers(config.parsed_args) s = @timed_str begin @@ -653,6 +696,7 @@ function get_simulation(config::AtmosConfig) surface_setup, sim_info, tracers.aerosol_names, + predicted_steady_state, ) end @info "Allocating cache (p): $s" diff --git a/src/topography/topography.jl b/src/topography/topography.jl index 57547749803..06e633ca4af 100644 --- a/src/topography/topography.jl +++ b/src/topography/topography.jl @@ -6,74 +6,33 @@ Given horizontal coordinates in lon-lat coordinates, compute and return the local elevation of the surface consistent with the test problem DCMIP-2-0-0. """ -function topography_dcmip200(coords) - λ, ϕ = coords.long, coords.lat +function topography_dcmip200(_, coord) + λ, ϕ = coord.long, coord.lat FT = eltype(λ) ϕₘ = FT(0) # degrees (equator) - λₘ = FT(3 / 2 * 180) # degrees - rₘ = @. FT(acos(sind(ϕₘ) * sind(ϕ) + cosd(ϕₘ) * cosd(ϕ) * cosd(λ - λₘ))) # Great circle distance (rads) - Rₘ = FT(3π / 4) # Moutain radius - ζₘ = FT(π / 16) # Mountain oscillation half-width + λₘ = 3 * FT(180) / 2 # degrees + rₘ = acos(sind(ϕₘ) * sind(ϕ) + cosd(ϕₘ) * cosd(ϕ) * cosd(λ - λₘ)) # Great circle distance (rads) + Rₘ = 3 * FT(π) / 4 # Moutain radius + ζₘ = FT(π) / 16 # Mountain oscillation half-width h₀ = FT(2000) - zₛ = @. ifelse( + zₛ = ifelse( rₘ < Rₘ, - FT(h₀ / 2) * (1 + cospi(rₘ / Rₘ)) * (cospi(rₘ / ζₘ))^2, + h₀ / 2 * (1 + cospi(rₘ / Rₘ)) * (cospi(rₘ / ζₘ))^2, FT(0), ) return zₛ end function generate_topography_warp(earth_spline) - function topography_earth(coords) - λ, Φ = coords.long, coords.lat + function topography_earth(_, coord) + λ, Φ = coord.long, coord.lat FT = eltype(λ) - @info "Load spline" - elevation = @. FT(earth_spline(λ, Φ)) - zₛ = @. ifelse(elevation > FT(0), elevation, FT(0)) - @info "Assign elevation" + elevation = FT(earth_spline(λ, Φ)) + zₛ = ifelse(elevation > 0, elevation, FT(0)) return zₛ end return topography_earth end -""" - topography_agnesi(x,z) -x = horizontal coordinate [m] -z = vertical coordinate [m] -h_c = 1 [m] -a_c = 10000 [m] -x_c = 120000 [m] -Generate a single mountain profile (Agnesi mountain) -for use with tests of gravity waves with topography. -""" -function topography_agnesi(coords) - x = coords.x - FT = eltype(x) - h_c = FT(1) - a_c = FT(10000) - x_c = FT(120000) - zₛ = @. h_c / (1 + ((x - x_c) / a_c)^2) - return zₛ -end - -""" - topography_schar(x,z) -x = horizontal coordinate [m] -z = vertical coordinate [m] -h_c = 250 [m] -a_c = 5000 [m] -x_c = 60000 [m] -Assumes [0, 120] km domain. -Generate a single mountain profile (Schar mountain) -for use with tests of gravity waves with topography. -""" -function topography_schar(coords) - x = coords.x - FT = eltype(x) - h_c = FT(250) - λ_c = FT(4000) - a_c = FT(5000) - x_c = FT(60000) - zₛ = @. h_c * exp(-((x - x_c) / a_c)^2) * (cospi((x - x_c) / λ_c))^2 - return zₛ -end +# Additional topography functions with analytic Fourier transforms can be found +# in src/analytic_solutions/analytic_solutions.jl. diff --git a/src/utils/common_spaces.jl b/src/utils/common_spaces.jl index 7feb5f38059..924347d3fa8 100644 --- a/src/utils/common_spaces.jl +++ b/src/utils/common_spaces.jl @@ -79,6 +79,7 @@ function make_hybrid_spaces( z_max, z_elem, z_stretch; + params = nothing, surface_warp = nothing, topo_smoothing = false, deep = false, @@ -104,7 +105,7 @@ function make_hybrid_spaces( hypsography = Hypsography.Flat() else topo_smoothing = parsed_args["topo_smoothing"] - z_surface = surface_warp(Fields.coordinate_field(h_space)) + z_surface = surface_warp.(params, Fields.coordinate_field(h_space)) if topo_smoothing Hypsography.diffuse_surface_elevation!(z_surface) end diff --git a/toml/plane_agnesi_mountain_test_stretched.toml b/toml/analytic_topography_high_sponge_test.toml similarity index 53% rename from toml/plane_agnesi_mountain_test_stretched.toml rename to toml/analytic_topography_high_sponge_test.toml index 7852943cccf..8dbb23cc7fc 100644 --- a/toml/plane_agnesi_mountain_test_stretched.toml +++ b/toml/analytic_topography_high_sponge_test.toml @@ -1,3 +1,5 @@ [alpha_rayleigh_w] value = 0.1 +[zd_rayleigh] +value = 21000 diff --git a/toml/analytic_topography_strong_sponge_test.toml b/toml/analytic_topography_strong_sponge_test.toml new file mode 100644 index 00000000000..942853c012a --- /dev/null +++ b/toml/analytic_topography_strong_sponge_test.toml @@ -0,0 +1,5 @@ +[alpha_rayleigh_w] +value = 0.2 + +[zd_rayleigh] +value = 13000 diff --git a/toml/plane_agnesi_mountain_test_uniform.toml b/toml/analytic_topography_test.toml similarity index 53% rename from toml/plane_agnesi_mountain_test_uniform.toml rename to toml/analytic_topography_test.toml index 7852943cccf..da5ffcc6203 100644 --- a/toml/plane_agnesi_mountain_test_uniform.toml +++ b/toml/analytic_topography_test.toml @@ -1,3 +1,5 @@ [alpha_rayleigh_w] value = 0.1 +[zd_rayleigh] +value = 13000 diff --git a/toml/analytic_topography_weak_high_sponge_test.toml b/toml/analytic_topography_weak_high_sponge_test.toml new file mode 100644 index 00000000000..213f56bc65b --- /dev/null +++ b/toml/analytic_topography_weak_high_sponge_test.toml @@ -0,0 +1,5 @@ +[alpha_rayleigh_w] +value = 0.05 + +[zd_rayleigh] +value = 21000 diff --git a/toml/analytic_topography_weak_sponge_test.toml b/toml/analytic_topography_weak_sponge_test.toml new file mode 100644 index 00000000000..373ccaf2a3e --- /dev/null +++ b/toml/analytic_topography_weak_sponge_test.toml @@ -0,0 +1,5 @@ +[alpha_rayleigh_w] +value = 0.05 + +[zd_rayleigh] +value = 13000 diff --git a/toml/plane_schar_mountain_test_stretched.toml b/toml/plane_schar_mountain_test_stretched.toml deleted file mode 100644 index 7852943cccf..00000000000 --- a/toml/plane_schar_mountain_test_stretched.toml +++ /dev/null @@ -1,3 +0,0 @@ -[alpha_rayleigh_w] -value = 0.1 - diff --git a/toml/plane_schar_mountain_test_uniform.toml b/toml/plane_schar_mountain_test_uniform.toml deleted file mode 100644 index 7852943cccf..00000000000 --- a/toml/plane_schar_mountain_test_uniform.toml +++ /dev/null @@ -1,3 +0,0 @@ -[alpha_rayleigh_w] -value = 0.1 -