diff --git a/Manifest.toml b/Manifest.toml index 94e39d56e..3a78e27c1 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,18 +1,19 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.9.2" +julia_version = "1.9.3" manifest_format = "2.0" -project_hash = "1df4fb4dcc221fff4ae2f940474ce7322fc6228d" +project_hash = "b8336a615981d764881d88f0d8193dd7f4667beb" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] -git-tree-sha1 = "cad4c758c0038eea30394b1b671526921ca85b21" +git-tree-sha1 = "d92ad398961a3ed262d8bf04a1a2b8340f915fef" uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" -version = "1.4.0" -weakdeps = ["ChainRulesCore"] +version = "1.5.0" +weakdeps = ["ChainRulesCore", "Test"] [deps.AbstractFFTs.extensions] AbstractFFTsChainRulesCoreExt = "ChainRulesCore" + AbstractFFTsTestExt = "Test" [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"] @@ -87,9 +88,9 @@ version = "0.1.2" [[deps.CUDA]] deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "Preferences", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "SpecialFunctions", "UnsafeAtomicsLLVM"] -git-tree-sha1 = "35160ef0f03b14768abfd68b830f8e3940e8e0dc" +git-tree-sha1 = "968c1365e2992824c3e7a794e30907483f8469a9" uuid = "052768ef-5323-5732-b1bb-66c8b64840ba" -version = "4.4.0" +version = "4.4.1" [[deps.CUDA_Driver_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] @@ -128,9 +129,9 @@ version = "0.2.4" [[deps.Compat]] deps = ["UUIDs"] -git-tree-sha1 = "4e88377ae7ebeaf29a047aa1ee40826e0b708a5d" +git-tree-sha1 = "e460f044ca8b99be31d35fe54fc33a5c33dd8ed7" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.7.0" +version = "4.9.0" weakdeps = ["Dates", "LinearAlgebra"] [deps.Compat.extensions] @@ -143,9 +144,9 @@ version = "1.0.5+0" [[deps.ConstructionBase]] deps = ["LinearAlgebra"] -git-tree-sha1 = "738fec4d684a9a6ee9598a8bfee305b26831f28c" +git-tree-sha1 = "fe2838a593b5f776e1597e086dcd47560d94e816" uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" -version = "1.5.2" +version = "1.5.3" [deps.ConstructionBase.extensions] ConstructionBaseIntervalSetsExt = "IntervalSets" @@ -173,9 +174,9 @@ version = "1.15.0" [[deps.DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "cf25ccb972fec4e4817764d01c82386ae94f77b4" +git-tree-sha1 = "3dbd312d370723b6bb43ba9d02fc36abade4518d" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.14" +version = "0.18.15" [[deps.DataValueInterfaces]] git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" @@ -186,6 +187,16 @@ version = "1.0.0" deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +[[deps.Distances]] +deps = ["LinearAlgebra", "Statistics", "StatsAPI"] +git-tree-sha1 = "b6def76ffad15143924a2199f72a5cd883a2e8a9" +uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +version = "0.10.9" +weakdeps = ["SparseArrays"] + + [deps.Distances.extensions] + DistancesSparseArraysExt = "SparseArrays" + [[deps.Distributed]] deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" @@ -207,9 +218,9 @@ uuid = "b305315f-e792-5b7a-8f41-49f472929428" version = "1.0.1" [[deps.ExprTools]] -git-tree-sha1 = "c1d06d129da9f55715c6c212866f5b1bddc5fa00" +git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec" uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" -version = "0.1.9" +version = "0.1.10" [[deps.FFTW]] deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] @@ -261,9 +272,9 @@ version = "1.3.1" [[deps.HDF5_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"] -git-tree-sha1 = "3b20c3ce9c14aedd0adca2bc8c882927844bd53d" +git-tree-sha1 = "10c72358aaaa5cd6bc7cc39b95e6eadf92f5a336" uuid = "0234f1f7-429e-5d53-9886-15a909be8d59" -version = "1.14.0+0" +version = "1.14.2+0" [[deps.IfElse]] git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" @@ -278,9 +289,9 @@ version = "0.2.1" [[deps.IntelOpenMP_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "0cb9352ef2e01574eeebdb102948a58740dcaf83" +git-tree-sha1 = "ad37c091f7d7daf900963171600d7c1c5c3ede32" uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" -version = "2023.1.0+0" +version = "2023.2.0+0" [[deps.InteractiveUtils]] deps = ["Markdown"] @@ -304,27 +315,27 @@ version = "1.0.0" [[deps.JLD2]] deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "Printf", "Reexport", "Requires", "TranscodingStreams", "UUIDs"] -git-tree-sha1 = "5df8278ad24772c0c6dbbeb97b162ccf29ced2a9" +git-tree-sha1 = "aa6ffef1fd85657f4999030c52eaeec22a279738" uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -version = "0.4.32" +version = "0.4.33" [[deps.JLLWrappers]] -deps = ["Preferences"] -git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" +deps = ["Artifacts", "Preferences"] +git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca" uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.4.1" +version = "1.5.0" [[deps.JSON3]] deps = ["Dates", "Mmap", "Parsers", "PrecompileTools", "StructTypes", "UUIDs"] -git-tree-sha1 = "5b62d93f2582b09e469b3099d839c2d2ebf5066d" +git-tree-sha1 = "95220473901735a0f4df9d1ca5b171b568b2daa3" uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" -version = "1.13.1" +version = "1.13.2" [[deps.KernelAbstractions]] deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "Requires", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"] -git-tree-sha1 = "6d08ca80b621635fed9cdfeb9a4280a574020bf3" +git-tree-sha1 = "4c5875e4c228247e1c2b087669846941fb6e0118" uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c" -version = "0.9.7" +version = "0.9.8" [deps.KernelAbstractions.extensions] EnzymeExt = "EnzymeCore" @@ -334,15 +345,15 @@ version = "0.9.7" [[deps.LLVM]] deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"] -git-tree-sha1 = "8695a49bfe05a2dc0feeefd06b4ca6361a018729" +git-tree-sha1 = "a9d2ce1d5007b1e8f6c5b89c5a31ff8bd146db5c" uuid = "929cbde3-209d-540e-8aea-75f648917ca0" -version = "6.1.0" +version = "6.2.1" [[deps.LLVMExtra_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] -git-tree-sha1 = "c35203c1e1002747da220ffc3c0762ce7754b08c" +git-tree-sha1 = "7ca6850ae880cc99b59b88517545f91a52020afa" uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab" -version = "0.0.23+0" +version = "0.0.25+0" [[deps.LLVMOpenMP_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -377,10 +388,10 @@ version = "1.10.2+0" uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" [[deps.Libiconv_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "c7cb1f5d892775ba13767a87c7ada0b980ea0a71" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175" uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" -version = "1.16.1+2" +version = "1.17.0+0" [[deps.LinearAlgebra]] deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] @@ -388,9 +399,9 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[deps.LogExpFunctions]] deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "c3ce8e7420b3a6e071e0fe4745f5d4300e37b13f" +git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa" uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.24" +version = "0.3.26" [deps.LogExpFunctions.extensions] LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" @@ -407,15 +418,15 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[deps.MKL_jll]] deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] -git-tree-sha1 = "154d7aaa82d24db6d8f7e4ffcfe596f40bff214b" +git-tree-sha1 = "eb006abbd7041c28e0d16260e50a24f8f9104913" uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" -version = "2023.1.0+0" +version = "2023.2.0+0" [[deps.MPI]] deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "PkgVersion", "PrecompileTools", "Requires", "Serialization", "Sockets"] -git-tree-sha1 = "98fc280a56d3b5cc218435f82df60e05771fefa6" +git-tree-sha1 = "df53d0e1e0dbebf2315f4cd35e13e52ad43416c2" uuid = "da04e1cc-30fd-572f-bb4f-1f8673147195" -version = "0.20.12" +version = "0.20.15" [deps.MPI.extensions] AMDGPUExt = "AMDGPU" @@ -433,9 +444,9 @@ version = "4.1.2+0" [[deps.MPIPreferences]] deps = ["Libdl", "Preferences"] -git-tree-sha1 = "d86a788b336e8ae96429c0c42740ccd60ac0dfcc" +git-tree-sha1 = "781916a2ebf2841467cda03b6f1af43e23839d85" uuid = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267" -version = "0.1.8" +version = "0.1.9" [[deps.MPItrampoline_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] @@ -445,9 +456,9 @@ version = "5.3.1+0" [[deps.MacroTools]] deps = ["Markdown", "Random"] -git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2" +git-tree-sha1 = "9ee1618cbf5240e6d4e0371d6f24065083f60c48" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.10" +version = "0.5.11" [[deps.Markdown]] deps = ["Base64"] @@ -488,10 +499,10 @@ uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" version = "1.2.0" [[deps.Oceananigans]] -deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"] -git-tree-sha1 = "c5120ad3b0d67ab23e9401063ecb9c14999f2613" +deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"] +git-tree-sha1 = "8746c619e15950de59160b2d99e61961cfa57fcc" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" -version = "0.85.0" +version = "0.87.2" [[deps.OffsetArrays]] deps = ["Adapt"] @@ -517,9 +528,9 @@ version = "4.1.5+0" [[deps.OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "cae3153c7f6cf3f069a853883fd1919a6e5bab5b" +git-tree-sha1 = "e78db7bd5c26fc5a6911b50a47ee302219157ea8" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.0.9+0" +version = "3.0.10+0" [[deps.OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] @@ -528,21 +539,21 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" version = "0.5.5+0" [[deps.OrderedCollections]] -git-tree-sha1 = "d321bf2de576bf25ec4d3e4360faca399afca282" +git-tree-sha1 = "2e73fe17cac3c62ad1aebe70d44c963c3cfdc3e3" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.6.0" +version = "1.6.2" [[deps.Parsers]] deps = ["Dates", "PrecompileTools", "UUIDs"] -git-tree-sha1 = "4b2e829ee66d4218e0cef22c0a64ee37cf258c29" +git-tree-sha1 = "716e24b21538abc91f6205fd1d8363f39b442851" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.7.1" +version = "2.7.2" [[deps.PencilArrays]] deps = ["Adapt", "JSON3", "LinearAlgebra", "MPI", "OffsetArrays", "Random", "Reexport", "StaticArrayInterface", "StaticArrays", "StaticPermutations", "Strided", "TimerOutputs", "VersionParsing"] -git-tree-sha1 = "c1d2b659ce423897de45e998f49f77e789e1859d" +git-tree-sha1 = "c30a7fb1e424ea572962dac493ad6c3f556695e0" uuid = "0e08944d-e94e-41b1-9406-dcf66b6a9d2e" -version = "0.18.1" +version = "0.19.0" [deps.PencilArrays.extensions] PencilArraysDiffEqExt = ["DiffEqBase"] @@ -554,9 +565,9 @@ version = "0.18.1" [[deps.PencilFFTs]] deps = ["AbstractFFTs", "FFTW", "LinearAlgebra", "MPI", "PencilArrays", "Reexport", "TimerOutputs"] -git-tree-sha1 = "af200cef52069f3428127b237919d7c69eb6b10d" +git-tree-sha1 = "bd69f3f0ee248cfb4241800aefb705b5ded592ff" uuid = "4a48f351-57a6-4416-9ec4-c37015456aae" -version = "0.15.0" +version = "0.15.1" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] @@ -565,15 +576,15 @@ version = "1.9.2" [[deps.PkgVersion]] deps = ["Pkg"] -git-tree-sha1 = "f6cf8e7944e50901594838951729a1861e668cb8" +git-tree-sha1 = "f9501cc0430a26bc3d156ae1b5b0c1b47af4d6da" uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688" -version = "0.3.2" +version = "0.3.3" [[deps.PrecompileTools]] deps = ["Preferences"] -git-tree-sha1 = "9673d39decc5feece56ef3940e5dafba15ba0f81" +git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f" uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -version = "1.1.2" +version = "1.2.0" [[deps.Preferences]] deps = ["TOML"] @@ -587,9 +598,9 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" [[deps.ProgressBars]] deps = ["Printf"] -git-tree-sha1 = "9d84c8646109eb8bc7a006d59b157c64d5155c81" +git-tree-sha1 = "b437cdb0385ed38312d91d9c00c20f3798b30256" uuid = "49802e3a-d2f1-5c88-81d8-b72133a6f568" -version = "1.5.0" +version = "1.5.1" [[deps.Quaternions]] deps = ["LinearAlgebra", "Random", "RealDot"] @@ -642,9 +653,9 @@ version = "1.3.0" [[deps.Roots]] deps = ["ChainRulesCore", "CommonSolve", "Printf", "Setfield"] -git-tree-sha1 = "de432823e8aab4dd1a985be4be768f95acf152d4" +git-tree-sha1 = "ff42754a57bb0d6dcfe302fd0d4272853190421f" uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" -version = "2.0.17" +version = "2.0.19" [deps.Roots.extensions] RootsForwardDiffExt = "ForwardDiff" @@ -686,12 +697,6 @@ git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac" uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" version = "1.1.1" -[[deps.SnoopPrecompile]] -deps = ["Preferences"] -git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c" -uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c" -version = "1.0.3" - [[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" @@ -701,9 +706,9 @@ uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[deps.SpecialFunctions]] deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "7beb031cf8145577fbccacd94b8a8f4ce78428d3" +git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.3.0" +version = "2.3.1" weakdeps = ["ChainRulesCore"] [deps.SpecialFunctions.extensions] @@ -711,15 +716,15 @@ weakdeps = ["ChainRulesCore"] [[deps.Static]] deps = ["IfElse"] -git-tree-sha1 = "dbde6766fc677423598138a5951269432b0fcc90" +git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.8.7" +version = "0.8.8" [[deps.StaticArrayInterface]] -deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "Requires", "SnoopPrecompile", "SparseArrays", "Static", "SuiteSparse"] -git-tree-sha1 = "33040351d2403b84afce74dae2e22d3f5b18edcb" +deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"] +git-tree-sha1 = "03fec6800a986d191f64f5c0996b59ed526eda25" uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718" -version = "1.4.0" +version = "1.4.1" weakdeps = ["OffsetArrays", "StaticArrays"] [deps.StaticArrayInterface.extensions] @@ -728,18 +733,18 @@ weakdeps = ["OffsetArrays", "StaticArrays"] [[deps.StaticArrays]] deps = ["LinearAlgebra", "Random", "StaticArraysCore"] -git-tree-sha1 = "fffc14c695c17bfdbfa92a2a01836cdc542a1e46" +git-tree-sha1 = "9cabadf6e7cd2349b6cf49f1915ad2028d65e881" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.6.1" +version = "1.6.2" weakdeps = ["Statistics"] [deps.StaticArrays.extensions] StaticArraysStatisticsExt = "Statistics" [[deps.StaticArraysCore]] -git-tree-sha1 = "1d5708d926c76a505052d0d24a846d5da08bc3a4" +git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d" uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" -version = "1.4.1" +version = "1.4.2" [[deps.StaticPermutations]] git-tree-sha1 = "193c3daa18ff3e55c1dae66acb6a762c4a3bdb0b" @@ -751,11 +756,17 @@ deps = ["LinearAlgebra", "SparseArrays"] uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" version = "1.9.0" +[[deps.StatsAPI]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "1ff449ad350c9c4cbc756624d6f8a8c3ef56d3ed" +uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" +version = "1.7.0" + [[deps.Strided]] deps = ["LinearAlgebra", "StridedViews", "TupleTools"] -git-tree-sha1 = "b32eadf6ac726a790567fdc872b63117712e16a8" +git-tree-sha1 = "137303f5e0a39f966b462c53ae2c5c6e34c4828b" uuid = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" -version = "2.0.1" +version = "2.0.3" [[deps.StridedViews]] deps = ["LinearAlgebra"] @@ -863,10 +874,10 @@ uuid = "81def892-9a0e-5fdd-b105-ffc91e053289" version = "1.3.0" [[deps.XML2_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "93c41695bc1c08c46c5899f4fe06d6ead504bb73" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"] +git-tree-sha1 = "04a51d15436a572301b5abbb9d099713327e9fc4" uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" -version = "2.10.3+0" +version = "2.10.4+0" [[deps.Zlib_jll]] deps = ["Libdl"] diff --git a/Project.toml b/Project.toml index a8c3ef66f..4761ddc56 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OceanBioME" uuid = "a49af516-9db8-4be4-be45-1dad61c5a376" authors = ["Jago Strong-Wright ", "John R Taylor ", "Si Chen "] -version = "0.5.0" +version = "0.6.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" @@ -12,12 +12,12 @@ Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" [compat] -Adapt = "3.4" -JLD2 = "0.4.31" +Adapt = "3" +JLD2 = "0.4" KernelAbstractions = "0.9" -Oceananigans = "0.84.1, 0.85" +Oceananigans = "0.84.1, 0.85, 0.86, 0.87" Roots = "2" -StructArrays = "0.6.15" +StructArrays = "0.4, 0.5, 0.6" julia = "1.9" [extras] diff --git a/docs/Project.toml b/docs/Project.toml index a6258f013..a85e13173 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -17,4 +17,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" CairoMakie = "0.10" DocumenterCitations = "1" Literate = "≥2.9.0" -Oceananigans = "0.84.1, 0.85" +Oceananigans = "0.84.1, 0.85, 0.86, 0.87" diff --git a/docs/src/model_components/biogeochemical/LOBSTER.md b/docs/src/model_components/biogeochemical/LOBSTER.md index d5686452b..848a57c73 100644 --- a/docs/src/model_components/biogeochemical/LOBSTER.md +++ b/docs/src/model_components/biogeochemical/LOBSTER.md @@ -13,15 +13,16 @@ julia> grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200)); julia> bgc_model = LOBSTER(; grid, carbonates = true) Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model (Float64) - Light Attenuation Model: - └── Two-band light attenuation model (Float64) Optional components: ├── Carbonates ✅ ├── Oxygen ❌ └── Variable Redfield Ratio ❌ Sinking Velocities: ├── sPOM: 0.0 to -3.47e-5 m/s - └── bPOM: 0.0 to -0.0023148148148148147 m/s + └── bPOM: 0.0 to -0.0023148148148148147 m/s + Light attenuation: Two-band light attenuation model (Float64) + Sediment: Nothing + Particles: Nothing ``` ## Model equations diff --git a/docs/src/model_components/biogeochemical/index.md b/docs/src/model_components/biogeochemical/index.md index b3a74f495..2e28f7f0a 100644 --- a/docs/src/model_components/biogeochemical/index.md +++ b/docs/src/model_components/biogeochemical/index.md @@ -5,7 +5,7 @@ Biogeochemical (BGC) models can be used within the [Oceananigans biogeochemistry For details of the BGC models currently implemented please see the following pages. ## Oceananigans setup -At the simplest level all that is required to setup a BGC model is to pass it to the Oceananigans model setup: +At the simplest level all that is required to setup an existing OceanBioME BGC model is to pass it to the Oceananigans model setup: ```julia model = NonhydrostaticModel(; grid, ..., @@ -19,3 +19,7 @@ MODEL_NAME(; grid, growth_rate = 10.0) This will set up the required tracers and auxiliary fields, and you may also set boundary conditions or additional forcing through the normal Oceananigans setup. Models usually have a default [light attenuation model](@ref light) specified, these may be substituted easily by passing different models as parameters as above. + +Our models are implemented in an abstract framework `Biogeochemistry` which contains `underlying_biogeochemistry`, `light_attenuation`, `sediment`, and `modifiers`. +This is automatically set up for existing BGC models, but may also be used to couple any BGC model with light attenuation and sediments. +See the [implementation page](@ref model_implementation) for some more information on how to couple other models. \ No newline at end of file diff --git a/docs/src/model_components/individuals/index.md b/docs/src/model_components/individuals/index.md index 95b39dc3b..be71f1e83 100644 --- a/docs/src/model_components/individuals/index.md +++ b/docs/src/model_components/individuals/index.md @@ -43,7 +43,7 @@ function update_lagrangian_particle_properties!(particles::GrowingParticles, mod return nothing end -nothing # hide +nothing #hide ``` In this example the particles will not move around, and are only integrated on a single thread. For a more comprehensive example see the [Sugar Kelp](@ref SLatissima) implementation. We then need to update the tracer tendencies to match the nutrients' uptake: @@ -67,7 +67,7 @@ function update_tendencies!(bgc, particles::GrowingParticles, model) return nothing end -nothing # hide +nothing #hide ``` Now we can just plug this into any biogeochemical model setup to have particles (currently [NPZD](@ref NPZD) and [LOBSTER](@ref LOBSTER)): diff --git a/docs/src/model_components/utils.md b/docs/src/model_components/utils.md index 8f542a928..1adfd85dc 100644 --- a/docs/src/model_components/utils.md +++ b/docs/src/model_components/utils.md @@ -22,20 +22,17 @@ simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) ## Negative tracer detection As a temporary measure we have implemented a callback to either detect negative tracers and either scale a conserved group, force them back to zero, or throw an error. Please see the numerical implementations' page for details. This can be set up by: ```julia -negativity_protection! = ScaleNegativeTracers(tracers = (:P, :Z, :N)) -simulation.callbacks[:neg] = Callback(negativity_protection!; callsite = UpdateStateCallsite()) +negativity_protection = ScaleNegativeTracers((:P, :Z, :N)) +biogeochemistry = Biogeochemistry(...; modifiers = negativity_protection) ``` You may also pass a scale factor for each component (e.g. in case they have different redfield ratios): ```julia -negativity_protection! = ScaleNegativeTracers(tracers = (:P, :Z, :D), scalefactors = (P = 1, Z = 1, D = 2)) -simulation.callbacks[:neg] = Callback(negativity_protection!; callsite = UpdateStateCallsite()) +negativity_protection = ScaleNegativeTracers((:P, :Z, :N); scalefactors = (1, 1, 2)) +biogeochemistry = Biogeochemistry(...; modifiers = negativity_protection) ``` Here you should carefully consider which tracers form a conserved group (if at all). Alternatively, force to zero by: ```julia -simulation.callbacks[:neg] = Callback(OceanBioME.no_negative_tracers!, callsite = UpdateStateCallsite()) +negativity_protection = ZeroNegativeTracers() +biogeochemistry = Biogeochemistry(...; modifiers = negativity_protection) ``` -or throw an error: -```julia -simulation.callbacks[:neg] = Callback(OceanBioME.error_on_neg!, callsite = UpdateStateCallsite()) -``` -The latter two both optionally take a named tuple of parameters which may include `exclude` which can be a tuple of tracer names (Symbols) which are allowed to be negative. +The latter optionally takes a named tuple of parameters that may include `exclude`, which can be a tuple of tracer names (Symbols) which are allowed to be negative. diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index 1d3d864ac..59fcd5185 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -18,20 +18,19 @@ For this example we are going to implement the simple Nutrient-Phytoplankton mod The first step is to import the abstract type from OceanBioME, some units from Oceananigans (for ease of parameter definition), and [`import`](https://stackoverflow.com/questions/27086159/what-is-the-difference-between-using-and-import-in-julia-when-building-a-mod) some functions from Oceananigans in order to add methods to: ```@example implementing -using OceanBioME: ContinuousFormBiogeochemistry +using OceanBioME: Biogeochemistry +using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry using Oceananigans.Units import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, required_biogeochemical_auxiliary_fields, - update_biogeochemical_state!, - biogeochemical_drift_velocity, - biogeochemical_auxiliary_fields + biogeochemical_drift_velocity ``` We then define our `struct` with the model parameters, as well as slots for the particles, light attenuation, and sediment models: ```@example implementing -@kwdef struct NutrientPhytoplankton{FT, LA, S, W, P} <:ContinuousFormBiogeochemistry{LA, S, P} +@kwdef struct NutrientPhytoplankton{FT, W} <: AbstractContinuousFormBiogeochemistry base_growth_rate :: FT = 1.27 / day # 1 / seconds nutrient_half_saturation :: FT = 0.025 * 1000 / 14 # mmol N / m³ light_half_saturation :: FT = 300.0 # micro einstein / m² / s @@ -41,10 +40,7 @@ We then define our `struct` with the model parameters, as well as slots for the mortality_rate :: FT = 0.15 / day # 1 / seconds crowding_mortality_rate :: FT = 0.004 / day / 1000 * 14 # 1 / seconds / mmol N / m³ - light_attenuation_model :: LA = nothing - sediment_model :: S = nothing sinking_velocity :: W = 2 / day - particles :: P = nothing end ``` @@ -54,8 +50,6 @@ Here, we use descriptive names for the parameters. Below, each of these paramete required_biogeochemical_tracers(::NutrientPhytoplankton) = (:N, :P, :T) required_biogeochemical_auxiliary_fields(::NutrientPhytoplankton) = (:PAR, ) - -biogeochemical_auxiliary_fields(bgc::NutrientPhytoplankton) = biogeochemical_auxiliary_fields(bgc.light_attenuation_model) ``` Next, we define the functions that specify how the phytoplankton ``P`` evolve. In the absence of advection and diffusion (both of which are handled by Oceananigans), we want the phytoplankton to evolve at the rate given by: @@ -126,10 +120,6 @@ Finally, we need to define some functions to allow us to update the time-depende using OceanBioME: BoxModel import OceanBioME.BoxModels: update_boxmodel_state! -function update_biogeochemical_state!(bgc::NutrientPhytoplankton, model) - update_PAR!(model, bgc.light_attenuation_model) -end - function update_boxmodel_state!(model::BoxModel{<:NutrientPhytoplankton, <:Any, <:Any, <:Any, <:Any, <:Any}) getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time) getproperty(model.values, :T) .= model.forcing.T(model.clock.time) @@ -228,7 +218,7 @@ import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisa @inline nitrogen_flux(i, j, k, grid, advection, bgc::NutrientPhytoplankton, tracers) = sinking_flux(i, j, k, grid, advection, Val(:P), bgc, tracers) -@inline carbon_flux(bgc::NutrientPhytoplankton, tracers, i, j) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * 6.56 +@inline carbon_flux(i, j, k, grid, advection, bgc::NutrientPhytoplankton, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * 6.56 @inline remineralisation_receiver(::NutrientPhytoplankton) = :N ``` @@ -254,9 +244,9 @@ grid = RectilinearGrid(topology = (Flat, Flat, Bounded), size = (32, ), extent = # setup the biogeochemical model -light_attenuation_model = TwoBandPhotosyntheticallyActiveRadiation(; grid, surface_PAR) +light_attenuation = TwoBandPhotosyntheticallyActiveRadiation(; grid, surface_PAR) -sediment_model = InstantRemineralisation(; grid) +sediment = InstantRemineralisation(; grid) sinking_velocity = ZFaceField(grid) @@ -264,7 +254,12 @@ w_sink(x, y, z) = 2 / day * tanh(z / 5) set!(sinking_velocity, w_sink) -biogeochemistry = NutrientPhytoplankton(; light_attenuation_model, sinking_velocity, sediment_model) +negative_tracer_scaling = ScaleNegativeTracers((:N, :P)) + +biogeochemistry = Biogeochemistry(NutrientPhytoplankton(; sinking_velocity); + light_attenuation, + sediment, + modifiers = negative_tracer_scaling) # put the model together @@ -284,15 +279,12 @@ simulation.output_writers[:tracers] = JLD2OutputWriter(model, model.tracers, schedule = TimeInterval(1day), overwrite_existing = true) -simulation.output_writers[:sediment] = JLD2OutputWriter(model, model.biogeochemistry.sediment_model.fields, +simulation.output_writers[:sediment] = JLD2OutputWriter(model, model.biogeochemistry.sediment.fields, indices = (:, :, 1), filename = "column_np_sediment.jld2", schedule = TimeInterval(1day), overwrite_existing = true) -scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:N, :P)) -simulation.callbacks[:nan_tendencies] = Callback(remove_NaN_tendencies!; callsite = TendencyCallsite()) - run!(simulation) ``` diff --git a/docs/src/quick_start.md b/docs/src/quick_start.md index 54eba8114..b27e43a7a 100644 --- a/docs/src/quick_start.md +++ b/docs/src/quick_start.md @@ -3,7 +3,7 @@ OceanBioME provides biogeochemical models to plug into [Oceananigans](https://github.com/CliMA/Oceananigans.jl), for example this code will run one month of a single column, 7 variable (P, Z, sPOM, bPOM, DOM, NO₃, NH₄) biogeochemical situation with constant forcing. First we need to check we have the required dependencies: -```@example quickstart +```julia using Pkg Pkg.add(["OceanBioME", "Oceananigans"]) ``` @@ -30,11 +30,15 @@ simulation.output_writers[:profiles] = JLD2OutputWriter(model, model.tracers, run!(simulation) ``` -This isn't quite as simple as it could be as it records the output so that we can visualize it: +We can then visualize it, first check the required packages are installed: -```@example quickstart +```julia Pkg.add("CairoMakie") +``` +and then load the data and plot: + +```@example quickstart using CairoMakie phytoplankton = FieldTimeSeries("quickstart.jld2", "P") diff --git a/examples/column.jl b/examples/column.jl index b28f72f62..44acb5d6a 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -40,14 +40,22 @@ nothing #hide # Define the grid. grid = RectilinearGrid(size=(1, 1, 50), extent=(20meters, 20meters, 200meters)) -# Specify the boundary conditions for DIC and O₂ based on the air-sea CO₂ and O₂ flux + +# ## Model +# First we define the biogeochemical model including carbonate chemistry +# and scaling of negative tracers(see discussion in the [positivity preservation](@ref pos-preservation)) +# and then setup the Oceananigans model with the boundary condition for the DIC based on the air-sea CO₂ flux. + +biogeochemistry = LOBSTER(; grid, + surface_phytosynthetically_active_radiation = PAR⁰, + carbonates = true, + scale_negatives = true) + CO₂_flux = GasExchange(; gas = :CO₂, temperature = temp, salinity = (args...) -> 35) model = NonhydrostaticModel(; grid, - closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), - biogeochemistry = LOBSTER(; grid, - surface_phytosynthetically_active_radiation = PAR⁰, - carbonates = true), + closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), + biogeochemistry, boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), )) set!(model, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2239.8, Alk = 2409.0) @@ -56,15 +64,14 @@ set!(model, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2239.8, Alk = 2 # Next we setup a simulation and add some callbacks that: # - Show the progress of the simulation # - Store the model and particles output -# - Prevent the tracers from going negative from numerical error (see discussion in the [positivity preservation](@ref pos-preservation) implementation section) simulation = Simulation(model, Δt = 3minutes, stop_time = 100days) progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", - iteration(sim), - prettytime(sim), - prettytime(sim.Δt), - prettytime(sim.run_wall_time)) + iteration(sim), + prettytime(sim), + prettytime(sim.Δt), + prettytime(sim.run_wall_time)) simulation.callbacks[:progress] = Callback(progress_message, TimeInterval(10days)) @@ -73,9 +80,6 @@ simulation.output_writers[:profiles] = JLD2OutputWriter(model, merge(model.trace filename = "$filename.jld2", schedule = TimeInterval(1day), overwrite_existing = true) - -scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)) -simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) nothing #hide # ## Run! @@ -95,6 +99,7 @@ bPOM = FieldTimeSeries("$filename.jld2", "bPOM") x, y, z = nodes(P) times = P.times +nothing #hide # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and # the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth @@ -102,10 +107,12 @@ times = P.times air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) +using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity + for (i, t) in enumerate(times) air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) - carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.sPOM.w[1, 1, grid.Nz-20] + - bPOM[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.bPOM.w[1, 1, grid.Nz-20]) * model.biogeochemistry.organic_redfield + carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + + bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), model.biogeochemistry) end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 0e8774a43..506befc75 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -17,7 +17,7 @@ # ## Model setup # First load the required packages -using OceanBioME +using OceanBioME using Oceananigans, Random, Printf, NetCDF, Interpolations, DataDeps using Oceananigans.Units @@ -28,7 +28,7 @@ nothing #hide # Loading the forcing data from our online copy dd = DataDep( "example_data", - "example data from subpolar re analysis and observational products", + "example data from subpolar re-analysis and observational products", "https://github.com/OceanBioME/OceanBioME_example_data/raw/main/subpolar.nc" ) register(dd) @@ -39,9 +39,9 @@ salinity = ncread(filename, "so") mld = ncread(filename, "mld") par = ncread(filename, "par") -temperature_itp = LinearInterpolation(times, temp) -salinity_itp = LinearInterpolation(times, salinity) -mld_itp = LinearInterpolation(times, mld) +temperature_itp = LinearInterpolation(times, temp) +salinity_itp = LinearInterpolation(times, salinity) +mld_itp = LinearInterpolation(times, mld) PAR_itp = LinearInterpolation(times, par) t_function(x, y, z, t) = temperature_itp(mod(t, 364days)) @@ -65,12 +65,17 @@ grid = RectilinearGrid(size = (1, 1, Nz), x = (0, 20meters), y = (0, 20meters), # ## Biogeochemical and Oceananigans model # Here we instantiate the LOBSTER model with carbonate chemistry and a surface flux of DIC (CO₂) + +biogeochemistry = LOBSTER(; grid, + surface_phytosynthetically_active_radiation = surface_PAR, + carbonates = true, + scale_negatives = true) + CO₂_flux = GasExchange(; gas = :CO₂, temperature = t_function, salinity = s_function) + model = NonhydrostaticModel(; grid, - closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), - biogeochemistry = LOBSTER(; grid, - surface_phytosynthetically_active_radiation = surface_PAR, - carbonates = true), + closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), + biogeochemistry, boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux),)) set!(model, P = 0.03, Z = 0.03, NO₃ = 11.0, NH₄ = 0.05, DIC = 2200.0, Alk = 2400.0) @@ -99,15 +104,11 @@ simulation.output_writers[:profiles] = JLD2OutputWriter(model, schedule = TimeInterval(1day), overwrite_existing = true) -# TODO: make tendency callback to force no NaNs in tendencies - -scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)) -simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) - wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, - max_change = 2.0, min_change = 0.5, + max_change = 1.5, min_change = 0.75, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) + simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) nothing #hide @@ -128,17 +129,20 @@ bPOM = FieldTimeSeries("$filename.jld2", "bPOM") x, y, z = nodes(P) times = P.times +nothing #hide # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and -# the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth +# the carbon export by computing how much carbon sinks below some arbitrary depth; here we use depth # that corresponds to `k = grid.Nz - 20`. air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) +using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity + for (i, t) in enumerate(times) air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], t_function(1, 1, 0, t), s_function(1, 1, 0, t)) - carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.sPOM.w[1, 1, grid.Nz-20] + - bPOM[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.bPOM.w[1, 1, grid.Nz-20]) * model.biogeochemistry.organic_redfield + carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + + bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), model.biogeochemistry) end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. diff --git a/examples/eady.jl b/examples/eady.jl index 15a729dfe..18f06783e 100644 --- a/examples/eady.jl +++ b/examples/eady.jl @@ -119,12 +119,6 @@ simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.tracers filename = "eady_turbulence_bgc", overwrite_existing = true) -# Prevent the tracer values going negative - this is especially important in this model while no positivity -# preserving diffusion is implemented. -scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)) -simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) -simulation.callbacks[:nan_tendencies] = Callback(remove_NaN_tendencies!; callsite = TendencyCallsite()) -simulation.callbacks[:abort_zeros] = Callback(zero_negative_tracers!; callsite = UpdateStateCallsite()) nothing #hide # Run the simulation @@ -141,6 +135,7 @@ times = ζ.times xζ, yζ, zζ = nodes(ζ) xc, yc, zc = nodes(P) +nothing #hide # and plot. diff --git a/examples/kelp.jl b/examples/kelp.jl index e8e3eea2b..cb91ff6d6 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -68,6 +68,7 @@ biogeochemistry = LOBSTER(; grid, surface_phytosynthetically_active_radiation = PAR⁰, carbonates = true, variable_redfield = true, + scale_negatives = true, particles) model = NonhydrostaticModel(; grid, @@ -103,8 +104,6 @@ simulation.output_writers[:particles] = JLD2OutputWriter(model, (; particles), schedule = TimeInterval(1day), overwrite_existing = true) -scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON)) -simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) nothing #hide # ## Run! @@ -123,6 +122,7 @@ bPOC = FieldTimeSeries("$filename.jld2", "bPOC") x, y, z = nodes(P) times = P.times +nothing #hide # We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and # the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth @@ -130,10 +130,12 @@ times = P.times air_sea_CO₂_flux = zeros(length(times)) carbon_export = zeros(length(times)) +using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity + for (i, t) in enumerate(times) air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) - carbon_export[i] = sPOC[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.sPOM.w[1, 1, grid.Nz-20] + - bPOC[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.bPOM.w[1, 1, grid.Nz-20] + carbon_export[i] = sPOC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOC)).w[1, 1, grid.Nz-20] + + bPOC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOC)).w[1, 1, grid.Nz-20] end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index e295ab63c..d1fcc9767 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -3,20 +3,23 @@ module Sediments export SimpleMultiG, InstantRemineralisation using KernelAbstractions -using OceanBioME: ContinuousFormBiogeochemistry, BoxModelGrid + +using OceanBioME: Biogeochemistry, BoxModelGrid + using Oceananigans using Oceananigans.Architectures: device using Oceananigans.Utils: launch! using Oceananigans.Advection: advective_tracer_flux_z using Oceananigans.Units: day -using Oceananigans.Fields: CenterField, Face +using Oceananigans.Fields: ConstantField using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity using Oceananigans.Grids: zspacing using Oceananigans.Operators: volume -using Oceananigans.Fields: Center, ConstantField using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, immersed_cell + import Adapt: adapt_structure, adapt +import Oceananigans.Biogeochemistry: update_tendencies! abstract type AbstractSediment end abstract type FlatSediment <: AbstractSediment end @@ -25,9 +28,7 @@ abstract type FlatSediment <: AbstractSediment end sediment_fields(::AbstractSediment) = () -@inline update_sediment_tendencies!(bgc, sediment, model) = nothing - -function update_sediment_tendencies!(bgc, sediment::FlatSediment, model) +function update_tendencies!(bgc, sediment::FlatSediment, model) arch = model.grid.architecture for (i, tracer) in enumerate(sediment_tracers(sediment)) @@ -36,7 +37,7 @@ function update_sediment_tendencies!(bgc, sediment::FlatSediment, model) launch!(arch, model.grid, :xy, _calculate_tendencies!, - bgc.sediment_model, bgc, model.grid, model.advection, model.tracers, model.timestepper) + sediment, bgc, model.grid, model.advection, model.tracers, model.timestepper) return nothing end @@ -50,6 +51,10 @@ end @inline carbon_flux() = 0 @inline remineralisation_receiver() = nothing +@inline nitrogen_flux(i, j, k, grid, advection, bgc::Biogeochemistry, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc.underlying_biogeochemistry, tracers) +@inline carbon_flux(i, j, k, grid, advection, bgc::Biogeochemistry, tracers) = carbon_flux(i, j, k, grid, advection, bgc.underlying_biogeochemistry, tracers) +@inline remineralisation_receiver(bgc::Biogeochemistry) = remineralisation_receiver(bgc.underlying_biogeochemistry) + @inline advection_scheme(advection, val_tracer) = advection @inline advection_scheme(advection::NamedTuple, val_tracer::Val{T}) where T = advection[T] diff --git a/src/Boundaries/Sediments/coupled_timesteppers.jl b/src/Boundaries/Sediments/coupled_timesteppers.jl index 3942a7031..61c2f085d 100644 --- a/src/Boundaries/Sediments/coupled_timesteppers.jl +++ b/src/Boundaries/Sediments/coupled_timesteppers.jl @@ -1,5 +1,4 @@ using Oceananigans: NonhydrostaticModel, prognostic_fields, HydrostaticFreeSurfaceModel -using OceanBioME: ContinuousFormBiogeochemistry using OceanBioME.Boundaries.Sediments: AbstractSediment using Oceananigans.TimeSteppers: ab2_step_field!, rk3_substep_field!, stage_Δt using Oceananigans.Utils: work_layout, launch! @@ -9,7 +8,7 @@ using Oceananigans.Architectures: AbstractArchitecture import Oceananigans.TimeSteppers: ab2_step!, rk3_substep! -@inline function ab2_step!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:ContinuousFormBiogeochemistry{<:Any, <:FlatSediment}}, Δt, χ) +@inline function ab2_step!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Biogeochemistry{<:Any, <:Any, <:FlatSediment}}, Δt, χ) workgroup, worksize = work_layout(model.grid, :xyz) arch = model.architecture step_field_kernel! = ab2_step_field!(device(arch), workgroup, worksize) @@ -33,7 +32,7 @@ import Oceananigans.TimeSteppers: ab2_step!, rk3_substep! Δt) end - sediment = model.biogeochemistry.sediment_model + sediment = model.biogeochemistry.sediment for (i, field) in enumerate(sediment_fields(sediment)) launch!(arch, model.grid, :xy, ab2_step_flat_field!, @@ -46,14 +45,14 @@ import Oceananigans.TimeSteppers: ab2_step!, rk3_substep! return nothing end -@inline function ab2_step!(model::HydrostaticFreeSurfaceModel{<:Any, <:Any, <:AbstractArchitecture, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:ContinuousFormBiogeochemistry{<:Any, <:FlatSediment}}, Δt, χ) +@inline function ab2_step!(model::HydrostaticFreeSurfaceModel{<:Any, <:Any, <:AbstractArchitecture, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Biogeochemistry{<:Any, <:Any, <:FlatSediment}}, Δt, χ) # Step locally velocity and tracers @apply_regionally local_ab2_step!(model, Δt, χ) # blocking step for implicit free surface, non blocking for explicit ab2_step_free_surface!(model.free_surface, model, Δt, χ) - sediment = model.biogeochemistry.sediment_model + sediment = model.biogeochemistry.sediment arch = model.architecture for (i, field) in enumerate(sediment_fields(sediment)) @@ -75,7 +74,7 @@ end @inbounds u[i, j, 1] += Δt * ((one_point_five + χ) * Gⁿ[i, j, 1] - (oh_point_five + χ) * G⁻[i, j, 1]) end -function rk3_substep!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:ContinuousFormBiogeochemistry{<:Any, <:FlatSediment}}, Δt, γⁿ, ζⁿ) +function rk3_substep!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Biogeochemistry{<:Any, <:Any, <:FlatSediment}}, Δt, γⁿ, ζⁿ) workgroup, worksize = work_layout(model.grid, :xyz) arch = model.architecture substep_field_kernel! = rk3_substep_field!(device(arch), workgroup, worksize) @@ -98,7 +97,7 @@ function rk3_substep!(model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:A stage_Δt(Δt, γⁿ, ζⁿ)) end - sediment = model.biogeochemistry.sediment_model + sediment = model.biogeochemistry.sediment for (i, field) in enumerate(sediment_fields(sediment)) launch!(arch, model.grid, :xy, rk3_step_flat_field!, diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index 1ea6e8d09..62046ccba 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -101,3 +101,6 @@ sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_st @inbounds timestepper.Gⁿ[remineralisation_receiver(bgc)][i, j, 1] += flux * (1 - burial_efficiency) / Δz end + +summary(::InstantRemineralisation{FT}) where {FT} = string("Single-layer instant remineralisaiton ($FT)") +show(io::IO, model::InstantRemineralisation) = print(io, summary(model)) \ No newline at end of file diff --git a/src/Boundaries/gasexchange.jl b/src/Boundaries/gasexchange.jl index 8f18ab88b..badb7610a 100644 --- a/src/Boundaries/gasexchange.jl +++ b/src/Boundaries/gasexchange.jl @@ -236,7 +236,7 @@ function GasExchange(; gas, return FluxBoundaryCondition(gasexchange_function, field_dependencies = field_dependencies, parameters = gasexchange) end -# Hack this nicer format into Oceananigans boundary conditions because `typeof(gasexhcange) = DataType` and `BoundaryCondition` has a special use for `DataType` in the first argument so doesn't work +# Hack this nicer format into Oceananigans boundary conditions because `typeof(gasexchange) = DataType` and `BoundaryCondition` has a special use for `DataType` in the first argument so doesn't work @inline gasexchange_function(x, y, t, args...) = @inbounds args[end](x, y, t, args[1:end-1]...) @inline (gasexchange::GasExchange)(x, y, t, conc) = gasexchange(x, y, t, conc, gasexchange.temperature(x, y, 0.0, t), gasexchange.salinity(x, y, 0.0, t)) diff --git a/src/Light/2band.jl b/src/Light/2band.jl index 6f7b4bd32..35ade4e06 100644 --- a/src/Light/2band.jl +++ b/src/Light/2band.jl @@ -120,15 +120,13 @@ function TwoBandPhotosyntheticallyActiveRadiation(; grid, surface_PAR) end -function update_PAR!(model, PAR::TwoBandPhotosyntheticallyActiveRadiation) +function update_biogeochemical_state!(model, PAR::TwoBandPhotosyntheticallyActiveRadiation) arch = architecture(model.grid) launch!(arch, model.grid, :xy, update_TwoBandPhotosyntheticallyActiveRadiation!, PAR.field, model.grid, model.tracers.P, PAR.surface_PAR, model.clock.time, PAR) fill_halo_regions!(PAR.field, model.clock, fields(model)) end -required_PAR_fields(::TwoBandPhotosyntheticallyActiveRadiation) = (:PAR, ) - summary(::TwoBandPhotosyntheticallyActiveRadiation{FT}) where {FT} = string("Two-band light attenuation model ($FT)") show(io::IO, model::TwoBandPhotosyntheticallyActiveRadiation{FT}) where {FT} = print(io, summary(model)) diff --git a/src/Light/Light.jl b/src/Light/Light.jl index f34cfffad..2b36457da 100644 --- a/src/Light/Light.jl +++ b/src/Light/Light.jl @@ -23,13 +23,9 @@ using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid import Adapt: adapt_structure, adapt import Base: show, summary -import Oceananigans.Biogeochemistry: biogeochemical_auxiliary_fields +import Oceananigans.Biogeochemistry: biogeochemical_auxiliary_fields, update_biogeochemical_state!, required_biogeochemical_auxiliary_fields import Oceananigans.BoundaryConditions: _fill_top_halo! -# Fallback -update_PAR!(model, PAR) = nothing -required_PAR_fields(PAR) = () - include("2band.jl") include("morel.jl") diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index d728ade6a..69fdbba6c 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -43,23 +43,22 @@ module LOBSTERModel export LOBSTER -using OceanBioME: ContinuousFormBiogeochemistry - using Oceananigans.Units using Oceananigans.Fields: Field, TracerFields, CenterField, ZeroField -using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation, update_PAR!, required_PAR_fields -using OceanBioME: setup_velocity_fields, show_sinking_velocities +using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation +using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry, ScaleNegativeTracers using OceanBioME.BoxModels: BoxModel using OceanBioME.Boundaries.Sediments: sinking_flux +using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry + +import OceanBioME: redfield, conserved_tracers import OceanBioME.BoxModels: update_boxmodel_state! import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, required_biogeochemical_auxiliary_fields, - biogeochemical_drift_velocity, - update_biogeochemical_state!, - biogeochemical_auxiliary_fields + biogeochemical_drift_velocity import OceanBioME: maximum_sinking_velocity @@ -68,7 +67,7 @@ import Base: show, summary import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver -struct LOBSTER{FT, LA, S, B, W, P} <: ContinuousFormBiogeochemistry{LA, S, P} +struct LOBSTER{FT, B, W} <: AbstractContinuousFormBiogeochemistry phytoplankton_preference :: FT maximum_grazing_rate :: FT grazing_half_saturation :: FT @@ -99,15 +98,10 @@ struct LOBSTER{FT, LA, S, B, W, P} <: ContinuousFormBiogeochemistry{LA, S, P} disolved_organic_breakdown_rate :: FT zooplankton_calcite_dissolution :: FT - light_attenuation_model :: LA - sediment_model :: S - optionals :: B sinking_velocities :: W - particles :: P - function LOBSTER(phytoplankton_preference::FT, maximum_grazing_rate::FT, grazing_half_saturation::FT, @@ -138,53 +132,43 @@ struct LOBSTER{FT, LA, S, B, W, P} <: ContinuousFormBiogeochemistry{LA, S, P} disolved_organic_breakdown_rate::FT, zooplankton_calcite_dissolution::FT, - light_attenuation_model::LA, - sediment_model::S, - optionals::B, - sinking_velocities::W, - - particles::P) where {FT, LA, S, B, W, P} - - return new{FT, LA, S, B, W, P}(phytoplankton_preference, - maximum_grazing_rate, - grazing_half_saturation, - light_half_saturation, - nitrate_ammonia_inhibition, - nitrate_half_saturation, - ammonia_half_saturation, - maximum_phytoplankton_growthrate, - zooplankton_assimilation_fraction, - zooplankton_mortality, - zooplankton_excretion_rate, - phytoplankton_mortality, - small_detritus_remineralisation_rate, - large_detritus_remineralisation_rate, - phytoplankton_exudation_fraction, - nitrifcaiton_rate, - ammonia_fraction_of_exudate, - ammonia_fraction_of_excriment, - ammonia_fraction_of_detritus, - phytoplankton_redfield, - organic_redfield, - phytoplankton_chlorophyll_ratio, - organic_carbon_calcate_ratio, - respiraiton_oxygen_nitrogen_ratio, - nitrifcation_oxygen_nitrogen_ratio, - slow_sinking_mortality_fraction, - fast_sinking_mortality_fraction, - disolved_organic_breakdown_rate, - zooplankton_calcite_dissolution, - - light_attenuation_model, - sediment_model, - - optionals, - - sinking_velocities, - - particles) + sinking_velocities::W) where {FT, B, W} + + return new{FT, B, W}(phytoplankton_preference, + maximum_grazing_rate, + grazing_half_saturation, + light_half_saturation, + nitrate_ammonia_inhibition, + nitrate_half_saturation, + ammonia_half_saturation, + maximum_phytoplankton_growthrate, + zooplankton_assimilation_fraction, + zooplankton_mortality, + zooplankton_excretion_rate, + phytoplankton_mortality, + small_detritus_remineralisation_rate, + large_detritus_remineralisation_rate, + phytoplankton_exudation_fraction, + nitrifcaiton_rate, + ammonia_fraction_of_exudate, + ammonia_fraction_of_excriment, + ammonia_fraction_of_detritus, + phytoplankton_redfield, + organic_redfield, + phytoplankton_chlorophyll_ratio, + organic_carbon_calcate_ratio, + respiraiton_oxygen_nitrogen_ratio, + nitrifcation_oxygen_nitrogen_ratio, + slow_sinking_mortality_fraction, + fast_sinking_mortality_fraction, + disolved_organic_breakdown_rate, + zooplankton_calcite_dissolution, + + optionals, + + sinking_velocities) end end @@ -263,15 +247,16 @@ julia> grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200)); julia> model = LOBSTER(; grid) Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model (Float64) - Light Attenuation Model: - └── Two-band light attenuation model (Float64) Optional components: ├── Carbonates ❌ ├── Oxygen ❌ └── Variable Redfield Ratio ❌ Sinking Velocities: ├── sPOM: 0.0 to -3.47e-5 m/s - └── bPOM: 0.0 to -0.0023148148148148147 m/s + └── bPOM: 0.0 to -0.0023148148148148147 m/s + Light attenuation: Two-band light attenuation model (Float64) + Sediment: Nothing + Particles: Nothing ``` """ function LOBSTER(; grid, @@ -319,7 +304,10 @@ function LOBSTER(; grid, sinking_speeds = (sPOM = 3.47e-5, bPOM = 200/day), open_bottom::Bool = true, - particles::P = nothing) where {FT, LA, S, P} + scale_negatives = false, + + particles::P = nothing, + modifiers::M = nothing) where {FT, LA, S, P, M} if !isnothing(sediment_model) && !open_bottom @warn "You have specified a sediment model but not `open_bottom` which will not work as the tracer will settle in the bottom cell" @@ -329,57 +317,63 @@ function LOBSTER(; grid, optionals = Val((carbonates, oxygen, variable_redfield)) - return LOBSTER(phytoplankton_preference, - maximum_grazing_rate, - grazing_half_saturation, - light_half_saturation, - nitrate_ammonia_inhibition, - nitrate_half_saturation, - ammonia_half_saturation, - maximum_phytoplankton_growthrate, - zooplankton_assimilation_fraction, - zooplankton_mortality, - zooplankton_excretion_rate, - phytoplankton_mortality, - small_detritus_remineralisation_rate, - large_detritus_remineralisation_rate, - phytoplankton_exudation_fraction, - nitrifcaiton_rate, - ammonia_fraction_of_exudate, - ammonia_fraction_of_excriment, - ammonia_fraction_of_detritus, - phytoplankton_redfield, - organic_redfield, - phytoplankton_chlorophyll_ratio, - organic_carbon_calcate_ratio, - respiraiton_oxygen_nitrogen_ratio, - nitrifcation_oxygen_nitrogen_ratio, - slow_sinking_mortality_fraction, - fast_sinking_mortality_fraction, - disolved_organic_breakdown_rate, - zooplankton_calcite_dissolution, - - light_attenuation_model, - sediment_model, - - optionals, - - sinking_velocities, - - particles) + underlying_biogeochemistry = LOBSTER(phytoplankton_preference, + maximum_grazing_rate, + grazing_half_saturation, + light_half_saturation, + nitrate_ammonia_inhibition, + nitrate_half_saturation, + ammonia_half_saturation, + maximum_phytoplankton_growthrate, + zooplankton_assimilation_fraction, + zooplankton_mortality, + zooplankton_excretion_rate, + phytoplankton_mortality, + small_detritus_remineralisation_rate, + large_detritus_remineralisation_rate, + phytoplankton_exudation_fraction, + nitrifcaiton_rate, + ammonia_fraction_of_exudate, + ammonia_fraction_of_excriment, + ammonia_fraction_of_detritus, + phytoplankton_redfield, + organic_redfield, + phytoplankton_chlorophyll_ratio, + organic_carbon_calcate_ratio, + respiraiton_oxygen_nitrogen_ratio, + nitrifcation_oxygen_nitrogen_ratio, + slow_sinking_mortality_fraction, + fast_sinking_mortality_fraction, + disolved_organic_breakdown_rate, + zooplankton_calcite_dissolution, + + optionals, + + sinking_velocities) + + if scale_negatives + scaler = ScaleNegativeTracers(underlying_biogeochemistry) + modifiers = isnothing(modifiers) ? scaler : (modifiers..., scaler) + end + + return Biogeochemistry(underlying_biogeochemistry; + light_attenuation = light_attenuation_model, + sediment = sediment_model, + particles, + modifiers) end # wrote this functionally and it took 2.5x longer so even though this is long going to use this way instead -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, false, false)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, false)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :DIC, :Alk) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, false)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :O₂) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, false, true)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :sPOC, :bPOC, :DOC) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, false)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :DIC, :Alk, :O₂) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :DIC, :Alk, :sPOC, :bPOC, :DOC) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :O₂, :sPOC, :bPOC, :DOC) -@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, true)}, <:Any, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :DIC, :Alk, :O₂, :sPOC, :bPOC, :DOC) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(false, false, false)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(true, false, false)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :DIC, :Alk) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(false, true, false)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :O₂) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(false, false, true)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :sPOC, :bPOC, :DOC) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(true, true, false)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM, :DIC, :Alk, :O₂) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :DIC, :Alk, :sPOC, :bPOC, :DOC) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :O₂, :sPOC, :bPOC, :DOC) +@inline required_biogeochemical_tracers(::LOBSTER{<:Any, <:Val{(true, true, true)}, <:Any}) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON, :DIC, :Alk, :O₂, :sPOC, :bPOC, :DOC) -@inline required_biogeochemical_auxiliary_fields(model::LOBSTER) = required_PAR_fields(model.light_attenuation_model) +required_biogeochemical_auxiliary_fields(::LOBSTER) = (:PAR, ) const small_detritus = Union{Val{:sPON}, Val{:sPOC}} const large_detritus = Union{Val{:bPON}, Val{:bPOC}} @@ -401,11 +395,7 @@ const DOM = Union{Val{:DOM}, Val{:DON}} end end -function update_biogeochemical_state!(bgc::LOBSTER, model) - update_PAR!(model, bgc.light_attenuation_model) -end - -function update_boxmodel_state!(model::BoxModel{<:LOBSTER, <:Any, <:Any, <:Any, <:Any, <:Any}) +function update_boxmodel_state!(model::BoxModel{<:Biogeochemistry{<:LOBSTER}, <:Any, <:Any, <:Any, <:Any, <:Any}) getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time) end @@ -439,26 +429,19 @@ adapt_structure(to, lobster::LOBSTER) = lobster.fast_sinking_mortality_fraction, lobster.disolved_organic_breakdown_rate, lobster.zooplankton_calcite_dissolution, - adapt(to, lobster.light_attenuation_model), - adapt(to, lobster.sediment_model), lobster.optionals, - adapt(to, lobster.sinking_velocities), - adapt(to, lobster.particles)) - -summary(::LOBSTER{FT, LA, S, B, W, P}) where {FT, LA, S, B, W, P} = string("Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model ($FT)") -show(io::IO, model::LOBSTER{FT, LA, S, Val{B}, W, P}) where {FT, LA, S, B, W, P} = - print(io, summary(model), " \n", - " Light Attenuation Model: ", "\n", - " └── ", summary(model.light_attenuation_model), "\n", - " Optional components:", "\n", - " ├── Carbonates $(B[1] ? :✅ : :❌) \n", - " ├── Oxygen $(B[2] ? :✅ : :❌) \n", - " └── Variable Redfield Ratio $(B[3] ? :✅ : :❌)", "\n", - " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) + adapt(to, lobster.sinking_velocities)) -@inline maximum_sinking_velocity(bgc::LOBSTER) = maximum(abs, bgc.sinking_velocities.bPOM.w) +summary(::LOBSTER{FT, B, W}) where {FT, B, W} = string("Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model ($FT)") -@inline biogeochemical_auxiliary_fields(bgc::LOBSTER) = biogeochemical_auxiliary_fields(bgc.light_attenuation_model) +show(io::IO, model::LOBSTER{FT, Val{B}, W}) where {FT, B, W} = string(summary(model), " \n", + " Optional components:", "\n", + " ├── Carbonates $(B[1] ? :✅ : :❌) \n", + " ├── Oxygen $(B[2] ? :✅ : :❌) \n", + " └── Variable Redfield Ratio $(B[3] ? :✅ : :❌)", "\n", + " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) + +@inline maximum_sinking_velocity(bgc::LOBSTER) = maximum(abs, bgc.sinking_velocities.bPOM.w) include("fallbacks.jl") @@ -467,17 +450,32 @@ include("carbonate_chemistry.jl") include("oxygen_chemistry.jl") include("variable_redfield.jl") -const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, false, true)}, <:Any, <:Any}, - LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any}, - LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any}, - LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, true)}, <:Any, <:Any}} +const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Val{(false, false, true)}, <:Any}, + LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any}, + LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any}, + LOBSTER{<:Any, <:Val{(true, true, true)}, <:Any}} +@inline redfield(i, j, k, val_tracer_name, bgc::LOBSTER, tracers) = redfield(val_tracer_name, bgc) + +@inline redfield(::Val{:P}, bgc::LOBSTER) = (1 + bgc.organic_carbon_calcate_ratio) * bgc.phytoplankton_redfield +@inline redfield(::Val{:Z}, bgc::LOBSTER) = bgc.phytoplankton_redfield +@inline redfield(::Union{Val{:NO₃}, Val{:NH₄}, Val{:Alk}, Val{:O₂}}, bgc::LOBSTER) = 0 +@inline redfield(::Union{Val{:sPOM}, Val{:bPOM}, Val{:DOM}}, bgc::LOBSTER) = bgc.organic_redfield +@inline redfield(::Union{Val{:sPOC}, Val{:bPOC}, Val{:DOC}, Val{:DIC}}, bgc::LOBSTER) = 1 + +@inline redfield(i, j, k, ::Val{:sPON}, bgc::lobster_variable_redfield, tracers) = @inbounds tracers.sPOC[i, j, k] / tracers.sPON[i, j, k] +@inline redfield(i, j, k, ::Val{:bPON}, bgc::lobster_variable_redfield, tracers) = @inbounds tracers.bPOC[i, j, k] / tracers.bPON[i, j, k] +@inline redfield(i, j, k, ::Val{:DON}, bgc::lobster_variable_redfield, tracers) = @inbounds tracers.DOC[i, j, k] / tracers.DON[i, j, k] + +@inline redfield(::Val{:sPON}, bgc::lobster_variable_redfield, tracers) = tracers.sPOC / tracers.sPON +@inline redfield(::Val{:bPON}, bgc::lobster_variable_redfield, tracers) = tracers.bPOC / tracers.bPON +@inline redfield(::Val{:DON}, bgc::lobster_variable_redfield, tracers) = tracers.DOC / tracers.DON @inline nitrogen_flux(i, j, k, grid, advection, bgc::LOBSTER, tracers) = sinking_flux(i, j, k, grid, advection, Val(:sPOM), bgc, tracers) + sinking_flux(i, j, k, grid, advection, Val(:bPOM), bgc, tracers) -@inline carbon_flux(i, j, k, grid, advection, bgc::LOBSTER, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * bgc.organic_redfield +@inline carbon_flux(i, j, k, grid, advection, bgc::LOBSTER, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * redfield(Val(:sPOM), bgc) @inline nitrogen_flux(i, j, k, grid, advection, bgc::lobster_variable_redfield, tracers) = sinking_flux(i, j, k, grid, advection, Val(:sPON), bgc, tracers) + @@ -488,4 +486,7 @@ const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Any, <:Any, <:Val{(fals sinking_flux(i, j, k, grid, advection, Val(:bPOC), bgc, tracers) @inline remineralisation_receiver(::LOBSTER) = :NH₄ + +@inline conserved_tracers(::LOBSTER) = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM) +@inline conserved_tracers(::lobster_variable_redfield) = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON) end # module diff --git a/src/Models/AdvectedPopulations/LOBSTER/fallbacks.jl b/src/Models/AdvectedPopulations/LOBSTER/fallbacks.jl index bf9cc5d28..fd8902c26 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/fallbacks.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/fallbacks.jl @@ -12,22 +12,22 @@ # Carbonates and oxygen # We can't tell the difference between carbonates and oxygen, and variable redfields without specifying more # We're lucky that the optional groups have 2, 1, 3 extra variables each so this is the only non-unique conbination -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, false)}, <:Any, <:Any})(tracer::Val{:DIC}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, false)}, <:Any, <:Any})(tracer::Val{:Alk}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, false)}, <:Any, <:Any})(tracer::Val{:O₂}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, O₂, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, true, false)}, <:Any})(tracer::Val{:DIC}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, true, false)}, <:Any})(tracer::Val{:Alk}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, true, false)}, <:Any})(tracer::Val{:O₂}, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk, O₂, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, O₂, PAR) # Carbonates and variable redfield -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any})(tracer::Val{:DIC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOC / bgc. organic_redfield, bPOC / bgc. organic_redfield, DOC / bgc. organic_redfield, DIC, Alk, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any})(tracer::Val{:Alk}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOC / bgc. organic_redfield, bPOC / bgc. organic_redfield, DOC / bgc. organic_redfield, DIC, Alk, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any})(tracer::Val{:sPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any})(tracer::Val{:bPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, false, true)}, <:Any, <:Any})(tracer::Val{:DOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(Val(:DOM), x, y, z, t, NO₃, NH₄, P * bgc.phytoplankton_redfield, Z * bgc.phytoplankton_redfield, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any})(tracer::Val{:DIC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOC / bgc. organic_redfield, bPOC / bgc. organic_redfield, DOC / bgc. organic_redfield, DIC, Alk, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any})(tracer::Val{:Alk}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOC / bgc. organic_redfield, bPOC / bgc. organic_redfield, DOC / bgc. organic_redfield, DIC, Alk, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any})(tracer::Val{:sPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any})(tracer::Val{:bPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(true, false, true)}, <:Any})(::Val{:DOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, sPOC, bPOC, DOC, PAR) = bgc(Val(:DOM), x, y, z, t, NO₃, NH₄, P * bgc.phytoplankton_redfield, Z * bgc.phytoplankton_redfield, sPOC, bPOC, DOC, PAR) # Oxygen and variable redfield -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any})(tracer::Val{:O₂}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any})(tracer::Val{:sPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any})(tracer::Val{:bPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) -@inline (bgc::LOBSTER{<:Any, <:Any, <:Any, <:Val{(false, true, true)}, <:Any, <:Any})(tracer::Val{:DOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(Val(:DOM), x, y, z, t, NO₃, NH₄, P * bgc.phytoplankton_redfield, Z * bgc.phytoplankton_redfield, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any})(tracer::Val{:O₂}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any})(tracer::Val{:sPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any})(tracer::Val{:bPOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, sPOC, bPOC, DOC, PAR) +@inline (bgc::LOBSTER{<:Any, <:Val{(false, true, true)}, <:Any})(::Val{:DOC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, O₂, sPOC, bPOC, DOC, PAR) = bgc(Val(:DOM), x, y, z, t, NO₃, NH₄, P * bgc.phytoplankton_redfield, Z * bgc.phytoplankton_redfield, sPOC, bPOC, DOC, PAR) # Carbonates, oxygen, and variable redfield @inline (bgc::LOBSTER)(tracer::Val{:DIC}, x, y, z, t, NO₃, NH₄, P, Z, sPON, bPON, DON, DIC, Alk, O₂, sPOC, bPOC, DOC, PAR) = bgc(tracer, x, y, z, t, NO₃, NH₄, P, Z, sPOC / bgc. organic_redfield, bPOC / bgc. organic_redfield, DOC / bgc. organic_redfield, DIC, Alk, PAR) diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 704bbb0e6..ed892a72e 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -16,24 +16,24 @@ module NPZDModel export NutrientPhytoplanktonZooplanktonDetritus, NPZD -using OceanBioME: ContinuousFormBiogeochemistry +using OceanBioME: Biogeochemistry, ScaleNegativeTracers +using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry using Oceananigans.Units using Oceananigans.Fields: ZeroField -using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation, update_PAR!, required_PAR_fields +using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation using OceanBioME: setup_velocity_fields, show_sinking_velocities using OceanBioME.BoxModels: BoxModel using OceanBioME.Boundaries.Sediments: sinking_flux +import OceanBioME: redfield, conserved_tracers import OceanBioME.BoxModels: update_boxmodel_state! import Base: show, summary import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, required_biogeochemical_auxiliary_fields, - update_biogeochemical_state!, - biogeochemical_drift_velocity, - biogeochemical_auxiliary_fields + biogeochemical_drift_velocity import OceanBioME: maximum_sinking_velocity @@ -41,7 +41,7 @@ import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisa import Adapt: adapt_structure, adapt -struct NutrientPhytoplanktonZooplanktonDetritus{FT, LA, S, W, P} <: ContinuousFormBiogeochemistry{LA, S, P} +struct NutrientPhytoplanktonZooplanktonDetritus{FT, W} <: AbstractContinuousFormBiogeochemistry # phytoplankton initial_photosynthetic_slope :: FT # α, 1/(W/m²)/s base_maximum_growth :: FT # μ₀, 1/s @@ -59,17 +59,9 @@ struct NutrientPhytoplanktonZooplanktonDetritus{FT, LA, S, W, P} <: ContinuousFo # detritus remineralization_rate :: FT # rᵈⁿ, 1/s - # light attenuation - light_attenuation_model :: LA - - # sediment - sediment_model :: S - # sinking sinking_velocities :: W - particles :: P - function NutrientPhytoplanktonZooplanktonDetritus(initial_photosynthetic_slope::FT, base_maximum_growth::FT, nutrient_half_saturation::FT, @@ -84,34 +76,22 @@ struct NutrientPhytoplanktonZooplanktonDetritus{FT, LA, S, W, P} <: ContinuousFo remineralization_rate::FT, - light_attenuation_model::LA, - - sediment_model::S, - - sinking_velocities::W, - - particles::P) where {FT, LA, S, W, P} - return new{FT, LA, S, W, P}(initial_photosynthetic_slope, - base_maximum_growth, - nutrient_half_saturation, - base_respiration_rate, - phyto_base_mortality_rate, - - maximum_grazing_rate, - grazing_half_saturation, - assimulation_efficiency, - base_excretion_rate, - zoo_base_mortality_rate, - - remineralization_rate, - - light_attenuation_model, - - sediment_model, - - sinking_velocities, - - particles) + sinking_velocities::W) where {FT, W} + return new{FT, W}(initial_photosynthetic_slope, + base_maximum_growth, + nutrient_half_saturation, + base_respiration_rate, + phyto_base_mortality_rate, + + maximum_grazing_rate, + grazing_half_saturation, + assimulation_efficiency, + base_excretion_rate, + zoo_base_mortality_rate, + + remineralization_rate, + + sinking_velocities) end end @@ -166,11 +146,13 @@ julia> grid = RectilinearGrid(size=(20, 30), extent=(200, 200), topology=(Bounde julia> model = NutrientPhytoplanktonZooplanktonDetritus(; grid) Nutrient Phytoplankton Zooplankton Detritus model (Float64) - Light Attenuation Model: - └── Two-band light attenuation model (Float64) Sinking Velocities: ├── P: 0.0 to -2.9525462962962963e-6 m/s - └── D: 0.0 to -3.181597222222222e-5 m/s + └── D: 0.0 to -3.181597222222222e-5 m/s + Light attenuation: Two-band light attenuation model (Float64) + Sediment: Nothing + Particles: Nothing + ``` """ function NutrientPhytoplanktonZooplanktonDetritus(; grid, @@ -194,30 +176,38 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid, sinking_speeds = (P = 0.2551/day, D = 2.7489/day), open_bottom::Bool = true, + + scale_negatives = false, - particles::P = nothing) where {FT, LA, S, P} + particles::P = nothing, + modifiers::M = nothing) where {FT, LA, S, P, M} sinking_velocities = setup_velocity_fields(sinking_speeds, grid, open_bottom) - if !isnothing(sediment_model) && !open_bottom - @warn "You have specified a sediment model but not `open_bottom` which will not work as the tracer will settle in the bottom cell" + underlying_biogeochemistry = + NutrientPhytoplanktonZooplanktonDetritus(initial_photosynthetic_slope, + base_maximum_growth, + nutrient_half_saturation, + base_respiration_rate, + phyto_base_mortality_rate, + maximum_grazing_rate, + grazing_half_saturation, + assimulation_efficiency, + base_excretion_rate, + zoo_base_mortality_rate, + remineralization_rate, + sinking_velocities) + + if scale_negatives + scaler = ScaleNegativeTracers(underlying_biogeochemistry) + modifiers = isnothing(modifiers) ? scaler : (modifiers..., scaler) end - return NutrientPhytoplanktonZooplanktonDetritus(initial_photosynthetic_slope, - base_maximum_growth, - nutrient_half_saturation, - base_respiration_rate, - phyto_base_mortality_rate, - maximum_grazing_rate, - grazing_half_saturation, - assimulation_efficiency, - base_excretion_rate, - zoo_base_mortality_rate, - remineralization_rate, - light_attenuation_model, - sediment_model, - sinking_velocities, - particles) + return Biogeochemistry(underlying_biogeochemistry; + light_attenuation = light_attenuation_model, + sediment = sediment_model, + particles, + modifiers) end const NPZD = NutrientPhytoplanktonZooplanktonDetritus @@ -302,24 +292,16 @@ end end end -function update_biogeochemical_state!(bgc::NPZD, model) - update_PAR!(model, bgc.light_attenuation_model) -end - -function update_boxmodel_state!(model::BoxModel{<:NPZD, <:Any, <:Any, <:Any, <:Any, <:Any}) +function update_boxmodel_state!(model::BoxModel{<:Biogeochemistry{<:NPZD}, <:Any, <:Any, <:Any, <:Any, <:Any}) getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time) getproperty(model.values, :T) .= model.forcing.T(model.clock.time) end -summary(::NPZD{FT, LA, W, P}) where {FT, LA, W, P} = string("Nutrient Phytoplankton Zooplankton Detritus model ($FT)") -show(io::IO, model::NPZD{FT, LA, W, P}) where {FT, LA, W, P} = - print(io, summary(model), " \n", - " Light Attenuation Model: ", "\n", - " └── ", summary(model.light_attenuation_model), "\n", - " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) +summary(::NPZD{FT, W}) where {FT, W} = string("Nutrient Phytoplankton Zooplankton Detritus model ($FT)") +show(io::IO, model::NPZD) = string(summary(model), " \n", + " Sinking Velocities:", "\n", show_sinking_velocities(model.sinking_velocities)) @inline maximum_sinking_velocity(bgc::NPZD) = maximum(abs, bgc.sinking_velocities.D.w) -@inline biogeochemical_auxiliary_fields(bgc::NPZD) = biogeochemical_auxiliary_fields(bgc.light_attenuation_model) adapt_structure(to, npzd::NPZD) = NutrientPhytoplanktonZooplanktonDetritus(npzd.initial_photosynthetic_slope, @@ -336,17 +318,18 @@ adapt_structure(to, npzd::NPZD) = npzd.remineralization_rate, - adapt(to, npzd.light_attenuation_model), - adapt(to, npzd.sediment_model), + adapt(to, npzd.sinking_velocities)) - adapt(to, npzd.sinking_velocities), - - adapt(to, npzd.particles)) +@inline redfield(i, j, k, val_tracer_name, bgc::NPZD, tracers) = redfield(val_tracer_name, bgc) +@inline redfield(::Union{Val{:N}}, bgc::NPZD) = 0 +@inline redfield(::Union{Val{:P}, Val{:Z}, Val{:D}}, bgc::NPZD) = 6.56 @inline nitrogen_flux(i, j, k, grid, advection, bgc::NPZD, tracers) = sinking_flux(i, j, k, grid, advection, Val(:D), bgc, tracers) + sinking_flux(i, j, k, grid, advection, Val(:P), bgc, tracers) -@inline carbon_flux(i, j, k, grid, advection, bgc::NPZD, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * 6.56 +@inline carbon_flux(i, j, k, grid, advection, bgc::NPZD, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * redfield(Val(:P), bgc) @inline remineralisation_receiver(::NPZD) = :N + +@inline conserved_tracers(::NPZD) = (:N, :P, :Z, :D) end # module diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index 5b2a68ae2..1d8ead7c3 100644 --- a/src/Models/Individuals/SLatissima.jl +++ b/src/Models/Individuals/SLatissima.jl @@ -32,7 +32,7 @@ using Oceananigans.Grids: AbstractGrid using Oceananigans.Fields: fractional_indices, _interpolate, datatuple import Adapt: adapt_structure, adapt -import OceanBioME.Particles: update_particles_tendencies! +import Oceananigans.Biogeochemistry: update_tendencies! import Oceananigans.Models.LagrangianParticleTracking: update_lagrangian_particle_properties!, _advect_particles! @inline no_dynamics(args...) = nothing @@ -264,7 +264,7 @@ adapt_structure(to, kelp::SLatissima) = SLatissima(kelp.architecture, kelp.scalefactor, kelp.latitude) -function update_particles_tendencies!(bgc, particles::SLatissima, model) +function update_tendencies!(bgc, particles::SLatissima, model) num_particles = length(particles) workgroup = min(num_particles, 256) worksize = num_particles diff --git a/src/Models/Models.jl b/src/Models/Models.jl index af590ffc1..c50245e38 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -1,14 +1,3 @@ -using OceanBioME.Boundaries.Sediments: update_sediment_tendencies! -using OceanBioME.Particles: update_particles_tendencies! - -import Oceananigans.Biogeochemistry: update_tendencies! - include("AdvectedPopulations/LOBSTER/LOBSTER.jl") include("AdvectedPopulations/NPZD.jl") -include("Individuals/SLatissima.jl") - -@inline function update_tendencies!(bgc::ContinuousFormBiogeochemistry, model) - update_sediment_tendencies!(bgc, bgc.sediment_model, model) - update_particles_tendencies!(bgc, bgc.particles, model) - return nothing -end +include("Individuals/SLatissima.jl") \ No newline at end of file diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 432dcf2a3..a7b872642 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -5,7 +5,7 @@ between ocean biogeochemistry, carbonate chemistry, and physics. module OceanBioME # Biogeochemistry models -export LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, PISCES, NPZD +export Biogeochemistry, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus, NPZD, redfield # Macroalgae models export SLatissima @@ -17,7 +17,7 @@ export BoxModel, BoxModelGrid, SaveBoxModel, run!, set! export Particles # Light models -export TwoBandPhotosyntheticallyActiveRadiation, update_PAR! +export TwoBandPhotosyntheticallyActiveRadiation # Boundaries export Boundaries, Sediments, GasExchange, FlatSediment @@ -26,19 +26,146 @@ export Boundaries, Sediments, GasExchange, FlatSediment export column_advection_timescale, column_diffusion_timescale, sinking_advection_timescale, Budget # Positivity preservaiton utilities -export zero_negative_tracers!, error_on_neg!, warn_on_neg!, ScaleNegativeTracers, remove_NaN_tendencies! +export ScaleNegativeTracers, ZeroNegativeTracers # Oceananigans extensions export ColumnField, isacolumn using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry +using Adapt + +import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, + required_biogeochemical_auxiliary_fields, + update_biogeochemical_state!, + biogeochemical_drift_velocity, + biogeochemical_auxiliary_fields, + update_tendencies!, + biogeochemical_transition + +import Adapt: adapt_structure +import Base: show, summary + +struct Biogeochemistry{B, L, S, P, M} <: AbstractContinuousFormBiogeochemistry + underlying_biogeochemistry :: B + light_attenuation :: L + sediment :: S + particles :: P + modifiers :: M + + Biogeochemistry(underlying_biogeochemistry::B, + light_attenuation::L, + sediment::S, + particles::P, + modifiers::M) where {B, L, S, P, M} = + new{B, L, S, P, M}(underlying_biogeochemistry, + light_attenuation, + sediment, + particles, + modifiers) +end -abstract type ContinuousFormBiogeochemistry{LA, S, P} <: AbstractContinuousFormBiogeochemistry end +""" + Biogeochemistry(underlying_biogeochemistry; + light_attenuation = nothing, + sediment = nothing, + particles = nothing, + modifiers = nothing) + +Construct a biogeochemical model based on `underlying_biogeochemistry` which may have +a `light_attenuation` model, `sediment`, `particles`, and `modifiers`. + +Keyword Arguments +================= + +- `light_attenuation_model`: light attenuation model which integrated the attenuation of available light +- `sediment_model`: slot for `AbstractSediment` +- `particles`: slot for `BiogeochemicalParticles` +- `modifiers`: slot for components which modfiy the biogeochemistry when the tendencies have been calculated or when the state is updated +""" +Biogeochemistry(underlying_biogeochemistry; + light_attenuation = nothing, + sediment = nothing, + particles = nothing, + modifiers = nothing) = + Biogeochemistry(underlying_biogeochemistry, + light_attenuation, + sediment, + particles, + modifiers) + +required_biogeochemical_tracers(bgc::Biogeochemistry) = required_biogeochemical_tracers(bgc.underlying_biogeochemistry) + +required_biogeochemical_auxiliary_fields(bgc::Biogeochemistry) = required_biogeochemical_auxiliary_fields(bgc.underlying_biogeochemistry) + +biogeochemical_drift_velocity(bgc::Biogeochemistry, val_name) = biogeochemical_drift_velocity(bgc.underlying_biogeochemistry, val_name) + +biogeochemical_auxiliary_fields(bgc::Biogeochemistry) = merge(biogeochemical_auxiliary_fields(bgc.underlying_biogeochemistry), + biogeochemical_auxiliary_fields(bgc.light_attenuation)) + +adapt_structure(to, bgc::Biogeochemistry) = Biogeochemistry(adapt(to, bgc.underlying_biogeochemistry), + adapt(to, bgc.light_attenuation), + adapt(to, bgc.sediment), + adapt(to, bgc.particles), + adapt(to, bgc.modifiers)) + +function update_tendencies!(bgc::Biogeochemistry, model) + update_tendencies!(bgc, bgc.sediment, model) + update_tendencies!(bgc, bgc.particles, model) + update_tendencies!(bgc, bgc.modifiers, model) +end + +update_tendencies!(bgc, modifier, model) = nothing +update_tendencies!(bgc, modifiers::Tuple, model) = [update_tendencies!(bgc, modifier, model) for modifier in modifiers] + +@inline biogeochemical_transition(i, j, k, grid, bgc::Biogeochemistry, val_tracer_name, clock, fields) = + biogeochemical_transition(i, j, k, grid, bgc.underlying_biogeochemistry, val_tracer_name, clock, fields) + +function update_biogeochemical_state!(bgc::Biogeochemistry, model) + update_biogeochemical_state!(model, bgc.modifiers) + update_biogeochemical_state!(model, bgc.light_attenuation) +end + +update_biogeochemical_state!(model, modifiers::Tuple) = [update_biogeochemical_state!(model, modifier) for modifier in modifiers] struct BoxModelGrid end @inline maximum_sinking_velocity(bgc) = 0.0 +""" + redfield(i, j, k, val_tracer_name, bgc, tracers) + +Returns the redfield ratio of `tracer_name` from `bgc` at `i`, `j`, `k`. +""" +@inline redfield(i, j, k, val_tracer_name, bgc, tracers) = redfield(val_tracer_name, bgc, tracers) + +""" + redfield(val_tracer_name, bgc) + redfield(val_tracer_name, bgc, tracers) + +Returns the redfield ratio of `tracer_name` from `bgc` when it is constant across the domain. +""" +@inline redfield(val_tracer_name, bgc) = NaN + +# fallbacks +@inline redfield(i, j, k, val_tracer_name, bgc::Biogeochemistry, tracers) = redfield(i, j, k, val_tracer_name, bgc.underlying_biogeochemistry, tracers) +@inline redfield(val_tracer_name, bgc::Biogeochemistry) = redfield(val_tracer_name, bgc.underlying_biogeochemistry) +@inline redfield(val_tracer_name, bgc::Biogeochemistry, tracers) = redfield(val_tracer_name, bgc.underlying_biogeochemistry, tracers) +@inline redfield(val_tracer_name, bgc, tracers) = redfield(val_tracer_name, bgc) + +""" + conserved_tracers(model::UnderlyingBiogeochemicalModel) + +Returns the names of tracers which together are conserved in `model` +""" +conserved_tracers(model::Biogeochemistry) = conserved_tracers(model.underlying_biogeochemistry) + +summary(bgc::Biogeochemistry) = string("Biogeochemical model based on $(summary(bgc.underlying_biogeochemistry))") +show(io::IO, model::Biogeochemistry) = + print(io, show(model.underlying_biogeochemistry), " \n", + " Light attenuation: ", summary(model.light_attenuation), "\n", + " Sediment: ", summary(model.sediment), "\n", + " Particles: ", summary(model.particles)) + include("Utils/Utils.jl") include("Boundaries/Boundaries.jl") include("Light/Light.jl") diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index 8972bdfa7..a9f28472b 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -1,8 +1,9 @@ module Particles -using OceanBioME: ContinuousFormBiogeochemistry using Oceananigans: NonhydrostaticModel, HydrostaticFreeSurfaceModel +using OceanBioME: Biogeochemistry +import Oceananigans.Biogeochemistry: update_tendencies! import Oceananigans.Models.LagrangianParticleTracking: update_lagrangian_particle_properties!, step_lagrangian_particles! import Oceananigans.OutputWriters: fetch_output import Base: length, size, show, summary @@ -13,20 +14,19 @@ abstract type BiogeochemicalParticles end @inline step_lagrangian_particles!(::Nothing, model::NonhydrostaticModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, - <:ContinuousFormBiogeochemistry{<:Any, <:Any, <:BiogeochemicalParticles}}, + <:Biogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}}, Δt) = update_lagrangian_particle_properties!(model, model.biogeochemistry, Δt) @inline step_lagrangian_particles!(::Nothing, model::HydrostaticFreeSurfaceModel{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, - <:ContinuousFormBiogeochemistry{<:Any, <:Any, <:BiogeochemicalParticles}, + <:Biogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}, <:Any, <:Any, <:Any, <:Any, <:Any,}, Δt) = update_lagrangian_particle_properties!(model, model.biogeochemistry, Δt) -@inline update_lagrangian_particle_properties!(model, bgc::ContinuousFormBiogeochemistry{<:Any, <:Any, <:BiogeochemicalParticles}, Δt) = +@inline update_lagrangian_particle_properties!(model, bgc::Biogeochemistry{<:Any, <:Any, <:Any, <:BiogeochemicalParticles}, Δt) = update_lagrangian_particle_properties!(bgc.particles, model, bgc, Δt) update_lagrangian_particle_properties!(::BiogeochemicalParticles, model, bgc, Δt) = nothing -update_particles_tendencies!(bgc, particles, model) = nothing size(particles::BiogeochemicalParticles) = size(particles.x) length(particles::BiogeochemicalParticles) = length(particles.x) diff --git a/src/Utils/negative_tracers.jl b/src/Utils/negative_tracers.jl index 79492f75f..a49229adc 100644 --- a/src/Utils/negative_tracers.jl +++ b/src/Utils/negative_tracers.jl @@ -3,64 +3,28 @@ using KernelAbstractions using KernelAbstractions.Extras.LoopInfo: @unroll using Oceananigans.Utils: work_layout using Oceananigans.Architectures: device +using Oceananigans.Biogeochemistry: AbstractBiogeochemistry import Adapt: adapt_structure, adapt +import Oceananigans.Biogeochemistry: update_tendencies!, update_biogeochemical_state! """ - zero_negative_tracers!(sim; params = (exclude=(), )) + ZeroNegativeTracers(; exclude = ()) -Sets any tracers in `sim.model` which are negative to zero. Use like: -```julia -simulation.callbacks[:neg] = Callback(zero_negative_tracers!) -``` - -Tracers to exclude can be set in the parameters. +Construct a modifier that zeroes any negative tracers excluding those listed in `exclude`. !!! danger "Tracer conservation" This method is _not_ recommended as a way to preserve positivity of tracers since it does not conserve the total tracer. """ -function zero_negative_tracers!(model; params = (exclude=(), )) - @unroll for (tracer_name, tracer) in pairs(model.tracers) - if !(tracer_name in params.exclude) - parent(tracer) .= max.(0.0, parent(tracer)) - end - end -end -@inline zero_negative_tracers!(sim::Simulation; params = (exclude=(), )) = zero_negative_tracers!(sim.model; params) - -""" - error_on_neg(sim; params = (exclude=(), )) - -Throws an error if any tracers in `sim.model` are negative. Use like: -```julia -simulation.callbacks[:neg] = Callback(error_on_neg!) -``` - -Tracers to exclude can be set in the parameters. -""" -function error_on_neg!(sim; params = (exclude=(), )) - @unroll for (tracer_name, tracer) in pairs(sim.model.tracers) - if !(tracer_name in params.exclude) - if any(tracer .< 0.0) error("$tracer_name < 0") end - end - end +@kwdef struct ZeroNegativeTracers{E} + exclude :: E = () end -""" - warn_on_neg(sim; params = (exclude=(), )) - -Raises a warning if any tracers in `sim.model` are negative. Use like: -```julia -simulation.callbacks[:neg] = Callback(warn_on_neg!) -``` - -Tracers to exclude can be set in the parameters. -""" -function warn_on_neg!(sim; params = (exclude=(), )) - @unroll for (tracer_name, tracer) in pairs(sim.model.tracers) - if !(tracer_name in params.exclude) - if any(tracer .< 0.0) @warn "$tracer_name < 0" end +function update_biogeochemical_state!(model, zero::ZeroNegativeTracers) + @unroll for (tracer_name, tracer) in pairs(model.tracers) + if !(tracer_name in zero.exclude) + parent(tracer) .= max.(0.0, parent(tracer)) end end end @@ -69,32 +33,6 @@ end ##### Infastructure to rescale negative values ##### -@kernel function scale_for_negs!(fields, tracers, scalefactors) - i, j, k = @index(Global, NTuple) - t, p = 0.0, 0.0 - @unroll for (idx, tracer) in enumerate(tracers) - field = @inbounds fields[tracer][i, j, k] - scalefactor = @inbounds scalefactors[idx] - - t += field * scalefactor - if field > 0 - p += field * scalefactor - end - end - t < 0 && error("Cell total < 0, can not scale negative tracers.") - @unroll for tracer in tracers - field = @inbounds fields[tracer][i, j, k] - - if field > 0 - field *= t / p - else - field = 0 - end - - @inbounds fields[tracer][i, j, k] = field - end -end - struct ScaleNegativeTracers{FA, SA, W} tracers :: FA scalefactors :: SA @@ -111,54 +49,70 @@ adapt_structure(to, snt::ScaleNegativeTracers) = ScaleNegativeTracers(adapt(to, """ ScaleNegativeTracers(; tracers, scalefactors = ones(length(tracers)), warn = false) -Returns a callback that scales `tracers` so that none are negative. Use like: +Constructs a modifier to scale `tracers` so that none are negative. Use like: ```julia -negativity_protection! = ScaleNegativeTracers(; model, tracers = (:P, :Z, :N)) -simulation.callbacks[:neg] = Callback(negativity_protection!; callsite = UpdateStateCallsite()) +modifier = ScaleNegativeTracers((:P, :Z, :N)) +biogeochemistry = Biogeochemistry(...; modifier) ``` This method is better, though still imperfect, method to prevent numerical errors that lead to -negative tracer values compared to [`zero_negative_tracers!`](@ref). Please see [discussion in +negative tracer values compared to [`ZeroNegativeTracers`](@ref). Please see [discussion in github](https://github.com/OceanBioME/OceanBioME.jl/discussions/48). Future plans include implement a positivity-preserving timestepping scheme as the ideal alternative. + +If `warn` is true then scaling will raise a warning. """ -function ScaleNegativeTracers(; model, tracers, scalefactors = NamedTuple{tracers}(ones(length(tracers))), warn = false) +function ScaleNegativeTracers(tracers; scalefactors = NamedTuple{tracers}(ones(length(tracers))), warn = false) if length(scalefactors) != length(tracers) error("Incorrect number of scale factors provided") end - tracers = ntuple(n -> indexin([tracers[n]], [keys(model.tracers)...])[1], length(tracers)) - scalefactors = values(scalefactors) + return ScaleNegativeTracers(tracers, scalefactors, warn) +end + +""" + ScaleNegativeTracers(model::UnderlyingBiogeochemicalModel; warn = false) + +Construct a modifier to scale the conserved tracers in `model`. + +If `warn` is true then scaling will raise a warning. +""" +function ScaleNegativeTracers(model::AbstractBiogeochemistry; warn = false) + tracers = conserved_tracers(model) + scalefactors = NamedTuple{tracers}(ones(length(tracers))) return ScaleNegativeTracers(tracers, scalefactors, warn) end -@inline function (scale::ScaleNegativeTracers)(model) +function update_biogeochemical_state!(model, scale::ScaleNegativeTracers) workgroup, worksize = work_layout(model.grid, :xyz) dev = device(model.grid.architecture) scale_for_negs_kernel! = scale_for_negs!(dev, workgroup, worksize) - scale_for_negs_kernel!(values(model.tracers), scale.tracers, scale.scalefactors) + scale_for_negs_kernel!(model.tracers, scale.tracers, scale.scalefactors) end -@inline (scale::ScaleNegativeTracers)(sim::Simulation) = scale(sim.model) -@kernel function _remove_NaN_tendencies!(fields) +@kernel function scale_for_negs!(fields, tracers, scalefactors) i, j, k = @index(Global, NTuple) - for field in fields - if @inbounds isnan(field[i, j, k]) - @inbounds field[i, j, k] = 0.0 - end - end -end - -""" - remove_NaN_tendencies!(model) + t, p = 0.0, 0.0 + @unroll for (idx, tracer) in enumerate(tracers) + field = @inbounds fields[tracer][i, j, k] + scalefactor = @inbounds scalefactors[tracer] -Zeros any `NaN` value tendencies as a final protection against negative tracer run away. -""" -@inline function remove_NaN_tendencies!(model) - workgroup, worksize = work_layout(model.grid, :xyz) - remove_NaN_tendencies_kernel! = _remove_NaN_tendencies!(device(model.grid.architecture), workgroup, worksize) - remove_NaN_tendencies_kernel!(values(model.timestepper.Gⁿ)) + t += field * scalefactor + if field > 0 + p += field * scalefactor + end + end + t < 0 && error("Cell total < 0, can not scale negative tracers.") + @unroll for tracer in tracers + field = @inbounds fields[tracer][i, j, k] + + if field > 0 + field *= t / p + else + field = 0 + end - return nothing -end + @inbounds fields[tracer][i, j, k] = field + end +end \ No newline at end of file diff --git a/src/Utils/sinking_velocity_fields.jl b/src/Utils/sinking_velocity_fields.jl index e4d410eea..2a2805814 100644 --- a/src/Utils/sinking_velocity_fields.jl +++ b/src/Utils/sinking_velocity_fields.jl @@ -30,15 +30,18 @@ adapt_structure(to, velocities::NamedTuple{(:u, :v, :w), Tuple{AbstractField, Ab function show_sinking_velocities(sinking_velocities::NamedTuple{T, V}) where {T, V} str = "" - if length(T) == 1 - str = " └── $(T[1]): $(maximum_sinking(sinking_velocities[1])) to $(minimum_sinking(sinking_velocities[1])) m/s" + if length(T) == 0 + return str + elseif length(T) == 1 + return " └── $(T[1]): $(maximum_sinking(sinking_velocities[1])) to $(minimum_sinking(sinking_velocities[1])) m/s" else for idx in 1:length(T) - 1 str *= " ├── $(T[idx]): $(maximum_sinking(sinking_velocities[idx])) to $(minimum_sinking(sinking_velocities[idx])) m/s \n" end str *= " └── $(T[end]): $(maximum_sinking(sinking_velocities[end])) to $(minimum_sinking(sinking_velocities[end])) m/s" + + return str end - return str end maximum_sinking(velocity) = maximum(velocity) diff --git a/test/test_LOBSTER.jl b/test/test_LOBSTER.jl index c3ef4eeca..db25b6d49 100644 --- a/test/test_LOBSTER.jl +++ b/test/test_LOBSTER.jl @@ -1,18 +1,11 @@ using Test -using OceanBioME: LOBSTER +using OceanBioME: LOBSTER, conserved_tracers, redfield using Oceananigans, CUDA -ΣN(model, variable_redfield) = variable_redfield ? ( - sum(model.tracers.NO₃) + sum(model.tracers.NH₄) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.sPON) + sum(model.tracers.bPON) + sum(model.tracers.DON)) : ( - sum(model.tracers.NO₃) + sum(model.tracers.NH₄) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.sPOM) + sum(model.tracers.bPOM) + sum(model.tracers.DOM)) +ΣN(model, biogeochemistry) = sum([sum(model.tracers[tracer_name]) for tracer_name in conserved_tracers(biogeochemistry)]) -function ΣC(model, carbonates, variable_redfield) - # *will only be conserved if carbonates on* - if variable_redfield - OC = sum(model.tracers.sPOC .+ model.tracers.bPOC .+ model.tracers.DOC) - else - OC = sum(model.tracers.sPOM .+ model.tracers.bPOM .+ model.tracers.DOM) * model.biogeochemistry. organic_redfield - end +function ΣC(model, carbonates, biogeochemistry) + OC = sum([sum(model.tracers[tracer_name] * redfield(Val(tracer_name), biogeochemistry, model.tracers)) for tracer_name in conserved_tracers(biogeochemistry)]) if carbonates IC = sum(model.tracers.DIC) @@ -20,31 +13,21 @@ function ΣC(model, carbonates, variable_redfield) IC = 0.0 end - LC = sum(model.tracers.P * (1 + model.biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.phytoplankton_redfield - - return OC + IC + LC + return OC + IC end -ΣGⁿ(model, variable_redfield) = variable_redfield ? ( - sum(model.timestepper.Gⁿ.NO₃) + sum(model.timestepper.Gⁿ.NH₄) + sum(model.timestepper.Gⁿ.P) + sum(model.timestepper.Gⁿ.Z) + sum(model.timestepper.Gⁿ.sPON) + sum(model.timestepper.Gⁿ.bPON) + sum(model.timestepper.Gⁿ.DON)) : ( - sum(model.timestepper.Gⁿ.NO₃) + sum(model.timestepper.Gⁿ.NH₄) + sum(model.timestepper.Gⁿ.P) + sum(model.timestepper.Gⁿ.Z) + sum(model.timestepper.Gⁿ.sPOM) + sum(model.timestepper.Gⁿ.bPOM) + sum(model.timestepper.Gⁿ.DOM)) +ΣGⁿ(model, biogeochemistry) = sum([sum(model.timestepper.Gⁿ[tracer_name]) for tracer_name in conserved_tracers(biogeochemistry)]) -function ΣGᶜ(model, carbonates, variable_redfield) +function ΣGᶜ(model, biogeochemistry) # *will only be conserved if carbonates on* - if variable_redfield - OC = sum(model.timestepper.Gⁿ.sPOC .+ model.timestepper.Gⁿ.bPOC .+ model.timestepper.Gⁿ.DOC) - else - OC = sum(model.timestepper.Gⁿ.sPOM .+ model.timestepper.Gⁿ.bPOM .+ model.timestepper.Gⁿ.DOM) * model.biogeochemistry. organic_redfield - end + # Leaving this as is for now since as the `redfield` is based on value not tendency then it returns the wrong results if + # the redfield is rapidly changing + OC = sum(model.timestepper.Gⁿ.sPOC .+ model.timestepper.Gⁿ.bPOC .+ model.timestepper.Gⁿ.DOC) - if carbonates - IC = sum(model.timestepper.Gⁿ.DIC) - else - IC = 0.0 - end + IC = sum(model.timestepper.Gⁿ.DIC) - LC = sum(model.timestepper.Gⁿ.P * (1 + model.biogeochemistry.organic_carbon_calcate_ratio) .+ model.timestepper.Gⁿ.Z) * model.biogeochemistry.phytoplankton_redfield + LC = sum(model.timestepper.Gⁿ.P * (1 + biogeochemistry.organic_carbon_calcate_ratio) .+ model.timestepper.Gⁿ.Z) * biogeochemistry.phytoplankton_redfield return OC + IC + LC end @@ -59,7 +42,7 @@ function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open end # correct tracers and auxiliary fields have been setup, and order has not changed - required_tracers = variable_redfield ? (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON) : (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM) + required_tracers = conserved_tracers(model.biogeochemistry) if carbonates required_tracers = (required_tracers..., :DIC, :Alk) end @@ -90,9 +73,9 @@ function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open model.tracers.bPON .= rand() model.tracers.DON .= rand() - model.tracers.sPOC .= model.tracers.sPON * model.biogeochemistry.phytoplankton_redfield - model.tracers.bPOC .= model.tracers.bPON * model.biogeochemistry.phytoplankton_redfield - model.tracers.DOC .= model.tracers.DON * model.biogeochemistry.phytoplankton_redfield + model.tracers.sPOC .= model.tracers.sPON * redfield(Val(:sPOM), model.biogeochemistry) + model.tracers.bPOC .= model.tracers.bPON * redfield(Val(:bPOM), model.biogeochemistry) + model.tracers.DOC .= model.tracers.DON * redfield(Val(:DOM), model.biogeochemistry) else model.tracers.sPOM .= rand() model.tracers.bPOM .= rand() @@ -105,25 +88,25 @@ function test_LOBSTER(grid, carbonates, oxygen, variable_redfield, sinking, open end end - ΣN₀ = CUDA.@allowscalar ΣN(model, variable_redfield) + ΣN₀ = CUDA.@allowscalar ΣN(model, model.biogeochemistry) - ΣC₀ = CUDA.@allowscalar ΣC(model, carbonates, variable_redfield) + ΣC₀ = CUDA.@allowscalar ΣC(model, carbonates, model.biogeochemistry) for _ in 1:n_timesteps time_step!(model, 1.0) end - ΣN₁ = CUDA.@allowscalar ΣN(model, variable_redfield) + ΣN₁ = CUDA.@allowscalar ΣN(model, model.biogeochemistry) CUDA.@allowscalar if !(sinking && open_bottom) #when we have open bottom sinking we won't conserve anything @test ΣN₀ ≈ ΣN₁ - @test ΣGⁿ(model, variable_redfield) ≈ 0.0 atol = 1e-15 # rtol=sqrt(eps) so is usually much larger than even this + @test ΣGⁿ(model, model.biogeochemistry) ≈ 0.0 atol = 1e-15 # rtol=sqrt(eps) so is usually much larger than even this if (carbonates && variable_redfield) - ΣC₁ = ΣC(model, carbonates, variable_redfield) - @test ΣC₀ ≈ ΣC₁# atol = 0.0001 # when we convert to and from + ΣC₁ = ΣC(model, carbonates, model.biogeochemistry) + @test ΣC₀ ≈ ΣC₁ - @test ΣGᶜ(model, carbonates, variable_redfield) ≈ 0.0 atol = 1e-15 + @test ΣGᶜ(model, model.biogeochemistry.underlying_biogeochemistry) ≈ 0.0 atol = 1e-15 end end return model diff --git a/test/test_light.jl b/test/test_light.jl index 7be4c12af..fdcf99169 100644 --- a/test/test_light.jl +++ b/test/test_light.jl @@ -1,6 +1,6 @@ using Oceananigans, Test using OceanBioME: TwoBandPhotosyntheticallyActiveRadiation, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus -using Oceananigans.Biogeochemistry: update_biogeochemical_state!, required_biogeochemical_tracers +using Oceananigans.Biogeochemistry: update_biogeochemical_state!, required_biogeochemical_tracers, biogeochemical_auxiliary_fields Pᵢ(x,y,z) = 2.5 + z @@ -17,7 +17,7 @@ function test_two_band(grid, bgc, model_type) update_biogeochemical_state!(model.biogeochemistry, model) - PAR_model = model.biogeochemistry.light_attenuation_model + PAR_model = model.biogeochemistry.light_attenuation kʳ = PAR_model.water_red_attenuation kᵇ = PAR_model.water_blue_attenuation @@ -37,7 +37,7 @@ function test_two_band(grid, bgc, model_type) expected_PAR = 100.0 .* [exp(- 0.5 * kʳ - ∫Chlʳ[1] * χʳ) + exp(- 0.5 * kᵇ - ∫Chlᵇ[1] * χᵇ), exp(- 1.5 * kʳ - ∫Chlʳ[2] * χʳ) + exp(- 1.5 * kᵇ - ∫Chlᵇ[2] * χᵇ)] ./ 2 - results_PAR = convert(Array, model.biogeochemistry.light_attenuation_model.field)[1, 1, 1:2] + results_PAR = convert(Array, biogeochemical_auxiliary_fields(biogeochemistry).PAR)[1, 1, 1:2] return all(results_PAR .≈ reverse(expected_PAR)) end @@ -49,7 +49,7 @@ archs = (CPU(), ) arch in archs, grid in (RectilinearGrid(arch; size = (2, 2, 2), extent = (2, 2, 2)), LatitudeLongitudeGrid(arch; size = (5, 5, 2), longitude = (-180, 180), latitude = (-85, 85), z = (-2, 0))), - bgc in (LOBSTER, NutrientPhytoplanktonZooplanktonDetritus) + bgc in (LOBSTER, NutrientPhytoplanktonZooplanktonDetritus) # this is now redundant since each model doesn't deal with the light separatly if !((model == NonhydrostaticModel) && ((grid isa LatitudeLongitudeGrid) | (grid isa OrthogonalSphericalShellGrid))) @info "Testing $bgc in $model on $grid..." diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 6e48e650c..73c7f940a 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -68,9 +68,9 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd buoyancy = nothing, tracers = nothing) - set_defaults!(model.biogeochemistry.sediment_model) + set_defaults!(model.biogeochemistry.sediment) - set_defaults!(biogeochemistry, model) + set_defaults!(biogeochemistry.underlying_biogeochemistry, model) simulation = Simulation(model, Δt = 50, stop_time = 1day) @@ -78,7 +78,7 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd simulation.callbacks[:intercept_tendencies] = Callback(intercept_tendencies!; callsite = TendencyCallsite(), parameters = intercepted_tendencies) - N₀ = total_nitrogen(biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(model.biogeochemistry.sediment_model) * Azᶠᶜᶜ(1, 1, 1, grid) + N₀ = total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) run!(simulation) @@ -86,14 +86,14 @@ function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAd @test any([any(intercepted_tendencies[tracer] .!= model.timestepper.Gⁿ[tracer]) for tracer in keys(model.tracers)]) # the sediment tendencies are being updated - @test all([any(tend .!= 0.0) for tend in model.biogeochemistry.sediment_model.tendencies.Gⁿ]) - @test all([any(tend .!= 0.0) for tend in model.biogeochemistry.sediment_model.tendencies.G⁻]) + @test all([any(tend .!= 0.0) for tend in model.biogeochemistry.sediment.tendencies.Gⁿ]) + @test all([any(tend .!= 0.0) for tend in model.biogeochemistry.sediment.tendencies.G⁻]) # the sediment values are being integrated initial_values = (N_fast = 0.0230, N_slow = 0.0807, C_fast = 0.5893, C_slow = 0.1677, N_ref = 0.0, C_ref = 0.0, N_storage = 0.0) - @test all([any(field .!= initial_values[name]) for (name, field) in pairs(model.biogeochemistry.sediment_model.fields)]) + @test all([any(field .!= initial_values[name]) for (name, field) in pairs(model.biogeochemistry.sediment.fields)]) - N₁ = total_nitrogen(biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(model.biogeochemistry.sediment_model) * Azᶠᶜᶜ(1, 1, 1, grid) + N₁ = total_nitrogen(biogeochemistry.underlying_biogeochemistry, model) * volume(1, 1, 1, grid, Center(), Center(), Center()) + total_nitrogen(biogeochemistry.sediment) * Azᶠᶜᶜ(1, 1, 1, grid) # conservations @test N₁ ≈ N₀ @@ -132,10 +132,10 @@ bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill # get rid of incompatible combinations run = ifelse((model == NonhydrostaticModel && (isa(grid, ImmersedBoundaryGrid) || isa(grid, LatitudeLongitudeGrid))) || (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) || - (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) + (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry.underlying_biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) if run - @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) in $(display_name(biogeochemistry)) on $(display_name(grid)) with $(model)" - @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry)), $(display_name(grid)), $(model)" test_flat_sediment(grid, biogeochemistry, model; timestepper) + @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry.underlying_biogeochemistry))" + @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry.underlying_biogeochemistry))" test_flat_sediment(grid, biogeochemistry, model; timestepper) end end end diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index d9738887a..ada407cd9 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -40,7 +40,7 @@ end initial_kelp_N = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.N .+ particles.structural_nitrogen)) ./ (14 * 0.001) initial_tracer_C = sum(model.tracers.sPOC) + sum(model.tracers.bPOC) + sum(model.tracers.DOC) + sum(model.tracers.DIC) + - sum(model.tracers.P * (1 + model.biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.phytoplankton_redfield + sum(model.tracers.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield initial_kelp_C = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.C .+ particles.structural_carbon)) ./ (12 * 0.001) @@ -54,7 +54,7 @@ end final_kelp_N = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.N .+ particles.structural_nitrogen)) ./ (14 * 0.001) final_tracer_C = sum(model.tracers.sPOC) + sum(model.tracers.bPOC) + sum(model.tracers.DOC) + sum(model.tracers.DIC) + - sum(model.tracers.P * (1 + model.biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.phytoplankton_redfield + sum(model.tracers.P * (1 + model.biogeochemistry.underlying_biogeochemistry.organic_carbon_calcate_ratio) .+ model.tracers.Z) * model.biogeochemistry.underlying_biogeochemistry.phytoplankton_redfield final_kelp_C = sum(particles.A .* particles.structural_dry_weight_per_area .* (particles.C .+ particles.structural_carbon)) ./ (12 * 0.001) # kelp is being integrated diff --git a/test/test_utils.jl b/test/test_utils.jl index 1245c8f69..86586a791 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -20,48 +20,35 @@ end function test_negative_scaling(arch) grid = RectilinearGrid(arch, size = (1, 1, 1), extent = (1, 1, 1)) - model = NonhydrostaticModel(; grid, tracers = (:A, :B)) + model = NonhydrostaticModel(; grid, biogeochemistry = NutrientPhytoplanktonZooplanktonDetritus(; grid, scale_negatives = true)) - set!(model, A = 1, B = -0.5) + set!(model, N = 2, P = -1) - simulation = Simulation(model, Δt = 1) + simulation = Simulation(model, Δt = 1e-10, stop_iteration = 1) - negativity_protection! = ScaleNegativeTracers(; model, tracers = (:A, :B)) + run!(simulation) - negativity_protection!(simulation) - - return (model.tracers.A[1, 1, 1] == 0.5) && (model.tracers.B[1, 1, 1] == 0.0) -end - -function test_negative_scaling(arch) - grid = RectilinearGrid(arch, size = (1, 1, 1), extent = (1, 1, 1)) - - model = NonhydrostaticModel(; grid, tracers = :A) - - model.timestepper.Gⁿ.A .= NaN - - remove_NaN_tendencies!(model) - - return model.timestepper.Gⁿ.A[1, 1, 1] == 0.0 + return (model.tracers.N[1, 1, 1] ≈ 1) && (model.tracers.P[1, 1, 1] ≈ 0.0) end function test_negative_zeroing(arch) grid = RectilinearGrid(arch, size = (1, 1, 1), extent = (1, 1, 1)) - model = NonhydrostaticModel(; grid, tracers = (:A, :B)) + model = NonhydrostaticModel(; grid, biogeochemistry = NutrientPhytoplanktonZooplanktonDetritus(; grid, modifiers = ZeroNegativeTracers(; exclude = (:Z, )))) - set!(model, A = -1, B = -0.5) + set!(model, N = 2, P = -1, Z = -1) - zero_negative_tracers!(model, params = (exclude = (:B, ), )) + simulation = Simulation(model, Δt = 1e-10, stop_iteration = 1) - return (model.tracers.A[1, 1, 1] == 0.0) && (model.tracers.B[1, 1, 1] == -0.5) + run!(simulation) + + return (model.tracers.N[1, 1, 1] ≈ 2) && (model.tracers.P[1, 1, 1] ≈ 0.0) && (model.tracers.Z[1, 1, 1] ≈ -1) end @testset "Test Utils" begin for arch in (CPU(), ) @test test_column_diffusion_timescale(arch) @test test_negative_scaling(arch) - @test test_negative_scaling(arch) @test test_negative_zeroing(arch) end end