From eadfba174891e684be2935c95dc1bc8cf28a75df Mon Sep 17 00:00:00 2001 From: kmdeck Date: Fri, 21 Jun 2024 10:33:11 -0700 Subject: [PATCH] jacobian wip --- .buildkite/Manifest.toml | 71 ++++++---- docs/Manifest.toml | 129 ++++++++++-------- .../fluxnet/US-Ha1/US-Ha1_simulation.jl | 4 +- .../fluxnet/US-MOz/US-MOz_simulation.jl | 4 +- .../fluxnet/US-NR1/US-NR1_simulation.jl | 4 +- .../fluxnet/US-Var/US-Var_simulation.jl | 4 +- .../integrated/fluxnet/fluxnet_simulation.jl | 4 +- .../performance/profile_allocations.jl | 4 +- src/shared_utilities/implicit_timestepping.jl | 69 ++++++---- src/standalone/Vegetation/Canopy.jl | 67 +++++++++ .../Vegetation/canopy_boundary_fluxes.jl | 16 ++- src/standalone/Vegetation/canopy_energy.jl | 92 +++++++++++-- src/standalone/Vegetation/component_models.jl | 41 +++++- test/standalone/Vegetation/canopy_model.jl | 8 ++ 14 files changed, 381 insertions(+), 136 deletions(-) diff --git a/.buildkite/Manifest.toml b/.buildkite/Manifest.toml index 0ffdeb970a..e7e54df2b4 100644 --- a/.buildkite/Manifest.toml +++ b/.buildkite/Manifest.toml @@ -227,9 +227,13 @@ version = "0.1.6" [[deps.BlockArrays]] deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"] -git-tree-sha1 = "9a9610fbe5779636f75229e423e367124034af41" +git-tree-sha1 = "5c0ffe1dff8cb7112de075f1b1cb32191675fcba" uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" -version = "0.16.43" +version = "1.1.0" +weakdeps = ["BandedMatrices"] + + [deps.BlockArrays.extensions] + BlockArraysBandedMatricesExt = "BandedMatrices" [[deps.Bzip2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -376,10 +380,12 @@ weakdeps = ["CUDA", "MPI"] ClimaCommsMPIExt = "MPI" [[deps.ClimaCore]] -deps = ["Adapt", "BandedMatrices", "BlockArrays", "ClimaComms", "CubedSphere", "DataStructures", "DocStringExtensions", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "KrylovKit", "LinearAlgebra", "MultiBroadcastFusion", "NVTX", "PkgVersion", "RecursiveArrayTools", "RootSolvers", "SparseArrays", "Static", "StaticArrays", "Statistics", "Unrolled"] -git-tree-sha1 = "ffd0b5afde1816c9e61697e9510fbd9aa01b62de" +deps = ["Adapt", "BandedMatrices", "BlockArrays", "ClimaComms", "CubedSphere", "DataStructures", "DocStringExtensions", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "KrylovKit", "LinearAlgebra", "MultiBroadcastFusion", "NVTX", "PkgVersion", "RecursiveArrayTools", "RootSolvers", "SparseArrays", "StaticArrays", "Statistics", "Unrolled"] +git-tree-sha1 = "d16d690d255c05cbdd448b2c694c2500141eef3b" +repo-rev = "kd/matrix_fields" +repo-url = "https://github.com/CliMA/ClimaCore.jl.git" uuid = "d414da3d-4745-48bb-8d80-42e94e092884" -version = "0.14.9" +version = "0.14.10" weakdeps = ["CUDA", "Krylov"] [deps.ClimaCore.extensions] @@ -518,6 +524,11 @@ git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" version = "0.3.0" +[[deps.CommonWorldInvalidations]] +git-tree-sha1 = "ae52d1c52048455e85a387fbee9be553ec2b68d0" +uuid = "f70d9fcc-98c5-4d4a-abd7-e4cdeebd8ca8" +version = "1.0.0" + [[deps.Compat]] deps = ["TOML", "UUIDs"] git-tree-sha1 = "b1c55339b7c6c350ee89f2c1604299660525b248" @@ -869,9 +880,9 @@ version = "3.3.10+0" [[deps.FLoops]] deps = ["BangBang", "Compat", "FLoopsBase", "InitialValues", "JuliaVariables", "MLStyle", "Serialization", "Setfield", "Transducers"] -git-tree-sha1 = "ffb97765602e3cbe59a0589d237bf07f245a8576" +git-tree-sha1 = "0a2e5873e9a5f54abb06418d57a8df689336a660" uuid = "cc61a311-1640-44b5-9fba-1b764f453329" -version = "0.2.1" +version = "0.2.2" [[deps.FLoopsBase]] deps = ["ContextVariablesX"] @@ -1421,10 +1432,14 @@ uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" version = "0.9.6" [[deps.KrylovKit]] -deps = ["ChainRulesCore", "GPUArraysCore", "LinearAlgebra", "Printf", "VectorInterface"] -git-tree-sha1 = "3f3a92bbe8f568b689a7f7bc193f7c717d793751" +deps = ["GPUArraysCore", "LinearAlgebra", "PackageExtensionCompat", "Printf", "VectorInterface"] +git-tree-sha1 = "3c2a016489c38f35160a246c91a3f3353c47bb68" uuid = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" -version = "0.7.1" +version = "0.8.1" +weakdeps = ["ChainRulesCore"] + + [deps.KrylovKit.extensions] + KrylovKitChainRulesCoreExt = "ChainRulesCore" [[deps.LAME_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -1506,15 +1521,23 @@ uuid = "10f19ff3-798f-405d-979b-55457f8fc047" version = "0.1.17" [[deps.LazyArrays]] -deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "MatrixFactorizations", "SparseArrays"] -git-tree-sha1 = "35079a6a869eecace778bcda8641f9a54ca3a828" +deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "SparseArrays"] +git-tree-sha1 = "46dd13736e33cc3bfc610f62b6c7f84b9c95539a" uuid = "5078a376-72f3-5289-bfd5-ec5146d43c02" -version = "1.10.0" -weakdeps = ["StaticArrays"] +version = "2.1.8" [deps.LazyArrays.extensions] + LazyArraysBandedMatricesExt = "BandedMatrices" + LazyArraysBlockArraysExt = "BlockArrays" + LazyArraysBlockBandedMatricesExt = "BlockBandedMatrices" LazyArraysStaticArraysExt = "StaticArrays" + [deps.LazyArrays.weakdeps] + BandedMatrices = "aae01518-5342-5314-be14-df237901396f" + BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" + BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + [[deps.LazyArtifacts]] deps = ["Artifacts", "Pkg"] uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" @@ -1627,9 +1650,9 @@ version = "2.8.0" [[deps.LinearSolve]] 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 = "7648cc20100504f4b453917aacc8520e9c0ecfb3" +git-tree-sha1 = "b2e2dba60642e07c062eb3143770d7e234316772" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -version = "2.30.1" +version = "2.30.2" [deps.LinearSolve.extensions] LinearSolveBandedMatricesExt = "BandedMatrices" @@ -1785,12 +1808,6 @@ git-tree-sha1 = "96ca8a313eb6437db5ffe946c457a401bbb8ce1d" uuid = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53" version = "0.5.7" -[[deps.MatrixFactorizations]] -deps = ["ArrayLayouts", "LinearAlgebra", "Printf", "Random"] -git-tree-sha1 = "6731e0574fa5ee21c02733e397beb133df90de35" -uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87" -version = "2.2.0" - [[deps.MaybeInplace]] deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "SparseArrays"] git-tree-sha1 = "1b9e613f2ca3b6cdcbfe36381e17ca2b66d4b3a1" @@ -1878,9 +1895,9 @@ version = "4.5.1" [[deps.NNlib]] deps = ["Adapt", "Atomix", "ChainRulesCore", "GPUArraysCore", "KernelAbstractions", "LinearAlgebra", "Pkg", "Random", "Requires", "Statistics"] -git-tree-sha1 = "1288e6db94d98f7d194454452176b82edb25b32c" +git-tree-sha1 = "333cd68c28bb57ac23c50acd1036ecd7dd78bf57" uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" -version = "0.9.19" +version = "0.9.20" [deps.NNlib.extensions] NNlibAMDGPUExt = "AMDGPU" @@ -2632,10 +2649,10 @@ uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15" version = "0.1.1" [[deps.Static]] -deps = ["IfElse"] -git-tree-sha1 = "d2fdac9ff3906e27f7a618d47b676941baa6c80c" +deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools"] +git-tree-sha1 = "0bbff21027dd8a107551847528127b62a35f7594" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.8.10" +version = "1.1.0" [[deps.StaticArrayInterface]] deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"] diff --git a/docs/Manifest.toml b/docs/Manifest.toml index f93e82c16c..c3defdee7d 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -123,9 +123,9 @@ version = "7.12.0" [[deps.ArrayLayouts]] deps = ["FillArrays", "LinearAlgebra"] -git-tree-sha1 = "8556500c18fcad8b4c44058e23fbc4a36143f6be" +git-tree-sha1 = "ce2ca959f932f5dad70697dd93133d1167cf1e4e" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" -version = "1.10.1" +version = "1.10.2" weakdeps = ["SparseArrays"] [deps.ArrayLayouts.extensions] @@ -221,9 +221,13 @@ version = "0.1.6" [[deps.BlockArrays]] deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"] -git-tree-sha1 = "9a9610fbe5779636f75229e423e367124034af41" +git-tree-sha1 = "5c0ffe1dff8cb7112de075f1b1cb32191675fcba" uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" -version = "0.16.43" +version = "1.1.0" +weakdeps = ["BandedMatrices"] + + [deps.BlockArrays.extensions] + BlockArraysBandedMatricesExt = "BandedMatrices" [[deps.Blosc_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Zlib_jll", "Zstd_jll"] @@ -277,9 +281,9 @@ version = "0.10.14" [[deps.CUDA]] deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "StaticArrays", "Statistics"] -git-tree-sha1 = "6e945e876652f2003e6ca74e19a3c45017d3e9f6" +git-tree-sha1 = "fdd9dfb67dfefd548f51000cc400bb51003de247" uuid = "052768ef-5323-5732-b1bb-66c8b64840ba" -version = "5.4.2" +version = "5.4.3" weakdeps = ["ChainRulesCore", "EnzymeCore", "SpecialFunctions"] [deps.CUDA.extensions] @@ -319,9 +323,9 @@ version = "1.0.5" [[deps.CairoMakie]] deps = ["CRC32c", "Cairo", "Cairo_jll", "Colors", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "PrecompileTools"] -git-tree-sha1 = "f84837ccd1411ba059bb0b752dab9c7f1b0b0826" +git-tree-sha1 = "e4da5095557f24713bae4c9f50e34ff4d3b959c0" uuid = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" -version = "0.12.4" +version = "0.12.5" [[deps.Cairo_jll]] deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] @@ -365,10 +369,12 @@ version = "0.6.3" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" [[deps.ClimaCore]] -deps = ["Adapt", "BandedMatrices", "BlockArrays", "ClimaComms", "CubedSphere", "DataStructures", "DocStringExtensions", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "KrylovKit", "LinearAlgebra", "MultiBroadcastFusion", "NVTX", "PkgVersion", "RecursiveArrayTools", "RootSolvers", "SparseArrays", "Static", "StaticArrays", "Statistics", "Unrolled"] -git-tree-sha1 = "ffd0b5afde1816c9e61697e9510fbd9aa01b62de" +deps = ["Adapt", "BandedMatrices", "BlockArrays", "ClimaComms", "CubedSphere", "DataStructures", "DocStringExtensions", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "KrylovKit", "LinearAlgebra", "MultiBroadcastFusion", "NVTX", "PkgVersion", "RecursiveArrayTools", "RootSolvers", "SparseArrays", "StaticArrays", "Statistics", "Unrolled"] +git-tree-sha1 = "d16d690d255c05cbdd448b2c694c2500141eef3b" +repo-rev = "kd/matrix_fields" +repo-url = "https://github.com/CliMA/ClimaCore.jl.git" uuid = "d414da3d-4745-48bb-8d80-42e94e092884" -version = "0.14.9" +version = "0.14.10" weakdeps = ["CUDA", "Krylov"] [deps.ClimaCore.extensions] @@ -501,6 +507,11 @@ git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" version = "0.3.0" +[[deps.CommonWorldInvalidations]] +git-tree-sha1 = "ae52d1c52048455e85a387fbee9be553ec2b68d0" +uuid = "f70d9fcc-98c5-4d4a-abd7-e4cdeebd8ca8" +version = "1.0.0" + [[deps.Compat]] deps = ["TOML", "UUIDs"] git-tree-sha1 = "b1c55339b7c6c350ee89f2c1604299660525b248" @@ -608,9 +619,9 @@ version = "0.1.2" [[deps.DelaunayTriangulation]] deps = ["EnumX", "ExactPredicates", "Random"] -git-tree-sha1 = "b0cb128d2e100646573e1da8565b02491fddb5ef" +git-tree-sha1 = "078c716cbb032242df18b960e8b1fec6b1b0b9f9" uuid = "927a84f5-c5f4-47a5-9785-b46e178433df" -version = "1.0.4" +version = "1.0.5" [[deps.DelimitedFiles]] deps = ["Mmap"] @@ -847,9 +858,9 @@ version = "3.3.10+0" [[deps.FLoops]] deps = ["BangBang", "Compat", "FLoopsBase", "InitialValues", "JuliaVariables", "MLStyle", "Serialization", "Setfield", "Transducers"] -git-tree-sha1 = "ffb97765602e3cbe59a0589d237bf07f245a8576" +git-tree-sha1 = "0a2e5873e9a5f54abb06418d57a8df689336a660" uuid = "cc61a311-1640-44b5-9fba-1b764f453329" -version = "0.2.1" +version = "0.2.2" [[deps.FLoopsBase]] deps = ["ContextVariablesX"] @@ -1049,10 +1060,10 @@ uuid = "46192b85-c4d5-4398-a991-12ede77f4527" version = "0.1.6" [[deps.GPUCompiler]] -deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "TimerOutputs", "UUIDs"] -git-tree-sha1 = "518ebd058c9895de468a8c255797b0c53fdb44dd" +deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Preferences", "Scratch", "Serialization", "TOML", "TimerOutputs", "UUIDs"] +git-tree-sha1 = "ab29216184312f99ff957b32cd63c2fe9c928b91" uuid = "61eb1bfa-7361-4325-ad38-22787b887f55" -version = "0.26.5" +version = "0.26.7" [[deps.GR]] deps = ["Artifacts", "Base64", "DelimitedFiles", "Downloads", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Preferences", "Printf", "Random", "Serialization", "Sockets", "TOML", "Tar", "Test", "p7zip_jll"] @@ -1288,9 +1299,9 @@ weakdeps = ["ClimaParams"] [[deps.IntelOpenMP_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "be50fe8df3acbffa0274a744f1a99d29c45a57f4" +git-tree-sha1 = "14eb2b542e748570b56446f4c50fbfb2306ebc45" uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" -version = "2024.1.0+0" +version = "2024.2.0+0" [[deps.InteractiveUtils]] deps = ["Markdown"] @@ -1430,10 +1441,14 @@ uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" version = "0.9.6" [[deps.KrylovKit]] -deps = ["ChainRulesCore", "GPUArraysCore", "LinearAlgebra", "Printf", "VectorInterface"] -git-tree-sha1 = "3f3a92bbe8f568b689a7f7bc193f7c717d793751" +deps = ["GPUArraysCore", "LinearAlgebra", "PackageExtensionCompat", "Printf", "VectorInterface"] +git-tree-sha1 = "3c2a016489c38f35160a246c91a3f3353c47bb68" uuid = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" -version = "0.7.1" +version = "0.8.1" +weakdeps = ["ChainRulesCore"] + + [deps.KrylovKit.extensions] + KrylovKitChainRulesCoreExt = "ChainRulesCore" [[deps.LAME_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -1449,9 +1464,9 @@ version = "3.0.0+1" [[deps.LLVM]] deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"] -git-tree-sha1 = "389aea28d882a40b5e1747069af71bdbd47a1cae" +git-tree-sha1 = "020abd49586480c1be84f57da0017b5d3db73f7c" uuid = "929cbde3-209d-540e-8aea-75f648917ca0" -version = "7.2.1" +version = "8.0.0" weakdeps = ["BFloat16s"] [deps.LLVM.extensions] @@ -1459,9 +1474,9 @@ weakdeps = ["BFloat16s"] [[deps.LLVMExtra_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] -git-tree-sha1 = "88b916503aac4fb7f701bb625cd84ca5dd1677bc" +git-tree-sha1 = "c2636c264861edc6d305e6b4d528f09566d24c5e" uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab" -version = "0.0.29+0" +version = "0.0.30+0" [[deps.LLVMLoopInfo]] git-tree-sha1 = "2e5c102cfc41f48ae4740c7eca7743cc7e7b75ea" @@ -1520,15 +1535,23 @@ uuid = "0e77f7df-68c5-4e49-93ce-4cd80f5598bf" version = "1.2.2" [[deps.LazyArrays]] -deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "MatrixFactorizations", "SparseArrays"] -git-tree-sha1 = "35079a6a869eecace778bcda8641f9a54ca3a828" +deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "SparseArrays"] +git-tree-sha1 = "46dd13736e33cc3bfc610f62b6c7f84b9c95539a" uuid = "5078a376-72f3-5289-bfd5-ec5146d43c02" -version = "1.10.0" -weakdeps = ["StaticArrays"] +version = "2.1.8" [deps.LazyArrays.extensions] + LazyArraysBandedMatricesExt = "BandedMatrices" + LazyArraysBlockArraysExt = "BlockArrays" + LazyArraysBlockBandedMatricesExt = "BlockBandedMatrices" LazyArraysStaticArraysExt = "StaticArrays" + [deps.LazyArrays.weakdeps] + BandedMatrices = "aae01518-5342-5314-be14-df237901396f" + BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" + BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + [[deps.LazyArtifacts]] deps = ["Artifacts", "Pkg"] uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" @@ -1641,9 +1664,9 @@ version = "2.8.0" [[deps.LinearSolve]] 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 = "7648cc20100504f4b453917aacc8520e9c0ecfb3" +git-tree-sha1 = "b2e2dba60642e07c062eb3143770d7e234316772" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -version = "2.30.1" +version = "2.30.2" [deps.LinearSolve.extensions] LinearSolveBandedMatricesExt = "BandedMatrices" @@ -1678,9 +1701,9 @@ version = "2.30.1" [[deps.Literate]] deps = ["Base64", "IOCapture", "JSON", "REPL"] -git-tree-sha1 = "596df2daea9c27da81eee63ef2cf101baf10c24c" +git-tree-sha1 = "eef2e1fc1dc38af90a18eb16e519e06d1fd10c2a" uuid = "98b081ad-f1c9-55d3-8b20-4c87d4299306" -version = "2.18.0" +version = "2.19.0" [[deps.LogExpFunctions]] deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] @@ -1726,9 +1749,9 @@ version = "1.9.4+0" [[deps.MKL_jll]] deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"] -git-tree-sha1 = "80b2833b56d466b3858d565adcd16a4a05f2089b" +git-tree-sha1 = "f046ccd0c6db2832a9f639e2c669c6fe867e5f4f" uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" -version = "2024.1.0+0" +version = "2024.2.0+0" [[deps.MLStyle]] git-tree-sha1 = "bc38dff0548128765760c79eb7388a4b37fae2c8" @@ -1767,15 +1790,15 @@ version = "0.5.13" [[deps.Makie]] deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "Dates", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun", "Unitful"] -git-tree-sha1 = "57a1a2b3d12e04f9e9fb77d61cd12571d5541c5f" +git-tree-sha1 = "863b9e666b5a099c8835e85476a5834f9d77c4c1" uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" -version = "0.21.4" +version = "0.21.5" [[deps.MakieCore]] deps = ["ColorTypes", "GeometryBasics", "IntervalSets", "Observables"] -git-tree-sha1 = "638bc817096742e8302f7b0b972ee5701fe00e97" +git-tree-sha1 = "c1c950560397ee68ad7302ee0e3efa1b07466a2f" uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b" -version = "0.8.3" +version = "0.8.4" [[deps.ManualMemory]] git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd" @@ -1803,12 +1826,6 @@ git-tree-sha1 = "1865d0b8a2d91477c8b16b49152a32764c7b1f5f" uuid = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53" version = "0.6.0" -[[deps.MatrixFactorizations]] -deps = ["ArrayLayouts", "LinearAlgebra", "Printf", "Random"] -git-tree-sha1 = "6731e0574fa5ee21c02733e397beb133df90de35" -uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87" -version = "2.2.0" - [[deps.MaybeInplace]] deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "SparseArrays"] git-tree-sha1 = "1b9e613f2ca3b6cdcbfe36381e17ca2b66d4b3a1" @@ -1890,9 +1907,9 @@ version = "7.8.3" [[deps.NNlib]] deps = ["Adapt", "Atomix", "ChainRulesCore", "GPUArraysCore", "KernelAbstractions", "LinearAlgebra", "Pkg", "Random", "Requires", "Statistics"] -git-tree-sha1 = "1288e6db94d98f7d194454452176b82edb25b32c" +git-tree-sha1 = "333cd68c28bb57ac23c50acd1036ecd7dd78bf57" uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" -version = "0.9.19" +version = "0.9.20" [deps.NNlib.extensions] NNlibAMDGPUExt = "AMDGPU" @@ -2346,9 +2363,9 @@ version = "0.6.12" [[deps.RecursiveArrayTools]] deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "SparseArrays", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"] -git-tree-sha1 = "76a35763102e16c97d7302f5dc67837d24802035" +git-tree-sha1 = "b450d967a770fb13d0e26358f58375e20361cf9c" uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" -version = "3.25.0" +version = "3.26.0" [deps.RecursiveArrayTools.extensions] RecursiveArrayToolsFastBroadcastExt = "FastBroadcast" @@ -2492,9 +2509,9 @@ version = "1.2.1" [[deps.SentinelArrays]] deps = ["Dates", "Random"] -git-tree-sha1 = "6bb314cb1aacfa37ef58e5a0ccf4a1ec0311f495" +git-tree-sha1 = "ff11acffdb082493657550959d4feb4b6149e73a" uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c" -version = "1.4.4" +version = "1.4.5" [[deps.Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -2647,10 +2664,10 @@ uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15" version = "0.1.1" [[deps.Static]] -deps = ["IfElse"] -git-tree-sha1 = "d2fdac9ff3906e27f7a618d47b676941baa6c80c" +deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools"] +git-tree-sha1 = "0bbff21027dd8a107551847528127b62a35f7594" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.8.10" +version = "1.1.0" [[deps.StaticArrayInterface]] deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"] diff --git a/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl b/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl index 6afcd4670d..db24273f77 100644 --- a/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl +++ b/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl @@ -20,7 +20,7 @@ h_stem = FT(14) # m t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow # Time step size: -dt = Float64(40) +dt = Float64(450) # Number of timesteps between saving output: -n = 45 +n = 4 diff --git a/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl b/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl index b229ab5024..810e6bd9ae 100644 --- a/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl +++ b/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl @@ -20,7 +20,7 @@ h_leaf = FT(9.5) # m t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow # Time step size: -dt = Float64(450) +dt = Float64(900) # Number of timesteps between saving output: -n = 4 +n = 2 diff --git a/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl b/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl index e27c0b4827..4644d7de9f 100644 --- a/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl +++ b/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl @@ -20,7 +20,7 @@ h_stem = FT(7.5) # m t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow # Time step size: -dt = Float64(40) +dt = Float64(450) # Number of timesteps between saving output: -n = 45 +n = 4 diff --git a/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl b/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl index 1ba85752a1..b05140d0d0 100644 --- a/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl +++ b/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl @@ -20,7 +20,7 @@ h_stem = FT(0) # m t0 = Float64(21 * 3600 * 24)# start day 21 of the year # Time step size: -dt = Float64(20) +dt = Float64(900) # Number of timesteps between saving output: -n = 45 +n = 1 diff --git a/experiments/integrated/fluxnet/fluxnet_simulation.jl b/experiments/integrated/fluxnet/fluxnet_simulation.jl index 818e9af174..e0554c63ae 100644 --- a/experiments/integrated/fluxnet/fluxnet_simulation.jl +++ b/experiments/integrated/fluxnet/fluxnet_simulation.jl @@ -13,7 +13,7 @@ timestepper = CTS.ARS343(); ode_algo = CTS.IMEXAlgorithm( timestepper, CTS.NewtonsMethod( - max_iters = 1, - update_j = CTS.UpdateEvery(CTS.NewNewtonIteration), + max_iters = 2, + update_j = CTS.UpdateEvery(CTS.NewTimeStep), ), ); diff --git a/experiments/integrated/performance/profile_allocations.jl b/experiments/integrated/performance/profile_allocations.jl index e9b898383c..b8ef791e67 100644 --- a/experiments/integrated/performance/profile_allocations.jl +++ b/experiments/integrated/performance/profile_allocations.jl @@ -326,7 +326,7 @@ imp_tendency! = make_imp_tendency(land); jacobian! = make_jacobian(land); # Set up timestepping and simulation callbacks -dt = Float64(150) +dt = Float64(900) tf = Float64(t0 + 100dt) saveat = Array(t0:dt:tf) updateat = Array(t0:dt:tf) @@ -334,7 +334,7 @@ timestepper = CTS.ARS343() ode_algo = CTS.IMEXAlgorithm( timestepper, CTS.NewtonsMethod( - max_iters = 1, + max_iters = 2, update_j = CTS.UpdateEvery(CTS.NewNewtonIteration), ), ); diff --git a/src/shared_utilities/implicit_timestepping.jl b/src/shared_utilities/implicit_timestepping.jl index c4343c038b..bf0a44a15a 100644 --- a/src/shared_utilities/implicit_timestepping.jl +++ b/src/shared_utilities/implicit_timestepping.jl @@ -1,5 +1,6 @@ using ClimaCore.MatrixFields import ClimaCore.MatrixFields: @name, ⋅ +using ClimaCore: Spaces import UnrolledUtilities import LinearAlgebra import LinearAlgebra: I @@ -108,51 +109,63 @@ to the `explicit_names` tuple. """ function ImplicitEquationJacobian(Y::ClimaCore.Fields.FieldVector) FT = eltype(Y) - center_space = axes(Y.soil.ϑ_l) - - # Construct a tridiagonal matrix that will be used as the Jacobian - tridiag_type = MatrixFields.TridiagonalMatrixRow{FT} - # Create a field containing a `TridiagonalMatrixRow` at each point - tridiag_field = Fields.Field(tridiag_type, center_space) - fill!(parent(tridiag_field), NaN) - # Only add jacobian blocks for fields that are in Y for this model - is_in_Y(name) = MatrixFields.has_field(Y, name) + is_in_Y(var) = MatrixFields.has_field(Y, var) # Define the implicit and explicit variables of any model we use - implicit_names = (@name(soil.ϑ_l), @name(soil.ρe_int)) + implicit_names = (@name(soil.ϑ_l), @name(soil.ρe_int), @name(canopy.energy.T)) explicit_names = ( @name(soilco2.C), @name(soil.θ_i), @name(canopy.hydraulics.ϑ_l), - @name(canopy.energy.T), ) # Filter out the variables that are not in this model's state, `Y` - available_implicit_names = - MatrixFields.unrolled_filter(is_in_Y, implicit_names) - available_explicit_names = - MatrixFields.unrolled_filter(is_in_Y, explicit_names) + available_implicit_vars = + MatrixFields.unrolled_filter(is_in_Y, implicit_vars) + available_explicit_vars = + MatrixFields.unrolled_filter(is_in_Y, explicit_vars) + + + function get_j_field( + space::Union{ + Spaces.FiniteDifferenceSpace, + Spaces.ExtrudedFiniteDifferenceSpace, + }, + ) + # Construct a tridiagonal matrix that will be used as the Jacobian + tridiag_type = MatrixFields.TridiagonalMatrixRow{FT} + # Create a field containing a `TridiagonalMatrixRow` at each point + tridiag_field = ClimaCore.Fields.Field(tridiag_type, space) + fill!(parent(tridiag_field), 0) + return tridiag_field + end + function get_j_field( + space::Union{Spaces.PointSpace, Spaces.SpectralElementSpace2D}, + ) + diag_type = MatrixFields.DiagonalMatrixRow{FT} + diag_field = ClimaCore.Fields.Field(diag_type, space) + fill!(parent(diag_field), 0) + return diag_field + end + + implicit_blocks = MatrixFields.unrolled_map( + var -> + (var, var) => get_j_field(axes(MatrixFields.get_field(Y, var))), + available_implicit_vars, + ) # For explicitly-stepped variables, use the negative identity matrix # Note: We have to use FT(-1) * I instead of -I because inv(-1) == -1.0, # which means that multiplying inv(-1) by a Float32 will yield a Float64. - identity_blocks = MatrixFields.unrolled_map( - name -> (name, name) => FT(-1) * I, - available_explicit_names, + explicit_blocks = MatrixFields.unrolled_map( + var -> (var, var) => FT(-1) * I, + available_explicit_vars, ) - # For implicitly-stepped variables, use a tridiagonal matrix - tridiagonal_blocks = MatrixFields.unrolled_map( - name -> - (name, name) => Fields.Field( - tridiag_type, - axes(MatrixFields.get_field(Y, name)), - ), - available_implicit_names, - ) - matrix = MatrixFields.FieldMatrix(identity_blocks..., tridiagonal_blocks...) + + matrix = MatrixFields.FieldMatrix(implicit_blocks..., explicit_blocks...) # Set up block diagonal solver for block Jacobian alg = MatrixFields.BlockDiagonalSolve() diff --git a/src/standalone/Vegetation/Canopy.jl b/src/standalone/Vegetation/Canopy.jl index b9d5c3a5db..8bdc4b4408 100644 --- a/src/standalone/Vegetation/Canopy.jl +++ b/src/standalone/Vegetation/Canopy.jl @@ -3,6 +3,9 @@ using DocStringExtensions using Thermodynamics using ClimaLand using ClimaCore +using ClimaCore.MatrixFields +import ClimaCore.MatrixFields: @name, ⋅ +import LinearAlgebra: I using ClimaLand: AbstractRadiativeDrivers, AbstractAtmosphericDrivers import ..Parameters as LP @@ -20,6 +23,8 @@ import ClimaLand: make_update_boundary_fluxes, make_update_aux, make_compute_exp_tendency, + make_compute_imp_tendency, + make_compute_jacobian, get_drivers using ClimaLand.Domains: Point, Plane, SphericalSurface @@ -605,6 +610,68 @@ function make_compute_exp_tendency( end return compute_exp_tendency! end + +""" + make_compute_imp_tendency(canopy::CanopyModel) + +Creates and returns the compute_imp_tendency! for the `CanopyModel`. +""" +function make_compute_imp_tendency( + canopy::CanopyModel{ + FT, + <:AutotrophicRespirationModel, + <:Union{BeerLambertModel, TwoStreamModel}, + <:Union{FarquharModel, OptimalityFarquharModel}, + <:MedlynConductanceModel, + <:PlantHydraulicsModel, + <:Union{PrescribedCanopyTempModel, BigLeafEnergyModel}, + }, +) where {FT} + components = canopy_components(canopy) + compute_imp_tendency_list = map( + x -> make_compute_imp_tendency(getproperty(canopy, x), canopy), + components, + ) + function compute_imp_tendency!(dY, Y, p, t) + for f! in compute_imp_tendency_list + f!(dY, Y, p, t) + end + + end + return compute_imp_tendency! +end + +""" + make_update_jacobian(canopy::CanopyModel) + +Creates and returns the update_jacobian! for the `CanopyModel`. +""" +function make_update_jacobian( + canopy::CanopyModel{ + FT, + <:AutotrophicRespirationModel, + <:Union{BeerLambertModel, TwoStreamModel}, + <:Union{FarquharModel, OptimalityFarquharModel}, + <:MedlynConductanceModel, + <:PlantHydraulicsModel, + <:Union{PrescribedCanopyTempModel, BigLeafEnergyModel}, + }, +) where {FT} + components = canopy_components(canopy) + update_jacobian_list = map( + x -> make_update_jacobian(getproperty(canopy, x), canopy), + components, + ) + function update_jacobian!(W, Y, p, dtγ, t) + for f! in update_jacobian_list + f!(W, Y, p, dtγ, t) + end + + end + return update_jacobian! +end + + function ClimaLand.get_drivers(model::CanopyModel) return (model.atmos, model.radiation) end diff --git a/src/standalone/Vegetation/canopy_boundary_fluxes.jl b/src/standalone/Vegetation/canopy_boundary_fluxes.jl index 9ee9518227..341d96987f 100644 --- a/src/standalone/Vegetation/canopy_boundary_fluxes.jl +++ b/src/standalone/Vegetation/canopy_boundary_fluxes.jl @@ -173,7 +173,8 @@ function canopy_boundary_fluxes!( shf .= canopy_tf.shf lhf .= canopy_tf.lhf r_ae .= canopy_tf.r_ae - + p.canopy.energy.∂LHF∂qc .= canopy_tf.∂LHF∂qc + p.canopy.energy.∂SHF∂Tc .= canopy_tf.∂SHF∂Tc # Transpiration is per unit ground area, not leaf area (mult by LAI) fa.:($i_end) .= PlantHydraulics.transpiration_per_ground_area( canopy.hydraulics.transpiration, @@ -334,7 +335,7 @@ function canopy_turbulent_fluxes_at_a_point( T_in::FT = Thermodynamics.air_temperature(thermo_params, ts_in) ΔT = T_in - T_sfc r_ae::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(sc)) - ρ_air::FT = Thermodynamics.air_density(thermo_params, ts_in) + ρ_air::FT = Thermodynamics.air_density(thermo_params, ts_sfc) ustar::FT = conditions.ustar r_b::FT = FT(1 / 0.01 * (ustar / 0.04)^(-1 / 2)) # CLM 5, tech note Equation 5.122 leaf_r_b = r_b / LAI @@ -344,5 +345,14 @@ function canopy_turbulent_fluxes_at_a_point( Ẽ = E / _ρ_liq H = -ρ_air * cp_m * ΔT / (canopy_r_b + r_ae) # CLM 5, tech note Equation 5.88, setting H_v = H and solving to remove T_s LH = _LH_v0 * E - return (lhf = LH, shf = H, vapor_flux = Ẽ, r_ae = r_ae) + ∂LHF∂qc = ρ_air * _LH_v0 / (leaf_r_b + leaf_r_stomata + r_ae) + ∂SHF∂Tc = ρ_air * cp_m / (canopy_r_b + r_ae) + return ( + lhf = LH, + shf = H, + vapor_flux = Ẽ, + r_ae = r_ae, + ∂LHF∂qc = ∂LHF∂qc, + ∂SHF∂Tc = ∂SHF∂Tc, + ) end diff --git a/src/standalone/Vegetation/canopy_energy.jl b/src/standalone/Vegetation/canopy_energy.jl index c3950225c7..21db571f0a 100644 --- a/src/standalone/Vegetation/canopy_energy.jl +++ b/src/standalone/Vegetation/canopy_energy.jl @@ -10,12 +10,28 @@ abstract type AbstractCanopyEnergyModel{FT} <: AbstractCanopyComponent{FT} end ClimaLand.name(model::AbstractCanopyEnergyModel) = :energy -ClimaLand.auxiliary_vars(model::AbstractCanopyEnergyModel) = - (:shf, :lhf, :fa_energy_roots, :r_ae) +ClimaLand.auxiliary_vars(model::AbstractCanopyEnergyModel) = ( + :shf, + :lhf, + :fa_energy_roots, + :r_ae, + :∂SHF∂Tc, + :∂LHF∂qc, + :∂LW_n∂Tc, + :∂qc∂Tc, +) ClimaLand.auxiliary_types(model::AbstractCanopyEnergyModel{FT}) where {FT} = - (FT, FT, FT, FT) -ClimaLand.auxiliary_domain_names(model::AbstractCanopyEnergyModel) = - (:surface, :surface, :surface, :surface) + (FT, FT, FT, FT, FT, FT, FT, FT) +ClimaLand.auxiliary_domain_names(model::AbstractCanopyEnergyModel) = ( + :surface, + :surface, + :surface, + :surface, + :surface, + :surface, + :surface, + :surface, +) """ PrescribedCanopyTempModel{FT} <: AbstractCanopyEnergyModel{FT} @@ -87,11 +103,11 @@ where the canopy temperature is modeled prognostically. canopy_temperature(model::BigLeafEnergyModel, canopy, Y, p, t) = Y.canopy.energy.T -function make_compute_exp_tendency( +function make_compute_imp_tendency( model::BigLeafEnergyModel{FT}, canopy, ) where {FT} - function compute_exp_tendency!(dY, Y, p, t) + function compute_imp_tendency!(dY, Y, p, t) area_index = p.canopy.hydraulics.area_index ac_canopy = model.parameters.ac_canopy # Energy Equation: @@ -117,7 +133,7 @@ function make_compute_exp_tendency( p.canopy.energy.lhf - p.canopy.energy.fa_energy_roots ) / c_per_ground_area end - return compute_exp_tendency! + return compute_imp_tendency! end """ @@ -154,3 +170,63 @@ function root_energy_flux_per_ground_area!( ) where {FT} fa_energy .= FT(0) end + +function ClimaLand.make_compute_jacobian( + model::BigLeafEnergyModel{FT}, + canopy, +) where {FT} + function compute_jacobian!(jacobian::ImplicitEquationJacobian, Y, p, dtγ, t) + (; matrix) = jacobian + + # The derivative of the residual with respect to the prognostic variable + ∂Tres∂T = matrix[@name(canopy.energy.T), @name(canopy.energy.T)] + ∂LHF∂qc = p.canopy.energy.∂LHF∂qc + ∂SHF∂Tc = p.canopy.energy.∂SHF∂Tc + ∂LW_n∂Tc = p.canopy.energy.∂LW_n∂Tc + ∂qc∂Tc = p.canopy.energy.∂qc∂Tc + ϵ_c = p.canopy.radiative_transfer.ϵ + area_index = p.canopy.hydraulics.area_index + ac_canopy = model.parameters.ac_canopy + earth_param_set = canopy.parameters.earth_param_set + _σ = FT(LP.Stefan(earth_param_set)) + @. ∂LW_n∂Tc = -2 * 4 * _σ * ϵ_c * Y.canopy.energy.T^3 # ≈ ϵ_ground = 1 + @. ∂qc∂Tc = partial_q_sat_partial_T_liq(p.drivers.P, Y.canopy.energy.T)# use atmos air pressure as approximation for surface air pressure + @. ∂Tres∂T = + dtγ * MatrixFields.DiagonalMatrixRow( + (∂LW_n∂Tc - ∂SHF∂Tc - ∂LHF∂qc * ∂qc∂Tc) / + (ac_canopy * max(area_index.leaf + area_index.stem, eps(FT))), + ) - (I,) + end + return compute_jacobian! +end + +""" + partial_q_sat_partial_T_liq(P::FT, T::FT) where {FT} + +Computes the quantity ∂q_sat∂T at temperature T and pressure P, +over liquid water. + +Uses the polynomial approximation from Flatau et al. (1992). +""" +function partial_q_sat_partial_T_liq(P::FT, T::FT) where {FT} + esat = FT( + 6.11213476e2 + + 4.44007856e1 * T + + 1.43064234 * T^2 + + 2.64461437e-2 * T^3 + + 3.05903558e-4 * T^4 + + 1.96237241e-6 * T^5 + + 8.92344772e-9 * T^6 - 3.73208410e-11 * T^7 + 2.09339997e-14 * T^8, + ) + desatdT = FT( + 4.44017302e1 + + 2.86064092 * T + + 7.94683137e-2 * T^2 + + 1.21211669e-3 * T^3 + + 1.03354611e-5 * T^4 + + 4.04125005e-8 * T^5 - 7.88037859e-11 * T^6 - 1.14596802e-12 * T^7 + + 3.81294516e-15 * T^8, + ) + + return FT(0.622) * P / (P - FT(0.378) * esat)^2 * desatdT +end diff --git a/src/standalone/Vegetation/component_models.jl b/src/standalone/Vegetation/component_models.jl index 3697e0ce7c..2b95ec0c79 100644 --- a/src/standalone/Vegetation/component_models.jl +++ b/src/standalone/Vegetation/component_models.jl @@ -123,11 +123,48 @@ may be needed and passed in (via the `canopy` model itself). The right hand side for the entire canopy model can make use of these functions for the individual components. """ -function ClimaLand.make_compute_exp_tendency(::AbstractCanopyComponent, canopy) - function compute_exp_tendency!(dY, Y, p, t) end +function ClimaLand.make_compute_exp_tendency( + component::AbstractCanopyComponent, + canopy, +) + function compute_exp_tendency!(dY, Y, p, t) + vars = prognostic_vars(component) + if vars != () + getproperty(getproperty(dY, name(canopy)), name(component)) .= 0 + end + end return compute_exp_tendency! end +""" + ClimaLand.make_compute_exp_tendency(component::AbstractCanopyComponent, canopy) + +Creates the compute_imp_tendency!(dY,Y,p,t) function for the canopy `component`. + +Since component models are not standalone models, other information +may be needed and passed in (via the `canopy` model itself). +The right hand side for the entire canopy model can make use of +these functions for the individual components. +""" +function ClimaLand.make_compute_imp_tendency( + component::AbstractCanopyComponent, + canopy, +) + function compute_imp_tendency!(dY, Y, p, t) + vars = prognostic_vars(component) + if vars != () + getproperty(getproperty(dY, name(canopy)), name(component)) .= 0 + end + + end + return compute_imp_tendency! +end + + +function make_update_jacobian(component::AbstractCanopyComponent, canopy) + function update_jacobian!(W, Y, p, dtγ, t) end + return update_jacobian! +end """ set_canopy_prescribed_field!(component::AbstractCanopyComponent, diff --git a/test/standalone/Vegetation/canopy_model.jl b/test/standalone/Vegetation/canopy_model.jl index ceb3aa1d17..4c77bb9850 100644 --- a/test/standalone/Vegetation/canopy_model.jl +++ b/test/standalone/Vegetation/canopy_model.jl @@ -719,6 +719,14 @@ for FT in (Float32, Float64) t0 = FT(0.0) set_initial_cache!(p, Y, t0) exp_tendency! = make_exp_tendency(canopy) + imp_tendency! = ClimaLand.make_imp_tendency(canopy) + tendency_jacobian! = ClimaLand.make_tendency_jacobian(canopy) + # set up jacobian info + jac_kwargs = (; + jac_prototype = ImplicitEquationJacobian(Y), + Wfact = tendency_jacobian!, + ) + dY = similar(Y) exp_tendency!(dY, Y, p, t0) turb_fluxes = ClimaLand.Canopy.canopy_turbulent_fluxes(