diff --git a/benchmark/Manifest.toml b/benchmark/Manifest.toml index d6f4771633..460c99cc3f 100644 --- a/benchmark/Manifest.toml +++ b/benchmark/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.9" manifest_format = "2.0" -project_hash = "1a00cfe2c9a8fb2818c5ca2754e33ad2254ad429" +project_hash = "42cecfba6d42b9c8f227e7cedee3d6880a3cef7c" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] @@ -72,10 +72,10 @@ version = "1.1.1" oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b" [[deps.BFloat16s]] -deps = ["LinearAlgebra", "Printf", "Random", "Test"] -git-tree-sha1 = "2c7cc21e8678eff479978a0a2ef5ce2f51b63dff" +deps = ["LinearAlgebra", "Printf", "Random"] +git-tree-sha1 = "3b642331600250f592719140c60cf12372b82d66" uuid = "ab4f0b2a-ad5b-11e8-123f-65d77653426b" -version = "0.5.0" +version = "0.5.1" [[deps.BSON]] git-tree-sha1 = "4c3e506685c527ac6a54ccc0c8c76fd6f91b42fb" @@ -109,9 +109,9 @@ version = "0.5.0" [[deps.CUDA]] deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "GPUToolbox", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "StaticArrays", "Statistics", "demumble_jll"] -git-tree-sha1 = "7f393da6f01f271d3791d06f2916db8cc7d999f3" +git-tree-sha1 = "acfe6dc14594afba70c002e3b4bfc2f1f8273347" uuid = "052768ef-5323-5732-b1bb-66c8b64840ba" -version = "5.7.0" +version = "5.7.2" [deps.CUDA.extensions] ChainRulesCoreExt = "ChainRulesCore" @@ -161,9 +161,9 @@ version = "3.29.0" [[deps.ColorTypes]] deps = ["FixedPointNumbers", "Random"] -git-tree-sha1 = "c7acce7a7e1078a20a285211dd73cd3941a871d6" +git-tree-sha1 = "67e11ee83a43eb71ddc950302c53bf33f0690dfe" uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" -version = "0.12.0" +version = "0.12.1" [deps.ColorTypes.extensions] StyledStringsExt = "StyledStrings" @@ -281,9 +281,9 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" [[deps.Dbus_jll]] deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl"] -git-tree-sha1 = "fc173b380865f70627d7dd1190dc2fce6cc105af" +git-tree-sha1 = "473e9afc9cf30814eb67ffa5f2db7df82c3ad9fd" uuid = "ee1fde0b-3d02-5ea6-8484-8dfef6360eab" -version = "1.14.10+0" +version = "1.16.2+0" [[deps.DelimitedFiles]] deps = ["Mmap"] @@ -310,10 +310,9 @@ deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" [[deps.DocStringExtensions]] -deps = ["LibGit2"] -git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" +git-tree-sha1 = "e7b7e6f178525d17c720ab9c081e4ef04429f860" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.9.3" +version = "0.9.4" [[deps.Downloads]] deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] @@ -343,6 +342,11 @@ git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec" uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" version = "0.1.10" +[[deps.ExpressionExplorer]] +git-tree-sha1 = "bf2e6a47b70dfb5d103f300ef83d950239f9fa50" +uuid = "21656369-7473-754a-2065-74616d696c43" +version = "1.1.2" + [[deps.FFMPEG]] deps = ["FFMPEG_jll"] git-tree-sha1 = "53ebe7511fa11d33bec688a9178fac4e49eeee00" @@ -362,10 +366,10 @@ uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" version = "1.8.1" [[deps.FFTW_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "4d81ed14783ec49ce9f2e168208a12ce1815aa25" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "6d6219a004b8cf1e0b4dbe27a2860b8e04eba0be" uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" -version = "3.3.10+3" +version = "3.3.11+0" [[deps.FileIO]] deps = ["Pkg", "Requires", "UUIDs"] @@ -388,9 +392,9 @@ version = "0.8.5" [[deps.Fontconfig_jll]] deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"] -git-tree-sha1 = "21fac3c77d7b5a9fc03b0ec503aa1a6392c34d2b" +git-tree-sha1 = "301b5d5d731a0654825f1f2e906990f7141a106b" uuid = "a3f928ae-7b40-5064-980b-68af3947d34b" -version = "2.15.0+0" +version = "2.16.0+0" [[deps.Format]] git-tree-sha1 = "9c68794ef81b08086aeb32eeaf33531668d5f5fc" @@ -399,15 +403,15 @@ version = "1.3.7" [[deps.FreeType2_jll]] deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"] -git-tree-sha1 = "786e968a8d2fb167f2e4880baba62e0e26bd8e4e" +git-tree-sha1 = "2c5512e11c791d1baed2049c5652441b28fc6a31" uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7" -version = "2.13.3+1" +version = "2.13.4+0" [[deps.FriBidi_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "846f7026a9decf3679419122b49f8a1fdb48d2d5" +git-tree-sha1 = "7a214fdac5ed5f59a22c2d9a885a16da1c74bbc7" uuid = "559328eb-81f9-559d-9380-de523a88c83c" -version = "1.0.16+0" +version = "1.0.17+0" [[deps.Future]] deps = ["Random"] @@ -433,14 +437,14 @@ version = "0.2.0" [[deps.GPUCompiler]] deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "PrecompileTools", "Preferences", "Scratch", "Serialization", "TOML", "TimerOutputs", "UUIDs"] -git-tree-sha1 = "199f213e40a7982e9138bc9edc3299419d510390" +git-tree-sha1 = "b08c164134dd0dbc76ff54e45e016cf7f30e16a4" uuid = "61eb1bfa-7361-4325-ad38-22787b887f55" -version = "1.2.0" +version = "1.3.2" [[deps.GPUToolbox]] -git-tree-sha1 = "52ad2902dda44ab527891317e144e7f718ab77d7" +git-tree-sha1 = "15d8b0f5a6dca9bf8c02eeaf6687660dafa638d0" uuid = "096a3bc2-3ced-46d0-87f4-dd12716f4bfc" -version = "0.1.0" +version = "0.2.0" [[deps.GR]] deps = ["Artifacts", "Base64", "DelimitedFiles", "Downloads", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Preferences", "Printf", "Qt6Wayland_jll", "Random", "Serialization", "Sockets", "TOML", "Tar", "Test", "p7zip_jll"] @@ -472,10 +476,10 @@ uuid = "c27321d9-0574-5035-807b-f59d2c89b15c" version = "1.3.1" [[deps.Graphite2_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "01979f9b37367603e2848ea225918a3b3861b606" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "8a6dbda1fd736d60cc477d99f2e7a042acfa46e8" uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472" -version = "1.3.14+1" +version = "1.3.15+0" [[deps.Grisu]] git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2" @@ -484,9 +488,9 @@ version = "1.0.2" [[deps.HTTP]] deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "PrecompileTools", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] -git-tree-sha1 = "c67b33b085f6e2faf8bf79a61962e7339a81129c" +git-tree-sha1 = "f93655dc73d7a0b4a368e3c0bce296ae035ad76e" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "1.10.15" +version = "1.10.16" [[deps.HarfBuzz_jll]] deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll"] @@ -505,12 +509,6 @@ git-tree-sha1 = "f93a9ce66cd89c9ba7a4695a47fd93b4c6bc59fa" uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8" version = "2.12.0+0" -[[deps.IncompleteLU]] -deps = ["LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "6c676e79f98abb6d33fa28122cad099f1e464afe" -uuid = "40713840-3770-5561-ab4c-a76e7d0d7895" -version = "0.2.1" - [[deps.Inflate]] git-tree-sha1 = "d1b1b796e47d94588b3757fe84fbf65a5ec4a80d" uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" @@ -562,15 +560,15 @@ version = "1.0.0" [[deps.JLD2]] deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "PrecompileTools", "Requires", "TranscodingStreams"] -git-tree-sha1 = "91d501cb908df6f134352ad73cde5efc50138279" +git-tree-sha1 = "1059c071429b4753c0c869b75c859c44ba09a526" uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -version = "0.5.11" +version = "0.5.12" [[deps.JLFzf]] -deps = ["Pipe", "REPL", "Random", "fzf_jll"] -git-tree-sha1 = "71b48d857e86bf7a1838c4736545699974ce79a2" +deps = ["REPL", "Random", "fzf_jll"] +git-tree-sha1 = "1d4015b1eb6dc3be7e6c400fbd8042fe825a6bac" uuid = "1019f520-868f-41f5-a6de-eb00f4b6a39c" -version = "0.1.9" +version = "0.1.10" [[deps.JLLWrappers]] deps = ["Artifacts", "Preferences"] @@ -584,18 +582,6 @@ git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" version = "0.21.4" -[[deps.JSON3]] -deps = ["Dates", "Mmap", "Parsers", "PrecompileTools", "StructTypes", "UUIDs"] -git-tree-sha1 = "1d322381ef7b087548321d3f878cb4c9bd8f8f9b" -uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" -version = "1.14.1" - - [deps.JSON3.extensions] - JSON3ArrowExt = ["ArrowTypes"] - - [deps.JSON3.weakdeps] - ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd" - [[deps.JpegTurbo_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] git-tree-sha1 = "eac1206917768cb54957c65a615460d87b455fc1" @@ -624,6 +610,12 @@ version = "0.9.34" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +[[deps.Krylov]] +deps = ["LinearAlgebra", "Printf", "SparseArrays"] +git-tree-sha1 = "b29d37ce30fa401a4563b18880ab91f979a29734" +uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" +version = "0.9.10" + [[deps.KrylovPreconditioners]] deps = ["Adapt", "KernelAbstractions", "LightGraphs", "LinearAlgebra", "Metis", "SparseArrays"] git-tree-sha1 = "517b267073c44f494357507e37650d8516adf100" @@ -692,9 +684,9 @@ version = "1.4.0" [[deps.Latexify]] deps = ["Format", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "OrderedCollections", "Requires"] -git-tree-sha1 = "cd714447457c660382fe634710fb56eb255ee42e" +git-tree-sha1 = "cd10d2cc78d34c0e2a3a36420ab607b611debfbb" uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" -version = "0.16.6" +version = "0.16.7" [deps.Latexify.extensions] DataFramesExt = "DataFrames" @@ -751,9 +743,9 @@ version = "3.2.2+2" [[deps.Libgcrypt_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll"] -git-tree-sha1 = "8be878062e0ffa2c3f67bb58a595375eda5de80b" +git-tree-sha1 = "d77592fa54ad343c5043b6f38a03f1a3c3959ffe" uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4" -version = "1.11.0+0" +version = "1.11.1+0" [[deps.Libglvnd_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll", "Xorg_libXext_jll"] @@ -775,9 +767,9 @@ version = "1.18.0+0" [[deps.Libmount_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "89211ea35d9df5831fca5d33552c02bd33878419" +git-tree-sha1 = "a31572773ac1b745e0343fe5e2c8ddda7a37e997" uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9" -version = "2.40.3+0" +version = "2.41.0+0" [[deps.Libtiff_jll]] deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"] @@ -787,9 +779,9 @@ version = "4.7.1+0" [[deps.Libuuid_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "e888ad02ce716b319e6bdb985d2ef300e7089889" +git-tree-sha1 = "321ccef73a96ba828cd51f2ab5b9f917fa73945a" uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" -version = "2.40.3+0" +version = "2.41.0+0" [[deps.LightGraphs]] deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] @@ -950,28 +942,31 @@ version = "3.1.1+0" [[deps.NaNMath]] deps = ["OpenLibm_jll"] -git-tree-sha1 = "cc0a5deefdb12ab3a096f00a6d42133af4560d71" +git-tree-sha1 = "9b8215b1ee9e78a293f99797cd31375471b2bcae" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" -version = "1.1.2" +version = "1.1.3" [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" version = "1.2.0" [[deps.Oceananigans]] -deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "GPUArrays", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "OffsetArrays", "OrderedCollections", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "StaticArrays", "Statistics", "StructArrays", "TimesDates"] -git-tree-sha1 = "abb9ab23cf36df4e24bb1837e3dc64efd5a404a4" +deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "GPUArrays", "Glob", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "Krylov", "KrylovPreconditioners", "LinearAlgebra", "Logging", "MPI", "OffsetArrays", "OrderedCollections", "Pkg", "Printf", "Random", "ReactantCore", "Rotations", "SeawaterPolynomials", "SparseArrays", "StaticArrays", "Statistics", "StructArrays", "TimesDates"] +path = ".." uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" -version = "0.96.0" +version = "0.96.20" [deps.Oceananigans.extensions] + OceananigansAMDGPUExt = "AMDGPU" OceananigansEnzymeExt = "Enzyme" OceananigansMakieExt = ["MakieCore", "Makie"] OceananigansMetalExt = "Metal" OceananigansNCDatasetsExt = "NCDatasets" + OceananigansOneAPIExt = "oneAPI" OceananigansReactantExt = ["Reactant", "KernelAbstractions", "ConstructionBase"] [deps.Oceananigans.weakdeps] + AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" @@ -979,6 +974,7 @@ version = "0.96.0" Metal = "dde4c033-4e86-420c-a63e-0dd931031962" NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" Reactant = "3c362404-f566-11ee-1572-e11a4b42c853" + oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b" [[deps.OffsetArrays]] git-tree-sha1 = "a414039192a155fb38c4599a60110f0018c6ec82" @@ -1007,9 +1003,9 @@ version = "0.8.1+4" [[deps.OpenMPI_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML", "Zlib_jll"] -git-tree-sha1 = "da913f03f17b449951e0461da960229d4a3d1a8c" +git-tree-sha1 = "047b66eb62f3cae59ed260ebb9075a32a04350f1" uuid = "fe0851c0-eecd-5654-98d4-656369965a5c" -version = "5.0.7+1" +version = "5.0.7+2" [[deps.OpenSSL]] deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"] @@ -1051,11 +1047,6 @@ git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" version = "2.8.1" -[[deps.Pipe]] -git-tree-sha1 = "6842804e7867b115ca9de748a0cf6b364523c16d" -uuid = "b98c9c47-44ae-5843-9183-064241ee97a0" -version = "1.3.0" - [[deps.Pixman_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"] git-tree-sha1 = "db76b1ecd5e9715f3d043cec13b2ec93ce015d53" @@ -1069,9 +1060,9 @@ version = "1.10.0" [[deps.PkgBenchmark]] deps = ["BenchmarkTools", "Dates", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Pkg", "Printf", "TerminalLoggers", "UUIDs"] -git-tree-sha1 = "e4a10b7cdb7ec836850e43a4cee196f4e7b02756" +git-tree-sha1 = "dc06838e215d9c072fbdb9e9e5497d7f476ccaf5" uuid = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d" -version = "0.2.12" +version = "0.2.13" [[deps.PkgVersion]] deps = ["Pkg"] @@ -1216,6 +1207,12 @@ git-tree-sha1 = "c6ec94d2aaba1ab2ff983052cf6a606ca5985902" uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143" version = "1.6.0" +[[deps.ReactantCore]] +deps = ["ExpressionExplorer", "MacroTools"] +git-tree-sha1 = "9189b6f86857f548a8981c99cc5e0df01baeb6cf" +uuid = "a3311ec8-5e00-46d5-b541-4f83e724a433" +version = "0.1.9" + [[deps.RealDot]] deps = ["LinearAlgebra"] git-tree-sha1 = "9f0a1b71baaf7650f4fa8a1d168c7fb6ee41f0c9" @@ -1278,9 +1275,9 @@ uuid = "6c6a2e73-6563-6170-7368-637461726353" version = "1.2.1" [[deps.SeawaterPolynomials]] -git-tree-sha1 = "78f965a2f0cd5250a20c9aba9979346dd2b35734" +git-tree-sha1 = "e2671e9abe2a2faa51dcecd9d911522931c16012" uuid = "d496a93d-167e-4197-9f49-d3af4ff8fe40" -version = "0.3.5" +version = "0.3.10" [[deps.SentinelArrays]] deps = ["Dates", "Random"] @@ -1376,9 +1373,9 @@ version = "0.4.1" [[deps.StructArrays]] deps = ["ConstructionBase", "DataAPI", "Tables"] -git-tree-sha1 = "5a3a31c41e15a1e042d60f2f4942adccba05d3c9" +git-tree-sha1 = "8ad2e38cbb812e29348719cc63580ec1dfeb9de4" uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" -version = "0.7.0" +version = "0.7.1" weakdeps = ["Adapt", "GPUArraysCore", "KernelAbstractions", "LinearAlgebra", "SparseArrays", "StaticArrays"] [deps.StructArrays.extensions] @@ -1388,12 +1385,6 @@ weakdeps = ["Adapt", "GPUArraysCore", "KernelAbstractions", "LinearAlgebra", "Sp StructArraysSparseArraysExt = "SparseArrays" StructArraysStaticArraysExt = "StaticArrays" -[[deps.StructTypes]] -deps = ["Dates", "UUIDs"] -git-tree-sha1 = "159331b30e94d7b11379037feeb9b690950cace8" -uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4" -version = "1.11.0" - [[deps.SuiteSparse_jll]] deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" @@ -1406,9 +1397,9 @@ version = "1.0.3" [[deps.TZJData]] deps = ["Artifacts"] -git-tree-sha1 = "7def47e953a91cdcebd08fbe76d69d2715499a7d" +git-tree-sha1 = "72df96b3a595b7aab1e101eb07d2a435963a97e2" uuid = "dc5dba14-91b3-4cab-a142-028a31da12f7" -version = "1.4.0+2025a" +version = "1.5.0+2025b" [[deps.TableTraits]] deps = ["IteratorInterfaceExtensions"] @@ -1428,10 +1419,10 @@ uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" version = "1.10.0" [[deps.TaylorSeries]] -deps = ["LinearAlgebra", "Markdown", "Requires", "SparseArrays"] -git-tree-sha1 = "ae73e40c647c0061697fa9708ee12cce385653e3" +deps = ["LinearAlgebra", "Markdown", "SparseArrays"] +git-tree-sha1 = "f80d024c30771a57bc173d08b566379d04af8055" uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" -version = "0.18.3" +version = "0.18.5" [deps.TaylorSeries.extensions] TaylorSeriesIAExt = "IntervalArithmetic" @@ -1495,9 +1486,9 @@ uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" version = "0.11.3" [[deps.URIs]] -git-tree-sha1 = "67db6cc7b3821e19ebe75791a9dd19c9b1188f2b" +git-tree-sha1 = "cbbebadbcc76c5ca1cc4b4f3b0614b3e603b5000" uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4" -version = "1.5.1" +version = "1.5.2" [[deps.UUIDs]] deps = ["Random", "SHA"] @@ -1577,21 +1568,21 @@ version = "2.13.6+1" [[deps.XSLT_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "XML2_jll", "Zlib_jll"] -git-tree-sha1 = "7d1671acbe47ac88e981868a078bd6b4e27c5191" +git-tree-sha1 = "82df486bfc568c29de4a207f7566d6716db6377c" uuid = "aed1982a-8fda-507f-9586-7b0439959a61" -version = "1.1.42+0" +version = "1.1.43+0" [[deps.XZ_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "56c6604ec8b2d82cc4cfe01aa03b00426aac7e1f" +git-tree-sha1 = "fee71455b0aaa3440dfdd54a9a36ccef829be7d4" uuid = "ffd25f8a-64ca-5728-b0f7-c24cf3aae800" -version = "5.6.4+1" +version = "5.8.1+0" [[deps.Xorg_libICE_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "326b4fea307b0b39892b3e85fa451692eda8d46c" +git-tree-sha1 = "a3ea76ee3f4facd7a64684f9af25310825ee3668" uuid = "f67eecfb-183a-506d-b269-f58e52b52d7c" -version = "1.1.1+0" +version = "1.1.2+0" [[deps.Xorg_libSM_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libICE_jll"] @@ -1607,9 +1598,9 @@ version = "1.8.6+3" [[deps.Xorg_libXau_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "e9216fdcd8514b7072b43653874fd688e4c6c003" +git-tree-sha1 = "aa1261ebbac3ccc8d16558ae6799524c450ed16b" uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec" -version = "1.0.12+0" +version = "1.0.13+0" [[deps.Xorg_libXcursor_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"] @@ -1727,9 +1718,9 @@ version = "2.39.0+0" [[deps.Xorg_xtrans_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "6dba04dbfb72ae3ebe5418ba33d087ba8aa8cb00" +git-tree-sha1 = "a63799ff68005991f9d9491b6e95bd3478d783cb" uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10" -version = "1.5.1+0" +version = "1.6.0+0" [[deps.Zlib_jll]] deps = ["Libdl"] diff --git a/benchmark/Project.toml b/benchmark/Project.toml index b7893ec2bf..ffb2dac5ee 100644 --- a/benchmark/Project.toml +++ b/benchmark/Project.toml @@ -10,7 +10,6 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" KrylovPreconditioners = "45d422c2-293f-44ce-8315-2cb988662dec" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" diff --git a/benchmark/benchmark_hydrostatic_model.jl b/benchmark/benchmark_hydrostatic_model.jl index 58381b858b..d3e70127eb 100644 --- a/benchmark/benchmark_hydrostatic_model.jl +++ b/benchmark/benchmark_hydrostatic_model.jl @@ -3,16 +3,16 @@ push!(LOAD_PATH, joinpath(@__DIR__, "..")) using Benchmarks using Oceananigans.TimeSteppers: update_state! -using Oceananigans.Diagnostics: accurate_cell_advection_timescale - using BenchmarkTools using CUDA using Oceananigans +using Oceananigans.Models.HydrostaticFreeSurfaceModels: DiagonallyDominantInversePreconditioner using Statistics +using Oceananigans.Solvers # Problem size -Nx = 256 -Ny = 128 +Nx = 512 +Ny = 256 function set_divergent_velocity!(model) # Create a divergent velocity @@ -33,31 +33,23 @@ function set_divergent_velocity!(model) return nothing end + +random_vector = - 0.5 .* rand(Nx, Ny) .- 0.5 +bottom_height(arch) = GridFittedBottom(Oceananigans.on_architecture(arch, random_vector)) +rgrid(arch) = RectilinearGrid(arch, size=(Nx, Ny, 1), extent=(1, 1, 1), halo = (3, 3, 3)) +lgrid(arch) = LatitudeLongitudeGrid(arch, size=(Nx, Ny, 1), longitude=(-180, 180), latitude=(-80, 80), z=(-1, 0), halo = (3, 3, 3)) grids = Dict( - (CPU, :RectilinearGrid) => RectilinearGrid(CPU(), size=(Nx, Ny, 1), extent=(1, 1, 1)), - (CPU, :LatitudeLongitudeGrid) => LatitudeLongitudeGrid(CPU(), size=(Nx, Ny, 1), longitude=(-180, 180), latitude=(-80, 80), z=(-1, 0), precompute_metrics=true), - (GPU, :RectilinearGrid) => RectilinearGrid(GPU(), size=(Nx, Ny, 1), extent=(1, 1, 1)), - (GPU, :LatitudeLongitudeGrid) => LatitudeLongitudeGrid(GPU(), size=(Nx, Ny, 1), longitude=(-160, 160), latitude=(-80, 80), z=(-1, 0), precompute_metrics=true), + (CPU, :RectilinearGrid) => rgrid(CPU()), + (CPU, :LatitudeLongitudeGrid) => lgrid(CPU()), + (CPU, :ImmersedRecGrid) => ImmersedBoundaryGrid(rgrid(CPU()), bottom_height(GPU())), + (CPU, :ImmersedLatGrid) => ImmersedBoundaryGrid(lgrid(CPU()), bottom_height(GPU())), + (GPU, :RectilinearGrid) => rgrid(GPU()), + (GPU, :LatitudeLongitudeGrid) => lgrid(GPU()), + (GPU, :ImmersedRecGrid) => ImmersedBoundaryGrid(rgrid(GPU()), bottom_height(GPU())), + (GPU, :ImmersedLatGrid) => ImmersedBoundaryGrid(lgrid(CPU()), bottom_height(GPU())) ) -# Cubed sphere cases maybe worth considering eventually -#= -using DataDeps - -ENV["DATADEPS_ALWAYS_ACCEPT"] = "true" - -dd = DataDep("cubed_sphere_510_grid", - "Conformal cubed sphere grid with 510×510 grid points on each face", - "https://engaging-web.mit.edu/~alir/cubed_sphere_grids/cs510/cubed_sphere_510_grid.jld2" -) - -DataDeps.register(dd) - -# (CPU, :ConformalCubedSphereGrid) => ConformalCubedSphereGrid(datadep"cubed_sphere_510_grid/cubed_sphere_510_grid.jld2", Nz=1, z=(-1, 0)), -# (GPU, :ConformalCubedSphereGrid) => ConformalCubedSphereGrid(datadep"cubed_sphere_510_grid/cubed_sphere_510_grid.jld2", Nz=1, z=(-1, 0), architecture=GPU()), -=# - free_surfaces = Dict( :ExplicitFreeSurface => ExplicitFreeSurface(), :PCGImplicitFreeSurface => ImplicitFreeSurface(solver_method = :PreconditionedConjugateGradient), @@ -74,11 +66,11 @@ function benchmark_hydrostatic_model(Arch, grid_type, free_surface_type) grid = grids[(Arch, grid_type)] model = HydrostaticFreeSurfaceModel(; grid, - momentum_advection = VectorInvariant(), - free_surface = free_surfaces[free_surface_type]) + momentum_advection = VectorInvariant(), + free_surface = free_surfaces[free_surface_type]) set_divergent_velocity!(model) - Δt = accurate_cell_advection_timescale(grid, model.velocities) / 2 + Δt = Oceananigans.Advection.cell_advection_timescale(grid, model.velocities) / 2 time_step!(model, Δt) # warmup trial = @benchmark begin @@ -90,15 +82,13 @@ end # Benchmark parameters -#architectures = has_cuda() ? [GPU] : [CPU] -architectures = [CPU] +architectures = has_cuda() ? [GPU, CPU] : [CPU] grid_types = [ - :RectilinearGrid, - :LatitudeLongitudeGrid, - # Uncomment when OrthogonalSphericalShellGrids of any size can be built natively without loading from file: - # :OrthogonalSphericalShellGrid, - # :ConformalCubedSphereGrid + :RectilinearGrid, + :LatitudeLongitudeGrid, + :ImmersedRecGrid, + :ImmersedLatGrid, ] free_surface_types = collect(keys(free_surfaces)) diff --git a/benchmark/benchmark_performance_issue.jl b/benchmark/benchmark_performance_issue.jl new file mode 100644 index 0000000000..36ef660e64 --- /dev/null +++ b/benchmark/benchmark_performance_issue.jl @@ -0,0 +1,128 @@ +push!(LOAD_PATH, joinpath(@__DIR__, "..")) + +using Benchmarks + +using Oceananigans.TimeSteppers: update_state! +using BenchmarkTools +using CUDA +using Oceananigans +using Oceananigans.Operators +using Oceananigans.TurbulenceClosures +using Statistics +using Oceananigans.Solvers +using SeawaterPolynomials.TEOS10 +using NVTX + +CUDA.device!(2) + +function ocean_benchmark(grid, closure) + momentum_advection = nothing # WENOVectorInvariant() + tracer_advection = WENO(order=7) + buoyancy = SeawaterBuoyancy(equation_of_state=TEOS10EquationOfState()) + coriolis = nothing # HydrostaticSphericalCoriolis() + free_surface = nothing # SplitExplicitFreeSurface(grid; substeps=70) + + model = HydrostaticFreeSurfaceModel(; grid, + momentum_advection, + tracer_advection, + buoyancy, + coriolis, + closure, + free_surface, + tracers = (:T, :S, :e)) + + @info "Model is built" + return model +end + +function benchmark_hydrostatic_model(Arch, grid_type, closure_type) + + grid = grids[(Arch, grid_type)] + model = ocean_benchmark(grid, closures[closure_type]) + + T = 0.0001 .* rand(size(model.grid)) .+ 20 + S = 0.0001 .* rand(size(model.grid)) .+ 35 + + set!(model.tracers.T, T) + set!(model.tracers.S, S) + + set!(model.velocities.u, (x, y, z) -> 1e-6 * rand()) + set!(model.velocities.v, (x, y, z) -> 1e-6 * rand()) + + Δt = 1 + for _ in 1:30 + time_step!(model, Δt) # warmup + end + + # Make sure we do not have any NaN or Inf values anywhere + fields = Oceananigans.fields(model) + for (key, field) in zip(propertynames(fields), fields) + arr = Array(interior(field)) + @assert all(isfinite.(arr)) || @show "Nan in $key" + @assert all(Array(arr) .< 1e10) || @show "Inf in $key" + @show key, extrema(arr) + end + + for _ in 1:10 + Oceananigans.TimeSteppers.compute_tendencies!(model, []) + end + + # NVTX.@range "compute tendencies" begin + trial = @benchmark begin + Oceananigans.TimeSteppers.compute_tendencies!($model, []) + end samples = 100 + + return trial +end + +# Problem size +Nx = 200 +Ny = 200 +Nz = 50 + +random_vector = - 5000 .* rand(Nx, Ny) + +bottom_height(arch) = GridFittedBottom(Oceananigans.on_architecture(arch, random_vector)) +lgrid(arch) = LatitudeLongitudeGrid(arch, size=(Nx, Ny, Nz), + longitude=(0, 360), + latitude=(-75, 75), + z=collect(range(-5000, 0, length=51)), + halo=(7, 7, 7)) + +grids = Dict( + (CPU, :LatitudeLongitudeGrid) => lgrid(CPU()), + (CPU, :ImmersedLatGrid) => ImmersedBoundaryGrid(lgrid(CPU()), bottom_height(CPU()); active_cells_map=true), + (GPU, :LatitudeLongitudeGrid) => lgrid(GPU()), + (GPU, :ImmersedLatGrid) => ImmersedBoundaryGrid(lgrid(GPU()), bottom_height(GPU()); active_cells_map=true), +) + +closures = Dict( + :DiffImplicit => VerticalScalarDiffusivity(TurbulenceClosures.VerticallyImplicitTimeDiscretization(), ν=1e-5 , κ=1.0), + :DiffExplicit => VerticalScalarDiffusivity(ν=1e-5, κ=1e-5), + :CATKEExplicit => Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities.CATKEVerticalDiffusivity(TurbulenceClosures.ExplicitTimeDiscretization()), + :CATKEImplicit => Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities.CATKEVerticalDiffusivity(), +) + +# Benchmark parameters + +architectures = has_cuda() ? [GPU] : [CPU] + +grid_types = [ +# :LatitudeLongitudeGrid, + :ImmersedLatGrid, +] + +closure_types = collect(keys(closures)) + +# Run and summarize benchmarks +# print_system_info() +suite = run_benchmarks(benchmark_hydrostatic_model; architectures, grid_types, closure_types) +df = benchmarks_dataframe(suite) +@show df2 = sort(df) + +for _ in 1:5 + suite = run_benchmarks(benchmark_hydrostatic_model; architectures, grid_types, closure_types) + @show df2 = sort(df) +end + +benchmarks_pretty_table(df, title="Hydrostatic model benchmarks") diff --git a/src/Advection/flux_form_advection.jl b/src/Advection/flux_form_advection.jl index c49910de5b..3bc0dfcc67 100644 --- a/src/Advection/flux_form_advection.jl +++ b/src/Advection/flux_form_advection.jl @@ -46,18 +46,18 @@ Adapt.adapt_structure(to, scheme::FluxFormAdvection{N, FT}) where {N, FT} = Adapt.adapt(to, scheme.y), Adapt.adapt(to, scheme.z)) -@inline _advective_tracer_flux_x(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_tracer_flux_x(i, j, k, grid, advection.x, args...) -@inline _advective_tracer_flux_y(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_tracer_flux_y(i, j, k, grid, advection.y, args...) -@inline _advective_tracer_flux_z(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_tracer_flux_z(i, j, k, grid, advection.z, args...) +@inline _advective_tracer_flux_x(i, j, k, grid, scheme::FluxFormAdvection, U, c) = _advective_tracer_flux_x(i, j, k, grid, scheme.x, U, c) +@inline _advective_tracer_flux_y(i, j, k, grid, scheme::FluxFormAdvection, V, c) = _advective_tracer_flux_y(i, j, k, grid, scheme.y, V, c) +@inline _advective_tracer_flux_z(i, j, k, grid, scheme::FluxFormAdvection, W, c) = _advective_tracer_flux_z(i, j, k, grid, scheme.z, W, c) -@inline _advective_momentum_flux_Uu(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Uu(i, j, k, grid, advection.x, args...) -@inline _advective_momentum_flux_Vu(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Vu(i, j, k, grid, advection.y, args...) -@inline _advective_momentum_flux_Wu(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Wu(i, j, k, grid, advection.z, args...) +@inline _advective_momentum_flux_Uu(i, j, k, grid, scheme::FluxFormAdvection, U, u) = _advective_momentum_flux_Uu(i, j, k, grid, scheme.x, U, u) +@inline _advective_momentum_flux_Vu(i, j, k, grid, scheme::FluxFormAdvection, V, u) = _advective_momentum_flux_Vu(i, j, k, grid, scheme.y, V, u) +@inline _advective_momentum_flux_Wu(i, j, k, grid, scheme::FluxFormAdvection, W, u) = _advective_momentum_flux_Wu(i, j, k, grid, scheme.z, W, u) -@inline _advective_momentum_flux_Uv(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Uv(i, j, k, grid, advection.x, args...) -@inline _advective_momentum_flux_Vv(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Vv(i, j, k, grid, advection.y, args...) -@inline _advective_momentum_flux_Wv(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Wv(i, j, k, grid, advection.z, args...) +@inline _advective_momentum_flux_Uv(i, j, k, grid, scheme::FluxFormAdvection, U, v) = _advective_momentum_flux_Uv(i, j, k, grid, scheme.x, U, v) +@inline _advective_momentum_flux_Vv(i, j, k, grid, scheme::FluxFormAdvection, V, v) = _advective_momentum_flux_Vv(i, j, k, grid, scheme.y, V, v) +@inline _advective_momentum_flux_Wv(i, j, k, grid, scheme::FluxFormAdvection, W, v) = _advective_momentum_flux_Wv(i, j, k, grid, scheme.z, W, v) -@inline _advective_momentum_flux_Uw(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Uw(i, j, k, grid, advection.x, args...) -@inline _advective_momentum_flux_Vw(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Vw(i, j, k, grid, advection.y, args...) -@inline _advective_momentum_flux_Ww(i, j, k, grid, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Ww(i, j, k, grid, advection.z, args...) +@inline _advective_momentum_flux_Uw(i, j, k, grid, scheme::FluxFormAdvection, U, w) = _advective_momentum_flux_Uw(i, j, k, grid, scheme.x, U, w) +@inline _advective_momentum_flux_Vw(i, j, k, grid, scheme::FluxFormAdvection, V, w) = _advective_momentum_flux_Vw(i, j, k, grid, scheme.y, V, w) +@inline _advective_momentum_flux_Ww(i, j, k, grid, scheme::FluxFormAdvection, W, w) = _advective_momentum_flux_Ww(i, j, k, grid, scheme.z, W, w) diff --git a/src/Advection/immersed_advective_fluxes.jl b/src/Advection/immersed_advective_fluxes.jl index abfc8a0223..cbc5da0e6b 100644 --- a/src/Advection/immersed_advective_fluxes.jl +++ b/src/Advection/immersed_advective_fluxes.jl @@ -38,25 +38,25 @@ end # dx(uu), dy(vu), dz(wu) # ccc, ffc, fcf -@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uu(i, j, k, ibg, args...)) -@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, args...) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vu(i, j, k, ibg, args...)) -@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, args...) = conditional_flux_fcf(i, j, k, ibg, zero(ibg), advective_momentum_flux_Wu(i, j, k, ibg, args...)) +@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, scheme, U, u) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uu(i, j, k, ibg, scheme, U, u)) +@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, scheme, V, u) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vu(i, j, k, ibg, scheme, V, u)) +@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, scheme, W, u) = conditional_flux_fcf(i, j, k, ibg, zero(ibg), advective_momentum_flux_Wu(i, j, k, ibg, scheme, W, u)) # dx(uv), dy(vv), dz(wv) # ffc, ccc, cff -@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, args...) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uv(i, j, k, ibg, args...)) -@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vv(i, j, k, ibg, args...)) -@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, args...) = conditional_flux_cff(i, j, k, ibg, zero(ibg), advective_momentum_flux_Wv(i, j, k, ibg, args...)) +@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, scheme, U, v) = conditional_flux_ffc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uv(i, j, k, ibg, scheme, U, v)) +@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, scheme, V, v) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vv(i, j, k, ibg, scheme, V, v)) +@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, scheme, W, v) = conditional_flux_cff(i, j, k, ibg, zero(ibg), advective_momentum_flux_Wv(i, j, k, ibg, scheme, W, v)) # dx(uw), dy(vw), dz(ww) # fcf, cff, ccc -@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, args...) = conditional_flux_fcf(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uw(i, j, k, ibg, args...)) -@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, args...) = conditional_flux_cff(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vw(i, j, k, ibg, args...)) -@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, args...) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Ww(i, j, k, ibg, args...)) +@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, scheme, U, w) = conditional_flux_fcf(i, j, k, ibg, zero(ibg), advective_momentum_flux_Uw(i, j, k, ibg, scheme, U, w)) +@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, scheme, V, w) = conditional_flux_cff(i, j, k, ibg, zero(ibg), advective_momentum_flux_Vw(i, j, k, ibg, scheme, V, w)) +@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, scheme, W, w) = conditional_flux_ccc(i, j, k, ibg, zero(ibg), advective_momentum_flux_Ww(i, j, k, ibg, scheme, W, w)) -@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, args...) = conditional_flux_fcc(i, j, k, ibg, zero(ibg), advective_tracer_flux_x(i, j, k, ibg, args...)) -@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, args...) = conditional_flux_cfc(i, j, k, ibg, zero(ibg), advective_tracer_flux_y(i, j, k, ibg, args...)) -@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, args...) = conditional_flux_ccf(i, j, k, ibg, zero(ibg), advective_tracer_flux_z(i, j, k, ibg, args...)) +@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, scheme, U, c) = conditional_flux_fcc(i, j, k, ibg, zero(ibg), advective_tracer_flux_x(i, j, k, ibg, scheme, U, c)) +@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, scheme, V, c) = conditional_flux_cfc(i, j, k, ibg, zero(ibg), advective_tracer_flux_y(i, j, k, ibg, scheme, V, c)) +@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, scheme, W, c) = conditional_flux_ccf(i, j, k, ibg, zero(ibg), advective_tracer_flux_z(i, j, k, ibg, scheme, W, c)) # Fallback for `nothing` advection @inline _advective_tracer_flux_x(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) @@ -64,27 +64,27 @@ end @inline _advective_tracer_flux_z(i, j, k, ibg::IBG, ::Nothing, args...) = zero(ibg) # Disambiguation for `FluxForm` momentum fluxes.... -@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Uu(i, j, k, ibg, advection.x, args...) -@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Vu(i, j, k, ibg, advection.y, args...) -@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Wu(i, j, k, ibg, advection.z, args...) +@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, scheme::FluxFormAdvection, U, u) = _advective_momentum_flux_Uu(i, j, k, ibg, scheme.x, U, u) +@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, scheme::FluxFormAdvection, V, u) = _advective_momentum_flux_Vu(i, j, k, ibg, scheme.y, V, u) +@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, scheme::FluxFormAdvection, W, u) = _advective_momentum_flux_Wu(i, j, k, ibg, scheme.z, W, u) -@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Uv(i, j, k, ibg, advection.x, args...) -@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Vv(i, j, k, ibg, advection.y, args...) -@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Wv(i, j, k, ibg, advection.z, args...) +@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, scheme::FluxFormAdvection, U, v) = _advective_momentum_flux_Uv(i, j, k, ibg, scheme.x, U, v) +@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, scheme::FluxFormAdvection, V, v) = _advective_momentum_flux_Vv(i, j, k, ibg, scheme.y, V, v) +@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, scheme::FluxFormAdvection, W, v) = _advective_momentum_flux_Wv(i, j, k, ibg, scheme.z, W, v) -@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Uw(i, j, k, ibg, advection.x, args...) -@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Vw(i, j, k, ibg, advection.y, args...) -@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = _advective_momentum_flux_Ww(i, j, k, ibg, advection.z, args...) +@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, scheme::FluxFormAdvection, U, w) = _advective_momentum_flux_Uw(i, j, k, ibg, scheme.x, U, w) +@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, scheme::FluxFormAdvection, V, w) = _advective_momentum_flux_Vw(i, j, k, ibg, scheme.y, V, w) +@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, scheme::FluxFormAdvection, W, w) = _advective_momentum_flux_Ww(i, j, k, ibg, scheme.z, W, w) # Disambiguation for `FluxForm` tracer fluxes.... -@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = - _advective_tracer_flux_x(i, j, k, ibg, advection.x, args...) +@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, scheme::FluxFormAdvection, U, c) = + _advective_tracer_flux_x(i, j, k, ibg, scheme.x, U, c) -@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = - _advective_tracer_flux_y(i, j, k, ibg, advection.y, args...) +@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, scheme::FluxFormAdvection, V, c) = + _advective_tracer_flux_y(i, j, k, ibg, scheme.y, V, c) -@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, advection::FluxFormAdvection, args...) = - _advective_tracer_flux_z(i, j, k, ibg, advection.z, args...) +@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, scheme::FluxFormAdvection, W, c) = + _advective_tracer_flux_z(i, j, k, ibg, scheme.z, W, c) ##### ##### "Boundary-aware" reconstruct @@ -161,13 +161,13 @@ for side in (:ᶜ, :ᶠ) near_z_boundary_bias = Symbol(:near_z_immersed_boundary_biased, side) @eval begin - @inline $near_x_boundary_symm(i, j, k, ibg, ::AbstractAdvectionScheme{0}, args...) = false - @inline $near_y_boundary_symm(i, j, k, ibg, ::AbstractAdvectionScheme{0}, args...) = false - @inline $near_z_boundary_symm(i, j, k, ibg, ::AbstractAdvectionScheme{0}, args...) = false + @inline $near_x_boundary_symm(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false + @inline $near_y_boundary_symm(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false + @inline $near_z_boundary_symm(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false - @inline $near_x_boundary_bias(i, j, k, ibg, ::AbstractAdvectionScheme{0}, args...) = false - @inline $near_y_boundary_bias(i, j, k, ibg, ::AbstractAdvectionScheme{0}, args...) = false - @inline $near_z_boundary_bias(i, j, k, ibg, ::AbstractAdvectionScheme{0}, args...) = false + @inline $near_x_boundary_bias(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false + @inline $near_y_boundary_bias(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false + @inline $near_z_boundary_bias(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false end for buffer in advection_buffers @@ -210,7 +210,7 @@ for bias in (:symmetric, :biased) @inline $alt1_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme::LOADV, args...) = $interp(i, j, k, ibg, scheme, args...) # Conditional high-order interpolation in Bounded directions - @inline $alt1_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme::HOADV, args...) = + @inline $alt1_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme::HOADV, args...) = ifelse($near_boundary(i, j, k, ibg, scheme), $alt2_interp(i, j, k, ibg, scheme.buffer_scheme, args...), $interp(i, j, k, ibg, scheme, args...)) diff --git a/src/Advection/momentum_advection_operators.jl b/src/Advection/momentum_advection_operators.jl index 2090d27144..fa334755b6 100644 --- a/src/Advection/momentum_advection_operators.jl +++ b/src/Advection/momentum_advection_operators.jl @@ -5,17 +5,17 @@ using Oceananigans.Fields: ZeroField ##### # Alternate names for advective fluxes -@inline _advective_momentum_flux_Uu(args...) = advective_momentum_flux_Uu(args...) -@inline _advective_momentum_flux_Vu(args...) = advective_momentum_flux_Vu(args...) -@inline _advective_momentum_flux_Wu(args...) = advective_momentum_flux_Wu(args...) +@inline _advective_momentum_flux_Uu(i, j, k, grid, scheme, U, u) = advective_momentum_flux_Uu(i, j, k, grid, scheme, U, u) +@inline _advective_momentum_flux_Vu(i, j, k, grid, scheme, V, u) = advective_momentum_flux_Vu(i, j, k, grid, scheme, V, u) +@inline _advective_momentum_flux_Wu(i, j, k, grid, scheme, W, u) = advective_momentum_flux_Wu(i, j, k, grid, scheme, W, u) -@inline _advective_momentum_flux_Uv(args...) = advective_momentum_flux_Uv(args...) -@inline _advective_momentum_flux_Vv(args...) = advective_momentum_flux_Vv(args...) -@inline _advective_momentum_flux_Wv(args...) = advective_momentum_flux_Wv(args...) +@inline _advective_momentum_flux_Uv(i, j, k, grid, scheme, U, v) = advective_momentum_flux_Uv(i, j, k, grid, scheme, U, v) +@inline _advective_momentum_flux_Vv(i, j, k, grid, scheme, V, v) = advective_momentum_flux_Vv(i, j, k, grid, scheme, V, v) +@inline _advective_momentum_flux_Wv(i, j, k, grid, scheme, W, v) = advective_momentum_flux_Wv(i, j, k, grid, scheme, W, v) -@inline _advective_momentum_flux_Uw(args...) = advective_momentum_flux_Uw(args...) -@inline _advective_momentum_flux_Vw(args...) = advective_momentum_flux_Vw(args...) -@inline _advective_momentum_flux_Ww(args...) = advective_momentum_flux_Ww(args...) +@inline _advective_momentum_flux_Uw(i, j, k, grid, scheme, U, w) = advective_momentum_flux_Uw(i, j, k, grid, scheme, U, w) +@inline _advective_momentum_flux_Vw(i, j, k, grid, scheme, V, w) = advective_momentum_flux_Vw(i, j, k, grid, scheme, V, w) +@inline _advective_momentum_flux_Ww(i, j, k, grid, scheme, W, w) = advective_momentum_flux_Ww(i, j, k, grid, scheme, W, w) const ZeroU = NamedTuple{(:u, :v, :w), Tuple{ZeroField, ZeroField, ZeroField}} diff --git a/src/Advection/tracer_advection_operators.jl b/src/Advection/tracer_advection_operators.jl index 667ae7ef52..051adbaf17 100644 --- a/src/Advection/tracer_advection_operators.jl +++ b/src/Advection/tracer_advection_operators.jl @@ -1,7 +1,7 @@ -@inline _advective_tracer_flux_x(args...) = advective_tracer_flux_x(args...) -@inline _advective_tracer_flux_y(args...) = advective_tracer_flux_y(args...) -@inline _advective_tracer_flux_z(args...) = advective_tracer_flux_z(args...) +@inline _advective_tracer_flux_x(i, j, k, grid, scheme, U, c) = advective_tracer_flux_x(i, j, k, grid, scheme, U, c) +@inline _advective_tracer_flux_y(i, j, k, grid, scheme, V, c) = advective_tracer_flux_y(i, j, k, grid, scheme, V, c) +@inline _advective_tracer_flux_z(i, j, k, grid, scheme, W, c) = advective_tracer_flux_z(i, j, k, grid, scheme, W, c) ##### ##### Fallback tracer fluxes! diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl index 83cd725fde..0040c66355 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl @@ -51,7 +51,8 @@ function ab2_step_velocities!(velocities, model, Δt, χ) model.closure, model.diffusivity_fields, nothing, - model.clock, + model.clock, + fields(model), Δt) end @@ -100,6 +101,7 @@ function ab2_step_tracers!(tracers, model, Δt, χ) model.diffusivity_fields, Val(tracer_index), model.clock, + fields(model), Δt) end end diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_rk3_step.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_rk3_step.jl index 11edf0f02a..767cd72b1c 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_rk3_step.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_rk3_step.jl @@ -81,7 +81,8 @@ function rk3_substep_velocities!(velocities, model, Δt, γⁿ, ζⁿ) model.closure, model.diffusivity_fields, nothing, - model.clock, + model.clock, + fields(model), Δt) end @@ -117,6 +118,7 @@ function rk3_substep_tracers!(tracers, model, Δt, γⁿ, ζⁿ) model.diffusivity_fields, Val(tracer_index), model.clock, + fields(model), Δt) end diff --git a/src/TimeSteppers/quasi_adams_bashforth_2.jl b/src/TimeSteppers/quasi_adams_bashforth_2.jl index f61fc09a0a..5f09fa5f45 100644 --- a/src/TimeSteppers/quasi_adams_bashforth_2.jl +++ b/src/TimeSteppers/quasi_adams_bashforth_2.jl @@ -147,6 +147,7 @@ function ab2_step!(model, Δt) model.diffusivity_fields, tracer_index, model.clock, + fields(model), Δt) end diff --git a/src/TimeSteppers/runge_kutta_3.jl b/src/TimeSteppers/runge_kutta_3.jl index 6c179ce250..b070667a9b 100644 --- a/src/TimeSteppers/runge_kutta_3.jl +++ b/src/TimeSteppers/runge_kutta_3.jl @@ -176,6 +176,7 @@ function rk3_substep!(model, Δt, γⁿ, ζⁿ) model.diffusivity_fields, tracer_index, model.clock, + fields(model), stage_Δt(Δt, γⁿ, ζⁿ)) end diff --git a/src/TimeSteppers/split_hydrostatic_runge_kutta_3.jl b/src/TimeSteppers/split_hydrostatic_runge_kutta_3.jl index 272b782d20..e766e52052 100644 --- a/src/TimeSteppers/split_hydrostatic_runge_kutta_3.jl +++ b/src/TimeSteppers/split_hydrostatic_runge_kutta_3.jl @@ -155,6 +155,7 @@ function split_rk3_substep!(model, Δt, γⁿ, ζⁿ) model.diffusivity_fields, tracer_index, model.clock, + fields(model), stage_Δt(Δt, γⁿ, ζⁿ)) end end diff --git a/src/TurbulenceClosures/abstract_scalar_biharmonic_diffusivity_closure.jl b/src/TurbulenceClosures/abstract_scalar_biharmonic_diffusivity_closure.jl index 11aaecc96b..b9e7cb6b92 100644 --- a/src/TurbulenceClosures/abstract_scalar_biharmonic_diffusivity_closure.jl +++ b/src/TurbulenceClosures/abstract_scalar_biharmonic_diffusivity_closure.jl @@ -60,12 +60,12 @@ const AVBD = AbstractScalarBiharmonicDiffusivity{<:VerticalFormulation} ##### Diffusive fluxes ##### -@inline diffusive_flux_x(i, j, k, grid, clo::AIBD, K, ::Val{id}, c, clk, fields, b) where id = κ_σᶠᶜᶜ(i, j, k, grid, clo, K, Val(id), clk, fields, biharmonic_mask_x, ∂xᶠᶜᶜ, ∇²ᶜᶜᶜ, c) -@inline diffusive_flux_y(i, j, k, grid, clo::AIBD, K, ::Val{id}, c, clk, fields, b) where id = κ_σᶜᶠᶜ(i, j, k, grid, clo, K, Val(id), clk, fields, biharmonic_mask_y, ∂yᶜᶠᶜ, ∇²ᶜᶜᶜ, c) -@inline diffusive_flux_z(i, j, k, grid, clo::AIBD, K, ::Val{id}, c, clk, fields, b) where id = κ_σᶜᶜᶠ(i, j, k, grid, clo, K, Val(id), clk, fields, biharmonic_mask_z, ∂zᶜᶜᶠ, ∇²ᶜᶜᶜ, c) -@inline diffusive_flux_x(i, j, k, grid, clo::AHBD, K, ::Val{id}, c, clk, fields, b) where id = κ_σᶠᶜᶜ(i, j, k, grid, clo, K, Val(id), clk, fields, biharmonic_mask_x, ∂x_∇²h_cᶠᶜᶜ, c) -@inline diffusive_flux_y(i, j, k, grid, clo::AHBD, K, ::Val{id}, c, clk, fields, b) where id = κ_σᶜᶠᶜ(i, j, k, grid, clo, K, Val(id), clk, fields, biharmonic_mask_y, ∂y_∇²h_cᶜᶠᶜ, c) -@inline diffusive_flux_z(i, j, k, grid, clo::AVBD, K, ::Val{id}, c, clk, fields, b) where id = κ_σᶜᶜᶠ(i, j, k, grid, clo, K, Val(id), clk, fields, biharmonic_mask_z, ∂³zᶜᶜᶠ, c) +@inline diffusive_flux_x(i, j, k, grid, clo::AIBD, K, id, c, clk, fields, b) = κ_σᶠᶜᶜ(i, j, k, grid, clo, K, id, clk, fields, biharmonic_mask_x, ∂xᶠᶜᶜ, ∇²ᶜᶜᶜ, c) +@inline diffusive_flux_y(i, j, k, grid, clo::AIBD, K, id, c, clk, fields, b) = κ_σᶜᶠᶜ(i, j, k, grid, clo, K, id, clk, fields, biharmonic_mask_y, ∂yᶜᶠᶜ, ∇²ᶜᶜᶜ, c) +@inline diffusive_flux_z(i, j, k, grid, clo::AIBD, K, id, c, clk, fields, b) = κ_σᶜᶜᶠ(i, j, k, grid, clo, K, id, clk, fields, biharmonic_mask_z, ∂zᶜᶜᶠ, ∇²ᶜᶜᶜ, c) +@inline diffusive_flux_x(i, j, k, grid, clo::AHBD, K, id, c, clk, fields, b) = κ_σᶠᶜᶜ(i, j, k, grid, clo, K, id, clk, fields, biharmonic_mask_x, ∂x_∇²h_cᶠᶜᶜ, c) +@inline diffusive_flux_y(i, j, k, grid, clo::AHBD, K, id, c, clk, fields, b) = κ_σᶜᶠᶜ(i, j, k, grid, clo, K, id, clk, fields, biharmonic_mask_y, ∂y_∇²h_cᶜᶠᶜ, c) +@inline diffusive_flux_z(i, j, k, grid, clo::AVBD, K, id, c, clk, fields, b) = κ_σᶜᶜᶠ(i, j, k, grid, clo, K, id, clk, fields, biharmonic_mask_z, ∂³zᶜᶜᶠ, c) ##### ##### Biharmonic-specific viscous operators diff --git a/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl b/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl index a30472c320..37ad1ece69 100644 --- a/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl +++ b/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl @@ -84,38 +84,41 @@ const AHD = AbstractScalarDiffusivity{<:Any, <:HorizontalFormulation} const ADD = AbstractScalarDiffusivity{<:Any, <:HorizontalDivergenceFormulation} const AVD = AbstractScalarDiffusivity{<:Any, <:VerticalFormulation} -@inline νᶜᶜᶜ(i, j, k, grid, closure::ASD, K, args...) = νᶜᶜᶜ(i, j, k, grid, viscosity_location(closure), viscosity(closure, K), args...) -@inline νᶠᶠᶜ(i, j, k, grid, closure::ASD, K, args...) = νᶠᶠᶜ(i, j, k, grid, viscosity_location(closure), viscosity(closure, K), args...) -@inline νᶠᶜᶠ(i, j, k, grid, closure::ASD, K, args...) = νᶠᶜᶠ(i, j, k, grid, viscosity_location(closure), viscosity(closure, K), args...) -@inline νᶜᶠᶠ(i, j, k, grid, closure::ASD, K, args...) = νᶜᶠᶠ(i, j, k, grid, viscosity_location(closure), viscosity(closure, K), args...) - -@inline κᶜᶜᶜ(i, j, k, grid, closure::ASD, K, id, args...) = κᶜᶜᶜ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), args...) -@inline κᶠᶜᶜ(i, j, k, grid, closure::ASD, K, id, args...) = κᶠᶜᶜ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), args...) -@inline κᶜᶠᶜ(i, j, k, grid, closure::ASD, K, id, args...) = κᶜᶠᶜ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), args...) -@inline κᶜᶜᶠ(i, j, k, grid, closure::ASD, K, id, args...) = κᶜᶜᶠ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), args...) +@inline νᶜᶜᶜ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶜᶜᶜ(i, j, k, grid, viscosity_location(closure), viscosity(closure, K), clk, fields) +@inline νᶠᶠᶜ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶠᶠᶜ(i, j, k, grid, viscosity_location(closure), viscosity(closure, K), clk, fields) +@inline νᶠᶜᶠ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶠᶜᶠ(i, j, k, grid, viscosity_location(closure), viscosity(closure, K), clk, fields) +@inline νᶜᶠᶠ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶜᶠᶠ(i, j, k, grid, viscosity_location(closure), viscosity(closure, K), clk, fields) + +@inline κᶜᶜᶜ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶜᶜᶜ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), clk, fields) +@inline κᶠᶜᶜ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶠᶜᶜ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), clk, fields) +@inline κᶜᶠᶜ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶜᶠᶜ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), clk, fields) +@inline κᶜᶜᶠ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶜᶜᶠ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), clk, fields) +@inline κᶠᶜᶠ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶠᶜᶠ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), clk, fields) +@inline κᶜᶠᶠ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶜᶠᶠ(i, j, k, grid, diffusivity_location(closure), diffusivity(closure, K, id), clk, fields) # Vertical and horizontal diffusivity -@inline νzᶜᶜᶜ(i, j, k, grid, closure::ASD, K, args...) = νᶜᶜᶜ(i, j, k, grid, closure, K, args...) -@inline νzᶠᶠᶜ(i, j, k, grid, closure::ASD, K, args...) = νᶠᶠᶜ(i, j, k, grid, closure, K, args...) -@inline νzᶠᶜᶠ(i, j, k, grid, closure::ASD, K, args...) = νᶠᶜᶠ(i, j, k, grid, closure, K, args...) -@inline νzᶜᶠᶠ(i, j, k, grid, closure::ASD, K, args...) = νᶜᶠᶠ(i, j, k, grid, closure, K, args...) -@inline νzᶜᶜᶜ(i, j, k, grid, closure::ASD, K, ::Nothing, args...) = νzᶜᶜᶜ(i, j, k, grid, closure, K, args...) -@inline νzᶠᶠᶜ(i, j, k, grid, closure::ASD, K, ::Nothing, args...) = νzᶠᶠᶜ(i, j, k, grid, closure, K, args...) -@inline νzᶠᶜᶠ(i, j, k, grid, closure::ASD, K, ::Nothing, args...) = νzᶠᶜᶠ(i, j, k, grid, closure, K, args...) -@inline νzᶜᶠᶠ(i, j, k, grid, closure::ASD, K, ::Nothing, args...) = νzᶜᶠᶠ(i, j, k, grid, closure, K, args...) - -@inline κzᶠᶜᶜ(i, j, k, grid, closure::ASD, K, id, args...) = κᶠᶜᶜ(i, j, k, grid, closure, K, id, args...) -@inline κzᶜᶠᶜ(i, j, k, grid, closure::ASD, K, id, args...) = κᶜᶠᶜ(i, j, k, grid, closure, K, id, args...) -@inline κzᶜᶜᶠ(i, j, k, grid, closure::ASD, K, id, args...) = κᶜᶜᶠ(i, j, k, grid, closure, K, id, args...) - -@inline νhᶜᶜᶜ(i, j, k, grid, closure::ASD, K, args...) = νᶜᶜᶜ(i, j, k, grid, closure, K, args...) -@inline νhᶠᶠᶜ(i, j, k, grid, closure::ASD, K, args...) = νᶠᶠᶜ(i, j, k, grid, closure, K, args...) -@inline νhᶠᶜᶠ(i, j, k, grid, closure::ASD, K, args...) = νᶠᶜᶠ(i, j, k, grid, closure, K, args...) -@inline νhᶜᶠᶠ(i, j, k, grid, closure::ASD, K, args...) = νᶜᶠᶠ(i, j, k, grid, closure, K, args...) - -@inline κhᶠᶜᶜ(i, j, k, grid, closure::ASD, K, id, args...) = κᶠᶜᶜ(i, j, k, grid, closure, K, id, args...) -@inline κhᶜᶠᶜ(i, j, k, grid, closure::ASD, K, id, args...) = κᶜᶠᶜ(i, j, k, grid, closure, K, id, args...) -@inline κhᶜᶜᶠ(i, j, k, grid, closure::ASD, K, id, args...) = κᶜᶜᶠ(i, j, k, grid, closure, K, id, args...) +@inline νzᶜᶜᶜ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶜᶜᶜ(i, j, k, grid, closure, K, clk, fields) +@inline νzᶠᶠᶜ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶠᶠᶜ(i, j, k, grid, closure, K, clk, fields) +@inline νzᶠᶜᶠ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶠᶜᶠ(i, j, k, grid, closure, K, clk, fields) +@inline νzᶜᶠᶠ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶜᶠᶠ(i, j, k, grid, closure, K, clk, fields) + +@inline νzᶜᶜᶜ(i, j, k, grid, closure::ASD, K, ::Nothing, clk, fields) = νzᶜᶜᶜ(i, j, k, grid, closure, K, clk, fields) +@inline νzᶠᶠᶜ(i, j, k, grid, closure::ASD, K, ::Nothing, clk, fields) = νzᶠᶠᶜ(i, j, k, grid, closure, K, clk, fields) +@inline νzᶠᶜᶠ(i, j, k, grid, closure::ASD, K, ::Nothing, clk, fields) = νzᶠᶜᶠ(i, j, k, grid, closure, K, clk, fields) +@inline νzᶜᶠᶠ(i, j, k, grid, closure::ASD, K, ::Nothing, clk, fields) = νzᶜᶠᶠ(i, j, k, grid, closure, K, clk, fields) + +@inline κzᶠᶜᶜ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶠᶜᶜ(i, j, k, grid, closure, K, id, clk, fields) +@inline κzᶜᶠᶜ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶜᶠᶜ(i, j, k, grid, closure, K, id, clk, fields) +@inline κzᶜᶜᶠ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶜᶜᶠ(i, j, k, grid, closure, K, id, clk, fields) + +@inline νhᶜᶜᶜ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶜᶜᶜ(i, j, k, grid, closure, K, clk, fields) +@inline νhᶠᶠᶜ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶠᶠᶜ(i, j, k, grid, closure, K, clk, fields) +@inline νhᶠᶜᶠ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶠᶜᶠ(i, j, k, grid, closure, K, clk, fields) +@inline νhᶜᶠᶠ(i, j, k, grid, closure::ASD, K, clk, fields) = νᶜᶠᶠ(i, j, k, grid, closure, K, clk, fields) + +@inline κhᶠᶜᶜ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶠᶜᶜ(i, j, k, grid, closure, K, id, clk, fields) +@inline κhᶜᶠᶜ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶜᶠᶜ(i, j, k, grid, closure, K, id, clk, fields) +@inline κhᶜᶜᶠ(i, j, k, grid, closure::ASD, K, id, clk, fields) = κᶜᶜᶠ(i, j, k, grid, closure, K, id, clk, fields) for (dir, Clo) in zip((:h, :z), (:AVD, :AHD)) for code in (:ᶜᶜᶜ, :ᶠᶠᶜ, :ᶠᶜᶠ, :ᶜᶠᶠ) @@ -136,25 +139,24 @@ end const F = Face const C = Center -@inline z_diffusivity(i, j, k, grid, ::F, ::C, ::C, closure::ASD, K, id, args...) = κzᶠᶜᶜ(i, j, k, grid, closure, K, id, args...) -@inline z_diffusivity(i, j, k, grid, ::C, ::F, ::C, closure::ASD, K, id, args...) = κzᶜᶠᶜ(i, j, k, grid, closure, K, id, args...) -@inline z_diffusivity(i, j, k, grid, ::C, ::C, ::F, closure::ASD, K, id, args...) = κzᶜᶜᶠ(i, j, k, grid, closure, K, id, args...) +@inline z_diffusivity(i, j, k, grid, ::F, ::C, ::C, closure::ASD, K, id, clk, fields) = κzᶠᶜᶜ(i, j, k, grid, closure, K, id, clk, fields) +@inline z_diffusivity(i, j, k, grid, ::C, ::F, ::C, closure::ASD, K, id, clk, fields) = κzᶜᶠᶜ(i, j, k, grid, closure, K, id, clk, fields) +@inline z_diffusivity(i, j, k, grid, ::C, ::C, ::F, closure::ASD, K, id, clk, fields) = κzᶜᶜᶠ(i, j, k, grid, closure, K, id, clk, fields) -@inline h_diffusivity(i, j, k, grid, ::F, ::C, ::C, closure::ASD, K, id, args...) = κhᶠᶜᶜ(i, j, k, grid, closure, K, id, args...) -@inline h_diffusivity(i, j, k, grid, ::C, ::F, ::C, closure::ASD, K, id, args...) = κhᶜᶠᶜ(i, j, k, grid, closure, K, id, args...) -@inline h_diffusivity(i, j, k, grid, ::C, ::C, ::F, closure::ASD, K, id, args...) = κhᶜᶜᶠ(i, j, k, grid, closure, K, id, args...) +@inline h_diffusivity(i, j, k, grid, ::F, ::C, ::C, closure::ASD, K, id, clk, fields) = κhᶠᶜᶜ(i, j, k, grid, closure, K, id, clk, fields) +@inline h_diffusivity(i, j, k, grid, ::C, ::F, ::C, closure::ASD, K, id, clk, fields) = κhᶜᶠᶜ(i, j, k, grid, closure, K, id, clk, fields) +@inline h_diffusivity(i, j, k, grid, ::C, ::C, ::F, closure::ASD, K, id, clk, fields) = κhᶜᶜᶠ(i, j, k, grid, closure, K, id, clk, fields) # "diffusivity" with "nothing" index => viscosity of course -@inline z_diffusivity(i, j, k, grid, ::C, ::C, ::C, closure::ASD, K, ::Nothing, args...) = νzᶜᶜᶜ(i, j, k, grid, closure, K, args...) -@inline z_diffusivity(i, j, k, grid, ::F, ::F, ::C, closure::ASD, K, ::Nothing, args...) = νzᶠᶠᶜ(i, j, k, grid, closure, K, args...) -@inline z_diffusivity(i, j, k, grid, ::F, ::C, ::F, closure::ASD, K, ::Nothing, args...) = νzᶠᶜᶠ(i, j, k, grid, closure, K, args...) -@inline z_diffusivity(i, j, k, grid, ::C, ::F, ::F, closure::ASD, K, ::Nothing, args...) = νzᶜᶠᶠ(i, j, k, grid, closure, K, args...) - -@inline h_diffusivity(i, j, k, grid, ::C, ::C, ::C, closure::ASD, K, ::Nothing, args...) = νhᶜᶜᶜ(i, j, k, grid, closure, K, args...) -@inline h_diffusivity(i, j, k, grid, ::F, ::F, ::C, closure::ASD, K, ::Nothing, args...) = νhᶠᶠᶜ(i, j, k, grid, closure, K, args...) -@inline h_diffusivity(i, j, k, grid, ::F, ::C, ::F, closure::ASD, K, ::Nothing, args...) = νhᶠᶜᶠ(i, j, k, grid, closure, K, args...) -@inline h_diffusivity(i, j, k, grid, ::C, ::F, ::F, closure::ASD, K, ::Nothing, args...) = νhᶜᶠᶠ(i, j, k, grid, closure, K, args...) +@inline z_diffusivity(i, j, k, grid, ::C, ::C, ::C, closure::ASD, K, ::Nothing, clk, fields) = νzᶜᶜᶜ(i, j, k, grid, closure, K, clk, fields) +@inline z_diffusivity(i, j, k, grid, ::F, ::F, ::C, closure::ASD, K, ::Nothing, clk, fields) = νzᶠᶠᶜ(i, j, k, grid, closure, K, clk, fields) +@inline z_diffusivity(i, j, k, grid, ::F, ::C, ::F, closure::ASD, K, ::Nothing, clk, fields) = νzᶠᶜᶠ(i, j, k, grid, closure, K, clk, fields) +@inline z_diffusivity(i, j, k, grid, ::C, ::F, ::F, closure::ASD, K, ::Nothing, clk, fields) = νzᶜᶠᶠ(i, j, k, grid, closure, K, clk, fields) +@inline h_diffusivity(i, j, k, grid, ::C, ::C, ::C, closure::ASD, K, ::Nothing, clk, fields) = νhᶜᶜᶜ(i, j, k, grid, closure, K, clk, fields) +@inline h_diffusivity(i, j, k, grid, ::F, ::F, ::C, closure::ASD, K, ::Nothing, clk, fields) = νhᶠᶠᶜ(i, j, k, grid, closure, K, clk, fields) +@inline h_diffusivity(i, j, k, grid, ::F, ::C, ::F, closure::ASD, K, ::Nothing, clk, fields) = νhᶠᶜᶠ(i, j, k, grid, closure, K, clk, fields) +@inline h_diffusivity(i, j, k, grid, ::C, ::F, ::F, closure::ASD, K, ::Nothing, clk, fields) = νhᶜᶠᶠ(i, j, k, grid, closure, K, clk, fields) # Horizontal viscous fluxes for isotropic diffusivities @inline ν_σᶜᶜᶜ(i, j, k, grid, closure, K, clock, fields, σᶜᶜᶜ, args...) = νᶜᶜᶜ(i, j, k, grid, closure, K, clock, fields) * σᶜᶜᶜ(i, j, k, grid, args...) @@ -208,9 +210,9 @@ const C = Center const AIDorAHD = Union{AID, AHD} const AIDorAVD = Union{AID, AVD} -@inline diffusive_flux_x(i, j, k, grid, cl::AIDorAHD, K, ::Val{id}, c, clk, fields, b) where id = - κhᶠᶜᶜ(i, j, k, grid, cl, K, Val(id), clk, fields) * ∂xᶠᶜᶜ(i, j, k, grid, c) -@inline diffusive_flux_y(i, j, k, grid, cl::AIDorAHD, K, ::Val{id}, c, clk, fields, b) where id = - κhᶜᶠᶜ(i, j, k, grid, cl, K, Val(id), clk, fields) * ∂yᶜᶠᶜ(i, j, k, grid, c) -@inline diffusive_flux_z(i, j, k, grid, cl::AIDorAVD, K, ::Val{id}, c, clk, fields, b) where id = - κzᶜᶜᶠ(i, j, k, grid, cl, K, Val(id), clk, fields) * ∂zᶜᶜᶠ(i, j, k, grid, c) +@inline diffusive_flux_x(i, j, k, grid, cl::AIDorAHD, K, id, c, clk, fields, b) = - κhᶠᶜᶜ(i, j, k, grid, cl, K, id, clk, fields) * ∂xᶠᶜᶜ(i, j, k, grid, c) +@inline diffusive_flux_y(i, j, k, grid, cl::AIDorAHD, K, id, c, clk, fields, b) = - κhᶜᶠᶜ(i, j, k, grid, cl, K, id, clk, fields) * ∂yᶜᶠᶜ(i, j, k, grid, c) +@inline diffusive_flux_z(i, j, k, grid, cl::AIDorAVD, K, id, c, clk, fields, b) = - κzᶜᶜᶠ(i, j, k, grid, cl, K, id, clk, fields) * ∂zᶜᶜᶠ(i, j, k, grid, c) ##### ##### Support for VerticallyImplicit @@ -224,10 +226,10 @@ const VITD = VerticallyImplicitTimeDiscretization @inline ivd_viscous_flux_vz(i, j, k, grid, closure::AVD, K, clock, fields, b) = zero(grid) # General functions (eg for vertically periodic) -@inline viscous_flux_uz(i, j, k, grid, ::VITD, closure::AIDorAVD, args...) = ivd_viscous_flux_uz(i, j, k, grid, closure, args...) -@inline viscous_flux_vz(i, j, k, grid, ::VITD, closure::AIDorAVD, args...) = ivd_viscous_flux_vz(i, j, k, grid, closure, args...) -@inline viscous_flux_wz(i, j, k, grid, ::VITD, closure::AIDorAVD, args...) = zero(grid) -@inline diffusive_flux_z(i, j, k, grid, ::VITD, closure::AIDorAVD, args...) = zero(grid) +@inline viscous_flux_uz(i, j, k, grid, ::VITD, closure::AIDorAVD, K, clock, fields, b) = ivd_viscous_flux_uz(i, j, k, grid, closure, K, clock, fields, b) +@inline viscous_flux_vz(i, j, k, grid, ::VITD, closure::AIDorAVD, K, clock, fields, b) = ivd_viscous_flux_vz(i, j, k, grid, closure, K, clock, fields, b) +@inline viscous_flux_wz(i, j, k, grid, ::VITD, closure::AIDorAVD, K, clock, fields, b) = zero(grid) +@inline diffusive_flux_z(i, j, k, grid, ::VITD, closure::AIDorAVD, K, id, c, clock, fields, b) = zero(grid) # Vertically bounded grids # @@ -237,27 +239,27 @@ const VITD = VerticallyImplicitTimeDiscretization # of boundary contributions (eg boundary contributions must be moved to the right # hand side of the tridiagonal system). -@inline function viscous_flux_uz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, args...) +@inline function viscous_flux_uz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, K, clk, fields, b) return ifelse((k == 1) | (k == grid.Nz+1), - viscous_flux_uz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...), - ivd_viscous_flux_uz(i, j, k, grid, closure, args...)) + viscous_flux_uz(i, j, k, grid, ExplicitTimeDiscretization(), closure, K, clk, fields, b), + ivd_viscous_flux_uz(i, j, k, grid, closure, K, clk, fields, b)) end -@inline function viscous_flux_vz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, args...) +@inline function viscous_flux_vz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, K, clk, fields, b) return ifelse((k == 1) | (k == grid.Nz+1), - viscous_flux_vz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...), - ivd_viscous_flux_vz(i, j, k, grid, closure, args...)) + viscous_flux_vz(i, j, k, grid, ExplicitTimeDiscretization(), closure, K, clk, fields, b), + ivd_viscous_flux_vz(i, j, k, grid, closure, K, clk, fields, b)) end -@inline function viscous_flux_wz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, args...) +@inline function viscous_flux_wz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, K, clk, fields, b) return ifelse((k == 1) | (k == grid.Nz+1), - viscous_flux_wz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...), + viscous_flux_wz(i, j, k, grid, ExplicitTimeDiscretization(), closure, K, clk, fields, b), zero(grid)) end -@inline function diffusive_flux_z(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, args...) +@inline function diffusive_flux_z(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, K, id, c, clk, fields, b) return ifelse((k == 1) | (k == grid.Nz+1), - diffusive_flux_z(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...), + diffusive_flux_z(i, j, k, grid, ExplicitTimeDiscretization(), closure, K, id, c, clk, fields, b), zero(grid)) end @@ -278,70 +280,70 @@ end # Number -@inline νᶜᶜᶜ(i, j, k, grid, loc, ν::Number, args...) = ν -@inline νᶠᶜᶠ(i, j, k, grid, loc, ν::Number, args...) = ν -@inline νᶜᶠᶠ(i, j, k, grid, loc, ν::Number, args...) = ν -@inline νᶠᶠᶜ(i, j, k, grid, loc, ν::Number, args...) = ν +@inline νᶜᶜᶜ(i, j, k, grid, loc, ν::Number, clk, fields) = ν +@inline νᶠᶜᶠ(i, j, k, grid, loc, ν::Number, clk, fields) = ν +@inline νᶜᶠᶠ(i, j, k, grid, loc, ν::Number, clk, fields) = ν +@inline νᶠᶠᶜ(i, j, k, grid, loc, ν::Number, clk, fields) = ν -@inline κᶠᶜᶜ(i, j, k, grid, loc, κ::Number, args...) = κ -@inline κᶜᶠᶜ(i, j, k, grid, loc, κ::Number, args...) = κ -@inline κᶜᶜᶠ(i, j, k, grid, loc, κ::Number, args...) = κ -@inline κᶠᶜᶠ(i, j, k, grid, loc, κ::Number, args...) = κ -@inline κᶜᶠᶠ(i, j, k, grid, loc, κ::Number, args...) = κ +@inline κᶠᶜᶜ(i, j, k, grid, loc, κ::Number, clk, fields) = κ +@inline κᶜᶠᶜ(i, j, k, grid, loc, κ::Number, clk, fields) = κ +@inline κᶜᶜᶠ(i, j, k, grid, loc, κ::Number, clk, fields) = κ +@inline κᶠᶜᶠ(i, j, k, grid, loc, κ::Number, clk, fields) = κ +@inline κᶜᶠᶠ(i, j, k, grid, loc, κ::Number, clk, fields) = κ # Array / Field at `Center, Center, Center` const Lᶜᶜᶜ = Tuple{Center, Center, Center} -@inline νᶜᶜᶜ(i, j, k, grid, ::Lᶜᶜᶜ, ν::AbstractArray, args...) = @inbounds ν[i, j, k] -@inline νᶠᶜᶠ(i, j, k, grid, ::Lᶜᶜᶜ, ν::AbstractArray, args...) = ℑxzᶠᵃᶠ(i, j, k, grid, ν) -@inline νᶜᶠᶠ(i, j, k, grid, ::Lᶜᶜᶜ, ν::AbstractArray, args...) = ℑyzᵃᶠᶠ(i, j, k, grid, ν) -@inline νᶠᶠᶜ(i, j, k, grid, ::Lᶜᶜᶜ, ν::AbstractArray, args...) = ℑxyᶠᶠᵃ(i, j, k, grid, ν) +@inline νᶜᶜᶜ(i, j, k, grid, ::Lᶜᶜᶜ, ν::AbstractArray, clk, fields) = @inbounds ν[i, j, k] +@inline νᶠᶜᶠ(i, j, k, grid, ::Lᶜᶜᶜ, ν::AbstractArray, clk, fields) = ℑxzᶠᵃᶠ(i, j, k, grid, ν) +@inline νᶜᶠᶠ(i, j, k, grid, ::Lᶜᶜᶜ, ν::AbstractArray, clk, fields) = ℑyzᵃᶠᶠ(i, j, k, grid, ν) +@inline νᶠᶠᶜ(i, j, k, grid, ::Lᶜᶜᶜ, ν::AbstractArray, clk, fields) = ℑxyᶠᶠᵃ(i, j, k, grid, ν) -@inline κᶠᶜᶜ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, args...) = ℑxᶠᵃᵃ(i, j, k, grid, κ) -@inline κᶜᶠᶜ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, args...) = ℑyᵃᶠᵃ(i, j, k, grid, κ) -@inline κᶜᶜᶠ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, args...) = ℑzᵃᵃᶠ(i, j, k, grid, κ) -@inline κᶠᶜᶠ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, args...) = ℑxzᶠᵃᶠ(i, j, k, grid, κ) -@inline κᶜᶠᶠ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, args...) = ℑyzᵃᶠᶠ(i, j, k, grid, κ) +@inline κᶠᶜᶜ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, clk, fields) = ℑxᶠᵃᵃ(i, j, k, grid, κ) +@inline κᶜᶠᶜ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, clk, fields) = ℑyᵃᶠᵃ(i, j, k, grid, κ) +@inline κᶜᶜᶠ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, clk, fields) = ℑzᵃᵃᶠ(i, j, k, grid, κ) +@inline κᶠᶜᶠ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, clk, fields) = ℑxzᶠᵃᶠ(i, j, k, grid, κ) +@inline κᶜᶠᶠ(i, j, k, grid, ::Lᶜᶜᶜ, κ::AbstractArray, clk, fields) = ℑyzᵃᶠᶠ(i, j, k, grid, κ) # Array / Field at `Center, Center, Face` const Lᶜᶜᶠ = Tuple{Center, Center, Face} -@inline νᶜᶜᶜ(i, j, k, grid, ::Lᶜᶜᶠ, ν::AbstractArray, args...) = ℑzᵃᵃᶜ(i, j, k, grid, ν) -@inline νᶠᶜᶠ(i, j, k, grid, ::Lᶜᶜᶠ, ν::AbstractArray, args...) = ℑxᶠᵃᵃ(i, j, k, grid, ν) -@inline νᶜᶠᶠ(i, j, k, grid, ::Lᶜᶜᶠ, ν::AbstractArray, args...) = ℑyᵃᶠᵃ(i, j, k, grid, ν) -@inline νᶠᶠᶜ(i, j, k, grid, ::Lᶜᶜᶠ, ν::AbstractArray, args...) = ℑxyzᶠᶠᶜ(i, j, k, grid, ν) +@inline νᶜᶜᶜ(i, j, k, grid, ::Lᶜᶜᶠ, ν::AbstractArray, clk, fields) = ℑzᵃᵃᶜ(i, j, k, grid, ν) +@inline νᶠᶜᶠ(i, j, k, grid, ::Lᶜᶜᶠ, ν::AbstractArray, clk, fields) = ℑxᶠᵃᵃ(i, j, k, grid, ν) +@inline νᶜᶠᶠ(i, j, k, grid, ::Lᶜᶜᶠ, ν::AbstractArray, clk, fields) = ℑyᵃᶠᵃ(i, j, k, grid, ν) +@inline νᶠᶠᶜ(i, j, k, grid, ::Lᶜᶜᶠ, ν::AbstractArray, clk, fields) = ℑxyzᶠᶠᶜ(i, j, k, grid, ν) -@inline κᶠᶜᶜ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, args...) = ℑxzᶠᵃᶠ(i, j, k, grid, κ) -@inline κᶜᶠᶜ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, args...) = ℑyzᵃᶠᶠ(i, j, k, grid, κ) -@inline κᶜᶜᶠ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, args...) = @inbounds κ[i, j, k] -@inline κᶠᶜᶠ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, args...) = ℑxᶠᵃᵃ(i, j, k, grid, κ) -@inline κᶜᶠᶠ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, args...) = ℑyᵃᶠᵃ(i, j, k, grid, κ) +@inline κᶠᶜᶜ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, clk, fields) = ℑxzᶠᵃᶠ(i, j, k, grid, κ) +@inline κᶜᶠᶜ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, clk, fields) = ℑyzᵃᶠᶠ(i, j, k, grid, κ) +@inline κᶜᶜᶠ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, clk, fields) = @inbounds κ[i, j, k] +@inline κᶠᶜᶠ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, clk, fields) = ℑxᶠᵃᵃ(i, j, k, grid, κ) +@inline κᶜᶠᶠ(i, j, k, grid, ::Lᶜᶜᶠ, κ::AbstractArray, clk, fields) = ℑyᵃᶠᵃ(i, j, k, grid, κ) # Function const c = Center() const f = Face() -@inline νᶜᶜᶜ(i, j, k, grid, loc, ν::F, clock, args...) where F<:Function = ν(node(i, j, k, grid, c, c, c)..., clock.time) -@inline νᶠᶜᶠ(i, j, k, grid, loc, ν::F, clock, args...) where F<:Function = ν(node(i, j, k, grid, f, c, f)..., clock.time) -@inline νᶜᶠᶠ(i, j, k, grid, loc, ν::F, clock, args...) where F<:Function = ν(node(i, j, k, grid, c, f, f)..., clock.time) -@inline νᶠᶠᶜ(i, j, k, grid, loc, ν::F, clock, args...) where F<:Function = ν(node(i, j, k, grid, f, f, c)..., clock.time) +@inline νᶜᶜᶜ(i, j, k, grid, loc, ν::Function, clk, fields) = ν(node(i, j, k, grid, c, c, c)..., clk.time) +@inline νᶠᶜᶠ(i, j, k, grid, loc, ν::Function, clk, fields) = ν(node(i, j, k, grid, f, c, f)..., clk.time) +@inline νᶜᶠᶠ(i, j, k, grid, loc, ν::Function, clk, fields) = ν(node(i, j, k, grid, c, f, f)..., clk.time) +@inline νᶠᶠᶜ(i, j, k, grid, loc, ν::Function, clk, fields) = ν(node(i, j, k, grid, f, f, c)..., clk.time) -@inline κᶜᶜᶜ(i, j, k, grid, loc, κ::F, clock, args...) where F<:Function = κ(node(i, j, k, grid, c, c, c)..., clock.time) -@inline κᶠᶜᶜ(i, j, k, grid, loc, κ::F, clock, args...) where F<:Function = κ(node(i, j, k, grid, f, c, c)..., clock.time) -@inline κᶜᶠᶜ(i, j, k, grid, loc, κ::F, clock, args...) where F<:Function = κ(node(i, j, k, grid, c, f, c)..., clock.time) -@inline κᶜᶜᶠ(i, j, k, grid, loc, κ::F, clock, args...) where F<:Function = κ(node(i, j, k, grid, c, c, f)..., clock.time) -@inline κᶠᶜᶠ(i, j, k, grid, loc, κ::F, clock, args...) where F<:Function = κ(node(i, j, k, grid, f, c, f)..., clock.time) -@inline κᶜᶠᶠ(i, j, k, grid, loc, κ::F, clock, args...) where F<:Function = κ(node(i, j, k, grid, c, f, f)..., clock.time) +@inline κᶜᶜᶜ(i, j, k, grid, loc, κ::Function, clk, fields) = κ(node(i, j, k, grid, c, c, c)..., clk.time) +@inline κᶠᶜᶜ(i, j, k, grid, loc, κ::Function, clk, fields) = κ(node(i, j, k, grid, f, c, c)..., clk.time) +@inline κᶜᶠᶜ(i, j, k, grid, loc, κ::Function, clk, fields) = κ(node(i, j, k, grid, c, f, c)..., clk.time) +@inline κᶜᶜᶠ(i, j, k, grid, loc, κ::Function, clk, fields) = κ(node(i, j, k, grid, c, c, f)..., clk.time) +@inline κᶠᶜᶠ(i, j, k, grid, loc, κ::Function, clk, fields) = κ(node(i, j, k, grid, f, c, f)..., clk.time) +@inline κᶜᶠᶠ(i, j, k, grid, loc, κ::Function, clk, fields) = κ(node(i, j, k, grid, c, f, f)..., clk.time) # "DiscreteDiffusionFunction" -@inline νᶜᶜᶜ(i, j, k, grid, loc, ν::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(ν, i, j, k, grid, (c, c, c), clock, fields) -@inline νᶠᶜᶠ(i, j, k, grid, loc, ν::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(ν, i, j, k, grid, (f, c, f), clock, fields) -@inline νᶜᶠᶠ(i, j, k, grid, loc, ν::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(ν, i, j, k, grid, (c, f, f), clock, fields) -@inline νᶠᶠᶜ(i, j, k, grid, loc, ν::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(ν, i, j, k, grid, (f, f, c), clock, fields) - -@inline κᶜᶜᶜ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(κ, i, j, k, grid, (c, c, c), clock, fields) -@inline κᶠᶜᶜ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(κ, i, j, k, grid, (f, c, c), clock, fields) -@inline κᶜᶠᶜ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(κ, i, j, k, grid, (c, f, c), clock, fields) -@inline κᶜᶜᶠ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(κ, i, j, k, grid, (c, c, f), clock, fields) -@inline κᶠᶜᶠ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(κ, i, j, k, grid, (f, c, f), clock, fields) -@inline κᶜᶠᶠ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clock, fields) = getdiffusivity(κ, i, j, k, grid, (c, f, f), clock, fields) +@inline νᶜᶜᶜ(i, j, k, grid, loc, ν::DiscreteDiffusionFunction, clk, fields) = getdiffusivity(ν, i, j, k, grid, (c, c, c), clk, fields) +@inline νᶠᶜᶠ(i, j, k, grid, loc, ν::DiscreteDiffusionFunction, clk, fields) = getdiffusivity(ν, i, j, k, grid, (f, c, f), clk, fields) +@inline νᶜᶠᶠ(i, j, k, grid, loc, ν::DiscreteDiffusionFunction, clk, fields) = getdiffusivity(ν, i, j, k, grid, (c, f, f), clk, fields) +@inline νᶠᶠᶜ(i, j, k, grid, loc, ν::DiscreteDiffusionFunction, clk, fields) = getdiffusivity(ν, i, j, k, grid, (f, f, c), clk, fields) + +@inline κᶜᶜᶜ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clk, fields) = getdiffusivity(κ, i, j, k, grid, (c, c, c), clk, fields) +@inline κᶠᶜᶜ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clk, fields) = getdiffusivity(κ, i, j, k, grid, (f, c, c), clk, fields) +@inline κᶜᶠᶜ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clk, fields) = getdiffusivity(κ, i, j, k, grid, (c, f, c), clk, fields) +@inline κᶜᶜᶠ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clk, fields) = getdiffusivity(κ, i, j, k, grid, (c, c, f), clk, fields) +@inline κᶠᶜᶠ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clk, fields) = getdiffusivity(κ, i, j, k, grid, (f, c, f), clk, fields) +@inline κᶜᶠᶠ(i, j, k, grid, loc, κ::DiscreteDiffusionFunction, clk, fields) = getdiffusivity(κ, i, j, k, grid, (c, f, f), clk, fields) diff --git a/src/TurbulenceClosures/closure_kernel_operators.jl b/src/TurbulenceClosures/closure_kernel_operators.jl index 8691a2a223..f60af9d2ba 100644 --- a/src/TurbulenceClosures/closure_kernel_operators.jl +++ b/src/TurbulenceClosures/closure_kernel_operators.jl @@ -19,30 +19,30 @@ using Oceananigans.Operators: Δy_qᶠᶜᶜ, Δx_qᶜᶠᶜ, Δx_qᶠᶜᶜ ##### Viscous flux divergences ##### -@inline function ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure::AbstractTurbulenceClosure, args...) +@inline function ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure::AbstractTurbulenceClosure, K, clock, fields, buoyancy) disc = time_discretization(closure) - return V⁻¹ᶠᶜᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Ax_qᶜᶜᶜ, _viscous_flux_ux, disc, closure, args...) + - δyᵃᶜᵃ(i, j, k, grid, Ay_qᶠᶠᶜ, _viscous_flux_uy, disc, closure, args...) + - δzᵃᵃᶜ(i, j, k, grid, Az_qᶠᶜᶠ, _viscous_flux_uz, disc, closure, args...)) + return V⁻¹ᶠᶜᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Ax_qᶜᶜᶜ, _viscous_flux_ux, disc, closure, K, clock, fields, buoyancy) + + δyᵃᶜᵃ(i, j, k, grid, Ay_qᶠᶠᶜ, _viscous_flux_uy, disc, closure, K, clock, fields, buoyancy) + + δzᵃᵃᶜ(i, j, k, grid, Az_qᶠᶜᶠ, _viscous_flux_uz, disc, closure, K, clock, fields, buoyancy)) end -@inline function ∂ⱼ_τ₂ⱼ(i, j, k, grid, closure::AbstractTurbulenceClosure, args...) +@inline function ∂ⱼ_τ₂ⱼ(i, j, k, grid, closure::AbstractTurbulenceClosure, K, clock, fields, buoyancy) disc = time_discretization(closure) - return V⁻¹ᶜᶠᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶠᶜ, _viscous_flux_vx, disc, closure, args...) + - δyᵃᶠᵃ(i, j, k, grid, Ay_qᶜᶜᶜ, _viscous_flux_vy, disc, closure, args...) + - δzᵃᵃᶜ(i, j, k, grid, Az_qᶜᶠᶠ, _viscous_flux_vz, disc, closure, args...)) + return V⁻¹ᶜᶠᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶠᶜ, _viscous_flux_vx, disc, closure, K, clock, fields, buoyancy) + + δyᵃᶠᵃ(i, j, k, grid, Ay_qᶜᶜᶜ, _viscous_flux_vy, disc, closure, K, clock, fields, buoyancy) + + δzᵃᵃᶜ(i, j, k, grid, Az_qᶜᶠᶠ, _viscous_flux_vz, disc, closure, K, clock, fields, buoyancy)) end -@inline function ∂ⱼ_τ₃ⱼ(i, j, k, grid, closure::AbstractTurbulenceClosure, args...) +@inline function ∂ⱼ_τ₃ⱼ(i, j, k, grid, closure::AbstractTurbulenceClosure, K, clock, fields, buoyancy) disc = time_discretization(closure) - return V⁻¹ᶜᶜᶠ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶜᶠ, _viscous_flux_wx, disc, closure, args...) + - δyᵃᶜᵃ(i, j, k, grid, Ay_qᶜᶠᶠ, _viscous_flux_wy, disc, closure, args...) + - δzᵃᵃᶠ(i, j, k, grid, Az_qᶜᶜᶜ, _viscous_flux_wz, disc, closure, args...)) + return V⁻¹ᶜᶜᶠ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶜᶠ, _viscous_flux_wx, disc, closure, K, clock, fields, buoyancy) + + δyᵃᶜᵃ(i, j, k, grid, Ay_qᶜᶠᶠ, _viscous_flux_wy, disc, closure, K, clock, fields, buoyancy) + + δzᵃᵃᶠ(i, j, k, grid, Az_qᶜᶜᶜ, _viscous_flux_wz, disc, closure, K, clock, fields, buoyancy)) end -@inline function ∇_dot_qᶜ(i, j, k, grid, closure::AbstractTurbulenceClosure, diffusivities, tracer_index, args...) +@inline function ∇_dot_qᶜ(i, j, k, grid, closure::AbstractTurbulenceClosure, K, id, c, clock, fields, buoyancy) disc = time_discretization(closure) - return V⁻¹ᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶜᶜ, _diffusive_flux_x, disc, closure, diffusivities, tracer_index, args...) + - δyᵃᶜᵃ(i, j, k, grid, Ay_qᶜᶠᶜ, _diffusive_flux_y, disc, closure, diffusivities, tracer_index, args...) + - δzᵃᵃᶜ(i, j, k, grid, Az_qᶜᶜᶠ, _diffusive_flux_z, disc, closure, diffusivities, tracer_index, args...)) + return V⁻¹ᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶜᶜ, _diffusive_flux_x, disc, closure, K, id, c, clock, fields, buoyancy) + + δyᵃᶜᵃ(i, j, k, grid, Ay_qᶜᶠᶜ, _diffusive_flux_y, disc, closure, K, id, c, clock, fields, buoyancy) + + δzᵃᵃᶜ(i, j, k, grid, Az_qᶜᶜᶠ, _diffusive_flux_z, disc, closure, K, id, c, clock, fields, buoyancy)) end diff --git a/src/TurbulenceClosures/immersed_diffusive_fluxes.jl b/src/TurbulenceClosures/immersed_diffusive_fluxes.jl index 5bff5ebba8..f58dc14efb 100644 --- a/src/TurbulenceClosures/immersed_diffusive_fluxes.jl +++ b/src/TurbulenceClosures/immersed_diffusive_fluxes.jl @@ -105,42 +105,42 @@ end @inline function _west_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields) Δ = Δx(index_left(i, LX), j, k, ibg, LX, LY, LZ) - κ = h_diffusivity(i, j, k, ibg, flip(LX), LY, LZ, closure, K, id, clock) + κ = h_diffusivity(i, j, k, ibg, flip(LX), LY, LZ, closure, K, id, clock, fields) ∇c = left_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields) return - κ * ∇c end @inline function _east_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields) Δ = Δx(index_right(i, LX), j, k, ibg, LX, LY, LZ) - κ = h_diffusivity(i, j, k, ibg, flip(LX), LY, LZ, closure, K, id, clock) + κ = h_diffusivity(i, j, k, ibg, flip(LX), LY, LZ, closure, K, id, clock, fields) ∇c = right_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields) return - κ * ∇c end @inline function _south_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields) Δ = Δy(i, index_left(j, LY), k, ibg, LX, LY, LZ) - κ = h_diffusivity(i, j, k, ibg, LX, flip(LY), LZ, closure, K, id, clock) + κ = h_diffusivity(i, j, k, ibg, LX, flip(LY), LZ, closure, K, id, clock, fields) ∇c = left_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields) return - κ * ∇c end @inline function _north_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields) Δ = Δy(i, index_right(j, LY), k, ibg, LX, LY, LZ) - κ = h_diffusivity(i, j, k, ibg, LX, flip(LY), LZ, closure, K, id, clock) + κ = h_diffusivity(i, j, k, ibg, LX, flip(LY), LZ, closure, K, id, clock, fields) ∇c = right_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields) return - κ * ∇c end @inline function _bottom_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields) Δ = Δz(i, j, index_left(k, LZ), ibg, LX, LY, LZ) - κ = z_diffusivity(i, j, k, ibg, LX, LY, flip(LZ), closure, K, id, clock) + κ = z_diffusivity(i, j, k, ibg, LX, LY, flip(LZ), closure, K, id, clock, fields) ∇c = left_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields) return - κ * ∇c end @inline function _top_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields) Δ = Δz(i, j, index_right(k, LZ), ibg, LX, LY, LZ) - κ = z_diffusivity(i, j, k, ibg, LX, LY, flip(LZ), closure, K, id, clock) + κ = z_diffusivity(i, j, k, ibg, LX, LY, flip(LZ), closure, K, id, clock, fields) ∇c = right_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields) return - κ * ∇c end diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl index 585f1c2036..790c165ab5 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl @@ -74,7 +74,9 @@ function time_step_catke_equation!(model) implicit_step!(e, implicit_solver, closure, diffusivity_fields, Val(tracer_index), - model.clock, Δτ) + model.clock, + fields(model), + Δτ) end return nothing @@ -152,6 +154,8 @@ end # = Lⁱ # # where ω = ϵ / e ∼ √e / ℓ. + + active = !inactive_cell(i, j, k, grid) @inbounds Le[i, j, k] = (wb⁻_e - ω + div_Jᵉ_e) * active @@ -175,7 +179,7 @@ end # See below. α = convert(FT, 1.5) + χ β = convert(FT, 0.5) + χ - + @inbounds begin total_Gⁿe = slow_Gⁿe[i, j, k] + fast_Gⁿe e[i, j, k] += Δτ * (α * total_Gⁿe - β * G⁻e[i, j, k]) * active diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl index c838820825..0653fab32d 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl @@ -83,11 +83,15 @@ function time_step_tke_dissipation_equations!(model) implicit_step!(e, implicit_solver, closure, model.diffusivity_fields, Val(e_index), - model.clock, Δτ) + model.clock, + fields(model), + Δτ) implicit_step!(ϵ, implicit_solver, closure, model.diffusivity_fields, Val(ϵ_index), - model.clock, Δτ) + model.clock, + fields(model), + Δτ) end return nothing diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl index b6e900e496..e946a89491 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl @@ -235,8 +235,8 @@ end κ_skew = get_tracer_κ(closure.κ_skew, grid, tracer_index) κ_symmetric = get_tracer_κ(closure.κ_symmetric, grid, tracer_index) - κ_skewᶠᶜᶜ = skew_diffusivity(i, j, k, grid, closure, κᶠᶜᶜ, issd_coefficient_loc, κ_skew, clock) - κ_symmetricᶠᶜᶜ = κᶠᶜᶜ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) + κ_skewᶠᶜᶜ = skew_diffusivity(i, j, k, grid, closure, κᶠᶜᶜ, issd_coefficient_loc, κ_skew, clock, fields) + κ_symmetricᶠᶜᶜ = κᶠᶜᶜ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock, fields) ∂x_c = ∂xᶠᶜᶜ(i, j, k, grid, c) @@ -265,8 +265,8 @@ end κ_skew = get_tracer_κ(closure.κ_skew, grid, tracer_index) κ_symmetric = get_tracer_κ(closure.κ_symmetric, grid, tracer_index) - κ_skewᶜᶠᶜ = skew_diffusivity(i, j, k, grid, closure, κᶜᶠᶜ, issd_coefficient_loc, κ_skew, clock) - κ_symmetricᶜᶠᶜ = κᶜᶠᶜ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) + κ_skewᶜᶠᶜ = skew_diffusivity(i, j, k, grid, closure, κᶜᶠᶜ, issd_coefficient_loc, κ_skew, clock, fields) + κ_symmetricᶜᶠᶜ = κᶜᶠᶜ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock, fields) ∂y_c = ∂yᶜᶠᶜ(i, j, k, grid, c) @@ -295,8 +295,8 @@ end κ_skew = get_tracer_κ(closure.κ_skew, grid, tracer_index) κ_symmetric = get_tracer_κ(closure.κ_symmetric, grid, tracer_index) - κ_skewᶜᶜᶠ = skew_diffusivity(i, j, k, grid, closure, κᶜᶜᶠ, issd_coefficient_loc, κ_skew, clock) - κ_symmetricᶜᶜᶠ = κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) + κ_skewᶜᶜᶠ = skew_diffusivity(i, j, k, grid, closure, κᶜᶜᶠ, issd_coefficient_loc, κ_skew, clock, fields) + κ_symmetricᶜᶜᶠ = κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock, fields) # Average... of... the gradient! ∂x_c = ℑxzᶜᵃᶠ(i, j, k, grid, ∂xᶠᶜᶜ, c) @@ -324,11 +324,11 @@ end @inline explicit_κ_∂z_c(i, j, k, grid, ::VerticallyImplicitTimeDiscretization, args...) = zero(grid) -@inline function κzᶜᶜᶠ(i, j, k, grid, closure::FlavorOfISSD, K, ::Val{id}, clock) where id +@inline function κzᶜᶜᶠ(i, j, k, grid, closure::FlavorOfISSD, K, ::Val{id}, clock, fields) where id closure = getclosure(i, j, closure) κ_symmetric = get_tracer_κ(closure.κ_symmetric, grid, id) ϵ_R₃₃ = @inbounds K.ϵ_R₃₃[i, j, k] # tapered 33 component of rotation tensor - return ϵ_R₃₃ * κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock) + return ϵ_R₃₃ * κᶜᶜᶠ(i, j, k, grid, issd_coefficient_loc, κ_symmetric, clock, fields) end @inline viscous_flux_ux(i, j, k, grid, closure::Union{ISSD, ISSDVector}, args...) = zero(grid) diff --git a/src/TurbulenceClosures/vertically_implicit_diffusion_solver.jl b/src/TurbulenceClosures/vertically_implicit_diffusion_solver.jl index dbf1e355d6..4219546f79 100644 --- a/src/TurbulenceClosures/vertically_implicit_diffusion_solver.jl +++ b/src/TurbulenceClosures/vertically_implicit_diffusion_solver.jl @@ -30,18 +30,18 @@ const F = Face @inline implicit_linear_coefficient(i, j, k, grid, args...) = zero(grid) # General implementation -@inline νzᶠᶜᶠ(i, j, k, grid, closure, K, args...) = zero(grid) -@inline νzᶜᶠᶠ(i, j, k, grid, closure, K, args...) = zero(grid) -@inline νzᶜᶜᶜ(i, j, k, grid, closure, K, args...) = zero(grid) -@inline κzᶜᶜᶠ(i, j, k, grid, closure, K, args...) = zero(grid) +@inline νzᶠᶜᶠ(i, j, k, grid, args...) = zero(grid) +@inline νzᶜᶠᶠ(i, j, k, grid, args...) = zero(grid) +@inline νzᶜᶜᶜ(i, j, k, grid, args...) = zero(grid) +@inline κzᶜᶜᶠ(i, j, k, grid, args...) = zero(grid) # Vertical momentum diffusivities: u, v, w -@inline ivd_diffusivity(i, j, k, grid, ::F, ::C, ::F, clo, K, id, clock) = νzᶠᶜᶠ(i, j, k, grid, clo, K, id, clock) * !inactive_node(i, j, k, grid, f, c, f) -@inline ivd_diffusivity(i, j, k, grid, ::C, ::F, ::F, clo, K, id, clock) = νzᶜᶠᶠ(i, j, k, grid, clo, K, id, clock) * !inactive_node(i, j, k, grid, c, f, f) -@inline ivd_diffusivity(i, j, k, grid, ::C, ::C, ::C, clo, K, id, clock) = νzᶜᶜᶜ(i, j, k, grid, clo, K, id, clock) * !inactive_node(i, j, k, grid, c, c, c) +@inline ivd_diffusivity(i, j, k, grid, ::F, ::C, ::F, clo, K, id, clk, fields) = νzᶠᶜᶠ(i, j, k, grid, clo, K, id, clk, fields) * !inactive_node(i, j, k, grid, f, c, f) +@inline ivd_diffusivity(i, j, k, grid, ::C, ::F, ::F, clo, K, id, clk, fields) = νzᶜᶠᶠ(i, j, k, grid, clo, K, id, clk, fields) * !inactive_node(i, j, k, grid, c, f, f) +@inline ivd_diffusivity(i, j, k, grid, ::C, ::C, ::C, clo, K, id, clk, fields) = νzᶜᶜᶜ(i, j, k, grid, clo, K, id, clk, fields) * !inactive_node(i, j, k, grid, c, c, c) # Tracer diffusivity -@inline ivd_diffusivity(i, j, k, grid, ::C, ::C, ::F, args...) = κzᶜᶜᶠ(i, j, k, grid, args...) * !inactive_node(i, j, k, grid, c, c, f) +@inline ivd_diffusivity(i, j, k, grid, ::C, ::C, ::F, clo, K, id, clk, fields) = κzᶜᶜᶠ(i, j, k, grid, clo, K, id, clk, fields) * !inactive_node(i, j, k, grid, c, c, f) ##### ##### Batched Tridiagonal solver for implicit diffusion @@ -61,9 +61,9 @@ implicit_diffusion_solver(::ExplicitTimeDiscretization, args...; kwargs...) = no @inline rcp_vertical_spacing(i, j, k, grid, ::Center, ::Center, ℓz) = Δr⁻¹(i, j, k, grid, c, c, ℓz) # Tracers and horizontal velocities at cell centers in z -@inline function ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ::Center, Δt, clock) +@inline function ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ::Center, Δt, clock, fields) closure_ij = getclosure(i, j, closure) - κᵏ⁺¹ = ivd_diffusivity(i, j, k+1, grid, ℓx, ℓy, f, closure_ij, K, id, clock) + κᵏ⁺¹ = ivd_diffusivity(i, j, k+1, grid, ℓx, ℓy, f, closure_ij, K, id, clock, fields) Δz⁻¹ᶜₖ = rcp_vertical_spacing(i, j, k, grid, ℓx, ℓy, c) Δz⁻¹ᶠₖ₊₁ = rcp_vertical_spacing(i, j, k+1, grid, ℓx, ℓy, f) du = - Δt * κᵏ⁺¹ * (Δz⁻¹ᶜₖ * Δz⁻¹ᶠₖ₊₁) @@ -71,10 +71,10 @@ implicit_diffusion_solver(::ExplicitTimeDiscretization, args...; kwargs...) = no return du * !peripheral_node(i, j, k+1, grid, ℓx, ℓy, f) end -@inline function ivd_lower_diagonal(i, j, k′, grid, closure, K, id, ℓx, ℓy, ::Center, Δt, clock) +@inline function ivd_lower_diagonal(i, j, k′, grid, closure, K, id, ℓx, ℓy, ::Center, Δt, clock, fields) k = k′ + 1 # Shift index to match LinearAlgebra.Tridiagonal indexing convenction closure_ij = getclosure(i, j, closure) - κᵏ = ivd_diffusivity(i, j, k, grid, ℓx, ℓy, f, closure_ij, K, id, clock) + κᵏ = ivd_diffusivity(i, j, k, grid, ℓx, ℓy, f, closure_ij, K, id, clock, fields) Δz⁻¹ᶜₖ = rcp_vertical_spacing(i, j, k, grid, ℓx, ℓy, c) Δz⁻¹ᶠₖ = rcp_vertical_spacing(i, j, k, grid, ℓx, ℓy, f) dl = - Δt * κᵏ * (Δz⁻¹ᶜₖ * Δz⁻¹ᶠₖ) @@ -90,19 +90,19 @@ end ##### Note: these coefficients are specific to vertically-bounded grids (and so is ##### the BatchedTridiagonalSolver). -@inline function ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ::Face, Δt, clock) +@inline function ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ::Face, Δt, clock, fields) closure_ij = getclosure(i, j, closure) - νᵏ = ivd_diffusivity(i, j, k, grid, ℓx, ℓy, c, closure_ij, K, id, clock) + νᵏ = ivd_diffusivity(i, j, k, grid, ℓx, ℓy, c, closure_ij, K, id, clock, fields) Δz⁻¹ᶜₖ = rcp_vertical_spacing(i, j, k, grid, ℓx, ℓy, c) Δz⁻¹ᶠₖ = rcp_vertical_spacing(i, j, k, grid, ℓx, ℓy, f) du = - Δt * νᵏ * (Δz⁻¹ᶜₖ * Δz⁻¹ᶠₖ) return du * !peripheral_node(i, j, k, grid, ℓx, ℓy, c) end -@inline function ivd_lower_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ::Face, Δt, clock) +@inline function ivd_lower_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ::Face, Δt, clock, fields) k′ = k + 2 # Shift to adjust for Tridiagonal indexing convention closure_ij = getclosure(i, j, closure) - νᵏ⁻¹ = ivd_diffusivity(i, j, k′-1, grid, ℓx, ℓy, c, closure_ij, K, id, clock) + νᵏ⁻¹ = ivd_diffusivity(i, j, k′-1, grid, ℓx, ℓy, c, closure_ij, K, id, clock, fields) Δz⁻¹ᶜₖ = rcp_vertical_spacing(i, j, k′, grid, ℓx, ℓy, c) Δz⁻¹ᶠₖ₋₁ = rcp_vertical_spacing(i, j, k′-1, grid, ℓx, ℓy, f) dl = - Δt * νᵏ⁻¹ * (Δz⁻¹ᶜₖ * Δz⁻¹ᶠₖ₋₁) @@ -110,16 +110,21 @@ end end ### Diagonal terms +@inline ivd_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) = + one(grid) - Δt * _implicit_linear_coefficient(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) - + _ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) - + _ivd_lower_diagonal(i, j, k-1, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) -@inline ivd_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock) = - one(grid) - Δt * _implicit_linear_coefficient(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock) - - _ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock) - - _ivd_lower_diagonal(i, j, k-1, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock) +# Fallback for single closure. These coefficients are extended for tupled closures in `closure_tuples.jl` +@inline _implicit_linear_coefficient(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) = + implicit_linear_coefficient(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) -@inline _implicit_linear_coefficient(args...) = implicit_linear_coefficient(args...) -@inline _ivd_upper_diagonal(args...) = ivd_upper_diagonal(args...) -@inline _ivd_lower_diagonal(args...) = ivd_lower_diagonal(args...) +@inline _ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) = + ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) + +@inline _ivd_lower_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) = + ivd_lower_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ℓz, Δt, clock, fields) ##### ##### Solver constructor @@ -163,9 +168,14 @@ end # Extend `get_coefficient` to retrieve `ivd_diagonal`, `_ivd_lower_diagonal` and `_ivd_upper_diagonal`. # Note that we use the "periphery-aware" upper and lower diagonals -@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionLowerDiagonal, p, ::ZDirection, args...) = _ivd_lower_diagonal(i, j, k, grid, args...) -@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionUpperDiagonal, p, ::ZDirection, args...) = _ivd_upper_diagonal(i, j, k, grid, args...) -@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionDiagonal, p, ::ZDirection, args...) = ivd_diagonal(i, j, k, grid, args...) +@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionLowerDiagonal, p, ::ZDirection, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) = + _ivd_lower_diagonal(i, j, k, grid, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) + +@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionUpperDiagonal, p, ::ZDirection, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) = + _ivd_upper_diagonal(i, j, k, grid, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) + +@inline get_coefficient(i, j, k, grid, ::VerticallyImplicitDiffusionDiagonal, p, ::ZDirection, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) = + ivd_diagonal(i, j, k, grid, clo, K, id, ℓx, ℓy, ℓz, Δt, clk, fields) ##### ##### Implicit step functions @@ -188,9 +198,9 @@ function implicit_step!(field::Field, diffusivity_fields, tracer_index, clock, - Δt; - kwargs...) - + fields, + Δt) + # Filter explicit closures for closure tuples if closure isa Tuple closure_tuple = closure @@ -203,9 +213,7 @@ function implicit_step!(field::Field, end LX, LY, LZ = location(field) - # Nullify tracer_index if `field` is not a tracer - (LX, LY, LZ) == (Center, Center, Center) || (tracer_index = nothing) return solve!(field, implicit_solver, field, # ivd_*_diagonal gets called with these args after (i, j, k, grid): - vi_closure, vi_diffusivity_fields, tracer_index, LX(), LY(), LZ(), Δt, clock; kwargs...) + vi_closure, vi_diffusivity_fields, tracer_index, LX(), LY(), LZ(), Δt, clock, fields) end