From dbf85f79db4177cb7d8428f718735d88d3d4c093 Mon Sep 17 00:00:00 2001 From: "Katherine M. Deck" Date: Thu, 5 Sep 2024 14:10:52 -0700 Subject: [PATCH] Regional land simulation --- .buildkite/Manifest.toml | 101 +-- .buildkite/longruns_gpu/pipeline.yml | 10 + NEWS.md | 1 + Project.toml | 2 +- experiments/long_runs/land_region.jl | 676 +++++++++++++++++++ experiments/standalone/Bucket/bucket_era5.jl | 4 +- src/shared_utilities/Domains.jl | 29 +- test/shared_utilities/domains.jl | 12 +- 8 files changed, 773 insertions(+), 62 deletions(-) create mode 100644 experiments/long_runs/land_region.jl diff --git a/.buildkite/Manifest.toml b/.buildkite/Manifest.toml index f5fd82a1ab..3f8e546bac 100644 --- a/.buildkite/Manifest.toml +++ b/.buildkite/Manifest.toml @@ -101,9 +101,9 @@ version = "0.4.0" [[deps.ArrayInterface]] deps = ["Adapt", "LinearAlgebra"] -git-tree-sha1 = "f54c23a5d304fb87110de62bace7777d59088c34" +git-tree-sha1 = "3640d077b6dafd64ceb8fd5c1ec76f7ca53bcf76" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "7.15.0" +version = "7.16.0" [deps.ArrayInterface.extensions] ArrayInterfaceBandedMatricesExt = "BandedMatrices" @@ -390,34 +390,35 @@ version = "0.4.5" [[deps.ClimaDiagnostics]] deps = ["Accessors", "ClimaComms", "ClimaCore", "Dates", "NCDatasets", "SciMLBase"] -git-tree-sha1 = "4f8abbf3af5a78b36bc40a33ef7e3fa1d3e8f138" +git-tree-sha1 = "0c60c2b40b0c5ae87690734e2aedea5a8c686e91" uuid = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f" -version = "0.2.4" +version = "0.2.5" [[deps.ClimaLand]] -deps = ["ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaUtilities", "DataFrames", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "IntervalSets", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics"] +deps = ["ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaUtilities", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics"] path = ".." uuid = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532" version = "0.14.3" [deps.ClimaLand.extensions] CreateParametersExt = "ClimaParams" - NeuralSnowExt = ["CSV", "HTTP", "Flux", "StatsBase", "cuDNN"] + NeuralSnowExt = ["CSV", "DataFrames", "HTTP", "Flux", "StatsBase", "cuDNN"] [deps.ClimaLand.weakdeps] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" + DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" cuDNN = "02a925ec-e4fe-4b08-9a7e-0d78e3d38ccd" [[deps.ClimaParams]] -deps = ["DocStringExtensions", "TOML", "Test"] -git-tree-sha1 = "56f4b28affba295c0cbc81d9684ddc759af04bfc" +deps = ["TOML"] +git-tree-sha1 = "b43ca371c435056129295445122ea87fd843b505" uuid = "5c42b081-d73a-476f-9059-fd94b934656c" -version = "0.10.13" +version = "0.10.14" [[deps.ClimaTimeSteppers]] deps = ["ClimaComms", "Colors", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "KernelAbstractions", "Krylov", "LinearAlgebra", "LinearOperators", "NVTX", "SciMLBase", "StaticArrays"] @@ -558,14 +559,14 @@ uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb" version = "2.4.2" [[deps.ConstructionBase]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "a33b7ced222c6165f624a3f2b55945fac5a598d9" +git-tree-sha1 = "76219f1ed5771adbb096743bff43fb5fdd4c1157" uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" -version = "1.5.7" -weakdeps = ["IntervalSets", "StaticArrays"] +version = "1.5.8" +weakdeps = ["IntervalSets", "LinearAlgebra", "StaticArrays"] [deps.ConstructionBase.extensions] ConstructionBaseIntervalSetsExt = "IntervalSets" + ConstructionBaseLinearAlgebraExt = "LinearAlgebra" ConstructionBaseStaticArraysExt = "StaticArrays" [[deps.ContextVariablesX]] @@ -597,10 +598,10 @@ uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" version = "4.1.1" [[deps.CubedSphere]] -deps = ["Elliptic", "FFTW", "Printf", "ProgressBars", "SpecialFunctions", "TaylorSeries", "Test"] -git-tree-sha1 = "10134667d7d3569b191a65801514271b8a93b292" +deps = ["Elliptic", "FFTW", "Printf", "ProgressBars", "SpecialFunctions", "TaylorSeries"] +git-tree-sha1 = "049f692019f52ad3b043edf20ef25ebd7b38eb94" uuid = "7445602f-e544-4518-8976-18f8e8ae6cdb" -version = "0.2.5" +version = "0.2.6" [[deps.DataAPI]] git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe" @@ -641,9 +642,9 @@ version = "0.1.2" [[deps.DelaunayTriangulation]] deps = ["AdaptivePredicates", "EnumX", "ExactPredicates", "Random"] -git-tree-sha1 = "b5f1c6532d2ea71e99b74231b0a3d53fba846ced" +git-tree-sha1 = "9903123ab7fc5e55053292aff04ff5d7aff92633" uuid = "927a84f5-c5f4-47a5-9785-b46e178433df" -version = "1.1.3" +version = "1.3.0" [[deps.DelimitedFiles]] deps = ["Mmap"] @@ -707,9 +708,9 @@ version = "1.15.1" [[deps.DifferentiationInterface]] deps = ["ADTypes", "Compat", "DocStringExtensions", "FillArrays", "LinearAlgebra", "PackageExtensionCompat", "SparseArrays", "SparseMatrixColorings"] -git-tree-sha1 = "e53e4eb7f8e7b25656ea17263b0c22f129091eeb" +git-tree-sha1 = "9b23f9a816790b8ab9914c3c86321a546e92cbe7" uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" -version = "0.5.15" +version = "0.5.17" [deps.DifferentiationInterface.extensions] DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore" @@ -1118,10 +1119,10 @@ uuid = "68eda718-8dee-11e9-39e7-89f7f65f511f" version = "0.4.2" [[deps.GeoInterface]] -deps = ["Extents"] -git-tree-sha1 = "9fff8990361d5127b770e3454488360443019bb3" +deps = ["Extents", "GeoFormatTypes"] +git-tree-sha1 = "5921fc0704e40c024571eca551800c699f86ceb4" uuid = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" -version = "1.3.5" +version = "1.3.6" [[deps.GeoInterfaceMakie]] deps = ["GeoInterface", "GeometryBasics", "MakieCore"] @@ -1590,9 +1591,9 @@ version = "0.1.17" [[deps.LazyArrays]] deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "SparseArrays"] -git-tree-sha1 = "507b423197fdd9e77b74aa2532c0a05eb7eb4004" +git-tree-sha1 = "360f6039babd6e4d6364eff0d4fc9120834a2d9a" uuid = "5078a376-72f3-5289-bfd5-ec5146d43c02" -version = "2.2.0" +version = "2.2.1" [deps.LazyArrays.extensions] LazyArraysBandedMatricesExt = "BandedMatrices" @@ -1717,10 +1718,10 @@ version = "2.8.0" LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" [[deps.LinearSolve]] -deps = ["ArrayInterface", "ChainRulesCore", "ConcreteStructs", "CpuId", "DocStringExtensions", "EnumX", "FastLapackInterface", "GPUArraysCore", "InteractiveUtils", "KLU", "Krylov", "LazyArrays", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "PrecompileTools", "Preferences", "RecursiveFactorization", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Sparspak", "StaticArraysCore", "UnPack"] -git-tree-sha1 = "6b7011b27e35453f68a8f851776271b31fce021c" +deps = ["ArrayInterface", "ChainRulesCore", "ConcreteStructs", "DocStringExtensions", "EnumX", "FastLapackInterface", "GPUArraysCore", "InteractiveUtils", "KLU", "Krylov", "LazyArrays", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "PrecompileTools", "Preferences", "RecursiveFactorization", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Sparspak", "StaticArraysCore", "UnPack"] +git-tree-sha1 = "6c5e4555ac2bc449a28604e184f759d18fc08420" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -version = "2.33.0" +version = "2.34.0" [deps.LinearSolve.extensions] LinearSolveBandedMatricesExt = "BandedMatrices" @@ -1966,10 +1967,10 @@ uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" version = "4.5.1" [[deps.NNlib]] -deps = ["Adapt", "Atomix", "ChainRulesCore", "GPUArraysCore", "KernelAbstractions", "LinearAlgebra", "Pkg", "Random", "Requires", "Statistics"] -git-tree-sha1 = "ae52c156a63bb647f80c26319b104e99e5977e51" +deps = ["Adapt", "Atomix", "ChainRulesCore", "GPUArraysCore", "KernelAbstractions", "LinearAlgebra", "Random", "Statistics"] +git-tree-sha1 = "4a83c2e01027a0bfcea28589222f2df60b2e20cb" uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" -version = "0.9.22" +version = "0.9.23" [deps.NNlib.extensions] NNlibAMDGPUExt = "AMDGPU" @@ -1977,12 +1978,14 @@ version = "0.9.22" NNlibCUDAExt = "CUDA" NNlibEnzymeCoreExt = "EnzymeCore" NNlibFFTWExt = "FFTW" + NNlibForwardDiffExt = "ForwardDiff" [deps.NNlib.weakdeps] AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" cuDNN = "02a925ec-e4fe-4b08-9a7e-0d78e3d38ccd" [[deps.NVTX]] @@ -2380,9 +2383,15 @@ version = "6.5.2+2" [[deps.QuadGK]] deps = ["DataStructures", "LinearAlgebra"] -git-tree-sha1 = "e237232771fdafbae3db5c31275303e056afaa9f" +git-tree-sha1 = "1d587203cf851a51bf1ea31ad7ff89eff8d625ea" uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" -version = "2.10.1" +version = "2.11.0" + + [deps.QuadGK.extensions] + QuadGKEnzymeExt = "Enzyme" + + [deps.QuadGK.weakdeps] + Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" [[deps.REPL]] deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] @@ -2549,9 +2558,9 @@ version = "0.1.0" [[deps.SciMLBase]] deps = ["ADTypes", "Accessors", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "Expronicon", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"] -git-tree-sha1 = "2c98f2723eea6a86fd4e956703f161ee1e6f0c3f" +git-tree-sha1 = "303a73db99326a8be43e695fbab9e076b02118ca" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -version = "2.50.6" +version = "2.52.2" [deps.SciMLBase.extensions] SciMLBaseChainRulesCoreExt = "ChainRulesCore" @@ -2585,9 +2594,9 @@ weakdeps = ["SparseArrays", "StaticArraysCore"] [[deps.SciMLStructures]] deps = ["ArrayInterface"] -git-tree-sha1 = "20ad3e7c137156c50c93c888d0f2bc5b7883c729" +git-tree-sha1 = "25514a6f200219cd1073e4ff23a6324e4a7efe64" uuid = "53ae85a6-f571-4167-b2af-e1d143709226" -version = "1.4.2" +version = "1.5.0" [[deps.Scratch]] deps = ["Dates"] @@ -2864,9 +2873,9 @@ weakdeps = ["ClimaParams"] [[deps.SymbolicIndexingInterface]] deps = ["Accessors", "ArrayInterface", "RuntimeGeneratedFunctions", "StaticArraysCore"] -git-tree-sha1 = "c9fce29fb41a10677e24f74421ebe31220b81ad0" +git-tree-sha1 = "988e04b34a4c3b824fb656f542473df99a4f610d" uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5" -version = "0.3.28" +version = "0.3.30" [[deps.TOML]] deps = ["Dates"] @@ -2892,13 +2901,21 @@ version = "1.10.0" [[deps.TaylorSeries]] deps = ["LinearAlgebra", "Markdown", "Requires", "SparseArrays"] -git-tree-sha1 = "1c7170668366821b0c4c4fe03ee78f8d6cf36e2c" +git-tree-sha1 = "90c9bc500f4c5cdd235c81503ec91b2048f06423" uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" -version = "0.16.0" -weakdeps = ["IntervalArithmetic"] +version = "0.17.8" [deps.TaylorSeries.extensions] TaylorSeriesIAExt = "IntervalArithmetic" + TaylorSeriesJLD2Ext = "JLD2" + TaylorSeriesRATExt = "RecursiveArrayTools" + TaylorSeriesSAExt = "StaticArrays" + + [deps.TaylorSeries.weakdeps] + IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" + JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" + RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [[deps.TensorCore]] deps = ["LinearAlgebra"] diff --git a/.buildkite/longruns_gpu/pipeline.yml b/.buildkite/longruns_gpu/pipeline.yml index e57a0641ed..84331200dd 100644 --- a/.buildkite/longruns_gpu/pipeline.yml +++ b/.buildkite/longruns_gpu/pipeline.yml @@ -41,6 +41,16 @@ steps: env: CLIMACOMMS_DEVICE: "CUDA" + - label: ":sun: California regional simulation" + command: + - julia --color=yes --project=.buildkite experiments/long_runs/land_region.jl + artifact_paths: "california_longrun_gpu/*png" + agents: + slurm_gpus: 1 + slurm_time: 03:00:00 + env: + CLIMACOMMS_DEVICE: "CUDA" + - label: "Soil" command: - julia --color=yes --project=.buildkite experiments/long_runs/soil.jl diff --git a/NEWS.md b/NEWS.md index 8fc556f3e3..2e40365078 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,7 @@ ClimaLand.jl Release Notes main -------- + - Add regional simulation example PR [#757](https://github.com/CliMA/ClimaLand.jl/pull/757) v0.14.3 -------- diff --git a/Project.toml b/Project.toml index 234a87e318..106b983cff 100644 --- a/Project.toml +++ b/Project.toml @@ -41,7 +41,7 @@ CSV = "0.10" CUDA = "5.2" ClimaComms = "0.6" ClimaCore = "0.14" -ClimaDiagnostics = "0.2" +ClimaDiagnostics = "0.2.5" ClimaParams = "0.10.2" ClimaUtilities = "0.1.6" DataFrames = "1.4" diff --git a/experiments/long_runs/land_region.jl b/experiments/long_runs/land_region.jl new file mode 100644 index 0000000000..35c1eceb81 --- /dev/null +++ b/experiments/long_runs/land_region.jl @@ -0,0 +1,676 @@ +# # Global run of land model + +# The code sets up and runs the soil/canopy model for 6 hours on a small region of the globe in Southern California, +# using ERA5 data. In this simulation, we have +# turned lateral flow off because horizontal boundary conditions and the +# land/sea mask are not yet supported by ClimaCore. + +# Simulation Setup +# Number of spatial elements: 10x10 in horizontal, 15 in vertical +# Soil depth: 50 m +# Simulation duration: 334 d +# Timestep: 900 s +# Timestepper: ARS111 +# Fixed number of iterations: 1 +# Jacobian update: every new timestep +# Atmos forcing update: every 3 hours +import SciMLBase +import ClimaComms +ClimaComms.@import_required_backends +import ClimaTimeSteppers as CTS +using ClimaCore +using ClimaUtilities.ClimaArtifacts +import Interpolations +using Insolation + +import ClimaDiagnostics +import ClimaAnalysis +import ClimaAnalysis.Visualize as viz +import ClimaUtilities + +import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput +import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput +import ClimaUtilities.Regridders: InterpolationsRegridder +import ClimaUtilities.ClimaArtifacts: @clima_artifact +import ClimaParams as CP + +using ClimaLand +using ClimaLand.Soil +using ClimaLand.Canopy +import ClimaLand +import ClimaLand.Parameters as LP + +using Statistics +using CairoMakie +using Dates +import NCDatasets + +const FT = Float64; + +regridder_type = :InterpolationsRegridder +context = ClimaComms.context() +device = ClimaComms.device() +device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu" +root_path = "california_longrun_$(device_suffix)" +diagnostics_outdir = joinpath(root_path, "regional_diagnostics") +outdir = + ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir) + +function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15)) + + earth_param_set = LP.LandParameters(FT) + radius = FT(6378.1e3) + depth = FT(50) + center_long, center_lat = FT(-117.59736), FT(34.23375) + delta_m = FT(100_000) # in meters + domain = ClimaLand.Domains.HybridBox(; + xlim = (delta_m, delta_m), + ylim = (delta_m, delta_m), + zlim = (-depth, FT(0)), + nelements = nelements, + npolynomial = 1, + longlat = (center_long, center_lat), + dz_tuple = FT.((10.0, 0.5)), + ) + surface_space = domain.space.surface + subsurface_space = domain.space.subsurface + + ref_time = DateTime(2021) + t_start = t0 + # Forcing data + era5_artifact_path = + ClimaLand.Artifacts.era5_land_forcing_data2021_folder_path(; context) # Precipitation: + precip = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "rf", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,), + ) + + snow_precip = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "sf", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,), + ) + + u_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "ws", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + ) + q_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "q", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + ) + P_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "sp", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + ) + + T_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "t2m", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + ) + h_atmos = FT(10) + + atmos = PrescribedAtmosphere( + precip, + snow_precip, + T_atmos, + u_atmos, + q_atmos, + P_atmos, + ref_time, + h_atmos, + earth_param_set, + ) + + # Prescribed radiation + SW_d = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "ssrd", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,), + ) + LW_d = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "strd", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,), + ) + + function zenith_angle( + t, + ref_time; + latitude = ClimaCore.Fields.coordinate_field(surface_space).lat, + longitude = ClimaCore.Fields.coordinate_field(surface_space).long, + insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params, + ) where {FT} + # This should be time in UTC + current_datetime = ref_time + Dates.Second(round(t)) + + # Orbital Data uses Float64, so we need to convert to our sim FT + d, δ, η_UTC = + FT.( + Insolation.helper_instantaneous_zenith_angle( + current_datetime, + ref_time, + insol_params, + ) + ) + + Insolation.instantaneous_zenith_angle.( + d, + δ, + η_UTC, + longitude, + latitude, + ).:1 + end + radiation = + PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time; θs = zenith_angle) + + soil_params_artifact_path = + ClimaLand.Artifacts.soil_params_artifact_folder_path(; context) + extrapolation_bc = ( + Interpolations.Periodic(), + Interpolations.Flat(), + Interpolations.Flat(), + ) + soil_params_mask = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc", + ), + "α", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + file_reader_kwargs = (; preprocess_func = (data) -> data > 0,), + ) + oceans_to_value(field, mask, value) = + mask == 1.0 ? field : eltype(field)(value) + + vg_α = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "vGalpha_map_gupta_etal2020_1.0x1.0x4.nc", + ), + "α", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + vg_n = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "vGn_map_gupta_etal2020_1.0x1.0x4.nc", + ), + "n", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + x = parent(vg_α) + μ = mean(log10.(x[x .> 0])) + vg_α .= oceans_to_value.(vg_α, soil_params_mask, 10.0^μ) + + x = parent(vg_n) + μ = mean(x[x .> 0]) + vg_n .= oceans_to_value.(vg_n, soil_params_mask, μ) + + vg_fields_to_hcm_field(α::FT, n::FT) where {FT} = + ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...) + hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n) + + θ_r = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "residual_map_gupta_etal2020_1.0x1.0x4.nc", + ), + "θ_r", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + + ν = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "porosity_map_gupta_etal2020_1.0x1.0x4.nc", + ), + "ν", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + K_sat = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "ksat_map_gupta_etal2020_1.0x1.0x4.nc", + ), + "Ksat", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + + x = parent(K_sat) + μ = mean(log10.(x[x .> 0])) + K_sat .= oceans_to_value.(K_sat, soil_params_mask, 10.0^μ) + + ν .= oceans_to_value.(ν, soil_params_mask, 1) + + θ_r .= oceans_to_value.(θ_r, soil_params_mask, 0) + + + S_s = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3), + soil_params_mask, + 1, + ) + ν_ss_om = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1), + soil_params_mask, + 0, + ) + ν_ss_quartz = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1), + soil_params_mask, + 0, + ) + ν_ss_gravel = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1), + soil_params_mask, + 0, + ) + PAR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2) + NIR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2) + soil_params = Soil.EnergyHydrologyParameters( + FT; + ν, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, + hydrology_cm, + K_sat, + S_s, + θ_r, + PAR_albedo = PAR_albedo, + NIR_albedo = NIR_albedo, + ) + + soil_params_mask_sfc = + ClimaLand.Domains.top_center_to_surface(soil_params_mask) + + # Read in f_max data and land sea mask + infile_path = ClimaLand.Artifacts.topmodel_data_path() + f_max = + SpaceVaryingInput(infile_path, "fmax", surface_space; regridder_type) + mask = SpaceVaryingInput( + infile_path, + "landsea_mask", + surface_space; + regridder_type, + ) + # Unsure how to handle two masks + f_max = oceans_to_value.(f_max, mask, FT(0.0)) + f_max = oceans_to_value.(f_max, soil_params_mask_sfc, FT(0.0)) + f_over = FT(3.28) # 1/m + R_sb = FT(1.484e-4 / 1000) # m/s + runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(; + f_over = f_over, + f_max = f_max, + R_sb = R_sb, + ) + + + # TwoStreamModel parameters + Ω = FT(0.69) + ld = FT(0.5) + α_PAR_leaf = FT(0.1) + τ_PAR_leaf = FT(0.05) + α_NIR_leaf = FT(0.45) + τ_NIR_leaf = FT(0.25) + + # Energy Balance model + ac_canopy = FT(2.5e4) # this will likely be 10x smaller! + + # Conductance Model + g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. + + clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context) + # vcmax is read in units of umol CO2/m^2/s and then converted to mol CO2/m^2/s + Vcmax25 = SpaceVaryingInput( + joinpath(clm_artifact_path, "vegetation_properties_map.nc"), + "vcmx25", + surface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + file_reader_kwargs = (; preprocess_func = (data) -> data / 1_000_000,), + ) + + # Plant Hydraulics and general plant parameters + SAI = FT(0.0) # m2/m2 + f_root_to_shoot = FT(3.5) + RAI = FT(1.0) + K_sat_plant = FT(5e-9) # m/s # seems much too small? + ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value is -4 MPa + Weibull_param = FT(4) # unitless, Holtzman's original c param value + a = FT(0.05 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity + conductivity_model = + Canopy.PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) + retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a) + plant_ν = FT(1.44e-4) + plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m + rooting_depth = FT(0.5) # from Natan + n_stem = 0 + n_leaf = 1 + h_stem = FT(0.0) + h_leaf = FT(1.0) + zmax = FT(0.0) + h_canopy = h_stem + h_leaf + compartment_midpoints = + n_stem > 0 ? [h_stem / 2, h_stem + h_leaf / 2] : [h_leaf / 2] + compartment_surfaces = + n_stem > 0 ? [zmax, h_stem, h_canopy] : [zmax, h_leaf] + + z0_m = FT(0.13) * h_canopy + z0_b = FT(0.1) * z0_m + + + soilco2_ps = Soil.Biogeochemistry.SoilCO2ModelParameters(FT) + + soil_args = (domain = domain, parameters = soil_params) + soil_model_type = Soil.EnergyHydrology{FT} + + # Soil microbes model + soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT} + + # soil microbes args + Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5)) + + # Set the soil CO2 BC to being atmospheric CO2 + soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC() + soilco2_bot_bc = Soil.Biogeochemistry.SoilCO2FluxBC((p, t) -> 0.0) # no flux + soilco2_sources = (Soil.Biogeochemistry.MicrobeProduction{FT}(),) + + soilco2_boundary_conditions = + (; top = soilco2_top_bc, bottom = soilco2_bot_bc) + + soilco2_args = (; + boundary_conditions = soilco2_boundary_conditions, + sources = soilco2_sources, + domain = domain, + parameters = soilco2_ps, + ) + + # Now we set up the canopy model, which we set up by component: + # Component Types + canopy_component_types = (; + autotrophic_respiration = Canopy.AutotrophicRespirationModel{FT}, + radiative_transfer = Canopy.TwoStreamModel{FT}, + photosynthesis = Canopy.FarquharModel{FT}, + conductance = Canopy.MedlynConductanceModel{FT}, + hydraulics = Canopy.PlantHydraulicsModel{FT}, + energy = Canopy.BigLeafEnergyModel{FT}, + ) + # Individual Component arguments + # Set up autotrophic respiration + autotrophic_respiration_args = + (; parameters = Canopy.AutotrophicRespirationParameters(FT)) + # Set up radiative transfer + radiative_transfer_args = (; + parameters = Canopy.TwoStreamParameters( + FT; + Ω, + α_PAR_leaf, + τ_PAR_leaf, + α_NIR_leaf, + τ_NIR_leaf, + ) + ) + # Set up conductance + conductance_args = + (; parameters = Canopy.MedlynConductanceParameters(FT; g1)) + # Set up photosynthesis + photosynthesis_args = (; + parameters = Canopy.FarquharParameters( + FT, + Canopy.C3(); + Vcmax25 = Vcmax25, + ) + ) + # Set up plant hydraulics + + # Note that we clip all values of LAI below 0.05 to zero. + # This is because we currently run into issues when LAI is + # of order eps(FT) in the SW radiation code. + # Please see Issue #644 + # or PR #645 for details. + # For now, this clipping is similar to what CLM does. + LAIfunction = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_lai_2021_0.9x1.25_clima.nc"), + "lai", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; + preprocess_func = (data) -> data > 0.05 ? data : 0.0, + ), + ) + ai_parameterization = + Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI) + + function root_distribution(z::T; rooting_depth = rooting_depth) where {T} + return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m + end + + plant_hydraulics_ps = Canopy.PlantHydraulics.PlantHydraulicsParameters(; + ai_parameterization = ai_parameterization, + ν = plant_ν, + S_s = plant_S_s, + root_distribution = root_distribution, + conductivity_model = conductivity_model, + retention_model = retention_model, + ) + plant_hydraulics_args = ( + parameters = plant_hydraulics_ps, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_midpoints = compartment_midpoints, + compartment_surfaces = compartment_surfaces, + ) + + energy_args = (parameters = Canopy.BigLeafEnergyParameters{FT}(ac_canopy),) + + # Canopy component args + canopy_component_args = (; + autotrophic_respiration = autotrophic_respiration_args, + radiative_transfer = radiative_transfer_args, + photosynthesis = photosynthesis_args, + conductance = conductance_args, + hydraulics = plant_hydraulics_args, + energy = energy_args, + ) + + # Other info needed + shared_params = Canopy.SharedCanopyParameters{FT, typeof(earth_param_set)}( + z0_m, + z0_b, + earth_param_set, + ) + + canopy_model_args = (; + parameters = shared_params, + domain = ClimaLand.obtain_surface_domain(domain), + ) + + # Integrated plant hydraulics and soil model + land_input = ( + atmos = atmos, + radiation = radiation, + runoff = runoff_model, + soil_organic_carbon = Csom, + ) + land = SoilCanopyModel{FT}(; + soilco2_type = soilco2_type, + soilco2_args = soilco2_args, + land_args = land_input, + soil_model_type = soil_model_type, + soil_args = soil_args, + canopy_component_types = canopy_component_types, + canopy_component_args = canopy_component_args, + canopy_model_args = canopy_model_args, + ) + + Y, p, cds = initialize(land) + + init_soil(ν, θ_r) = θ_r + (ν - θ_r) / 2 + Y.soil.ϑ_l .= init_soil.(ν, θ_r) + Y.soil.θ_i .= FT(0.0) + T = FT(276.85) + ρc_s = + Soil.volumetric_heat_capacity.( + Y.soil.ϑ_l, + Y.soil.θ_i, + soil_params.ρc_ds, + soil_params.earth_param_set, + ) + Y.soil.ρe_int .= + Soil.volumetric_internal_energy.( + Y.soil.θ_i, + ρc_s, + T, + soil_params.earth_param_set, + ) + Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air + Y.canopy.hydraulics.ϑ_l.:1 .= plant_ν + evaluate!(Y.canopy.energy.T, atmos.T, t0) + + set_initial_cache! = make_set_initial_cache(land) + exp_tendency! = make_exp_tendency(land) + imp_tendency! = ClimaLand.make_imp_tendency(land) + jacobian! = ClimaLand.make_jacobian(land) + set_initial_cache!(p, Y, t0) + + # set up jacobian info + jac_kwargs = + (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = jacobian!) + + prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction( + T_exp! = exp_tendency!, + T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), + dss! = ClimaLand.dss!, + ), + Y, + (t0, tf), + p, + ) + + updateat = Array(t0:(3600 * 3):tf) + drivers = ClimaLand.get_drivers(land) + updatefunc = ClimaLand.make_update_drivers(drivers) + + # ClimaDiagnostics + + nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(subsurface_space, outdir) + + diags = ClimaLand.default_diagnostics( + land, + t0, + ref_time; + output_writer = nc_writer, + output_vars = :short, + average_period = :monthly, + ) + + diagnostic_handler = + ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt) + + diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler) + + driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) + return prob, SciMLBase.CallbackSet(driver_cb, diag_cb) +end + +function setup_and_solve_problem(; greet = false) + + t0 = 0.0 + tf = 60 * 60.0 * 24 * 334 + Δt = 900.0 + nelements = (10, 10, 15) + if greet + @info "Run: Regional Soil-Canopy Model" + @info "Resolution: $nelements" + @info "Timestep: $Δt s" + @info "Duration: $(tf - t0) s" + end + + prob, cb = setup_prob(t0, tf, Δt; nelements) + + # Define timestepper and ODE algorithm + stepper = CTS.ARS111() + ode_algo = CTS.IMEXAlgorithm( + stepper, + CTS.NewtonsMethod( + max_iters = 3, + update_j = CTS.UpdateEvery(CTS.NewNewtonIteration), + ), + ) + SciMLBase.solve(prob, ode_algo; dt = Δt, callback = cb, adaptive = false) + return nothing +end + +setup_and_solve_problem(; greet = true); +# read in diagnostics and make some plots! +#### ClimaAnalysis #### +simdir = ClimaAnalysis.SimDir(outdir) +short_names = ["gpp", "ct", "swc", "si"] +for short_name in short_names + var = get(simdir; short_name) + times = ClimaAnalysis.times(var) + for t in times + fig = CairoMakie.Figure(size = (800, 600)) + kwargs = ClimaAnalysis.has_altitude(var) ? Dict(:z => 1) : Dict() + tmp = ClimaAnalysis.slice(var, time = t; kwargs...) + if ~all(isnan.(tmp.data)) + viz.heatmap2D!(fig, tmp) + CairoMakie.save(joinpath(root_path, "$(short_name)_$t.png"), fig) + end + end +end diff --git a/experiments/standalone/Bucket/bucket_era5.jl b/experiments/standalone/Bucket/bucket_era5.jl index d76ab40de6..0291f30a29 100644 --- a/experiments/standalone/Bucket/bucket_era5.jl +++ b/experiments/standalone/Bucket/bucket_era5.jl @@ -298,11 +298,11 @@ if regional_simulation ClimaCore.Geometry._coordinate( pt::ClimaCore.Geometry.LatLongPoint, ::Val{1}, - ) = ClimaCore.Geometry.LongPoint(pt.long) + ) = ClimaCore.Geometry.LatPoint(pt.lat) ClimaCore.Geometry._coordinate( pt::ClimaCore.Geometry.LatLongPoint, ::Val{2}, - ) = ClimaCore.Geometry.LatPoint(pt.lat) + ) = ClimaCore.Geometry.LongPoint(pt.long) else longpts = range(-180.0, 180.0, num_pts) diff --git a/src/shared_utilities/Domains.jl b/src/shared_utilities/Domains.jl index 76e5d5654f..9cc109afeb 100644 --- a/src/shared_utilities/Domains.jl +++ b/src/shared_utilities/Domains.jl @@ -267,22 +267,29 @@ function Plane(; @assert periodic == (false, false) radius_earth = FT(radius_earth) long, lat = longlat + dxlim = xlim # long bounds + dylim = ylim # lat bounds + # Now make x refer to lat, and y refer to long, + # for compatibility with ClimaCore xlim = - (long - xlim[1] / (2radius_earth), long + xlim[2] / (2radius_earth)) - ylim = - (lat - ylim[1] / (2radius_earth), lat + ylim[2] / (2radius_earth)) + (lat - dylim[1] / (2radius_earth), lat + dylim[2] / (2radius_earth)) + ylim = ( + long - dxlim[1] / (2radius_earth), + long + dxlim[2] / (2radius_earth), + ) @assert xlim[1] < xlim[2] @assert ylim[1] < ylim[2] + # NOTE: We have LatLong instead of the other way because of ClimaCore domain_x = ClimaCore.Domains.IntervalDomain( - ClimaCore.Geometry.LongPoint(xlim[1]), - ClimaCore.Geometry.LongPoint(xlim[2]); - boundary_names = (:west, :east), + ClimaCore.Geometry.LatPoint(xlim[1]), + ClimaCore.Geometry.LatPoint(xlim[2]); + boundary_names = (:north, :south), ) domain_y = ClimaCore.Domains.IntervalDomain( - ClimaCore.Geometry.LatPoint(ylim[1]), - ClimaCore.Geometry.LatPoint(ylim[2]); - boundary_names = (:north, :south), + ClimaCore.Geometry.LongPoint(ylim[1]), + ClimaCore.Geometry.LongPoint(ylim[2]); + boundary_names = (:west, :east), ) end plane = ClimaCore.Domains.RectangleDomain(domain_x, domain_y) @@ -922,10 +929,10 @@ end """ depth(space::Union{ClimaCore.Spaces.CenterFiniteDifferenceSpace, - ClimaCore.Spaces.CenterExtrudedFiniteDifferenceSpace}) + ClimaCore.Spaces.CenterExtrudedFiniteDifferenceSpace}) Returns the depth of the domain as a scalar. Note that these functions -will need to be modified upon the introduction of +will need to be modified upon the introduction of - topography at surface - depth to bedrock (topography at bottom) diff --git a/test/shared_utilities/domains.jl b/test/shared_utilities/domains.jl index 66db30e853..4b1b03d51f 100644 --- a/test/shared_utilities/domains.jl +++ b/test/shared_utilities/domains.jl @@ -172,17 +172,17 @@ FT = Float32 ClimaCore.Spaces.SpectralElementSpace2D # Plane latlong - dxlim = (FT(-50_000), FT(80_000)) - dylim = (FT(-30_000), FT(40_000)) + dxlim = (FT(50_000), FT(80_000)) + dylim = (FT(30_000), FT(40_000)) longlat = (FT(-118.14452), FT(34.14778)) radius_earth = FT(6.378e6) xlim_longlat = ( - longlat[1] - dxlim[1] / 2radius_earth, - longlat[1] + dxlim[2] / 2radius_earth, + longlat[2] - dylim[1] / 2radius_earth, + longlat[2] + dylim[2] / 2radius_earth, ) ylim_longlat = ( - longlat[2] - dylim[1] / 2radius_earth, - longlat[2] + dylim[2] / (2radius_earth), + longlat[1] - dxlim[1] / 2radius_earth, + longlat[1] + dxlim[2] / 2radius_earth, ) longlat_plane = Plane(;