diff --git a/.gitignore b/.gitignore index 2680fc9..cccfb88 100644 --- a/.gitignore +++ b/.gitignore @@ -32,3 +32,8 @@ docs/Manifest.toml *.json *.sha512 data + +# File generated by NVIDIA profiling tools +*.qdprep +*.nsys-rep +*.ncu-rep \ No newline at end of file diff --git a/Manifest.toml b/Manifest.toml index bd9a9c9..cc3de0d 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,8 +1,8 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.10.2" +julia_version = "1.10.4" manifest_format = "2.0" -project_hash = "c7a99ea0ec59ff7a94dd62ff00b2ed61ec5aa58c" +project_hash = "19484cf8dfad99c0f33810f2fcf7362f476ab5bf" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] @@ -52,14 +52,15 @@ version = "1.1.1" [[deps.ArrayInterface]] deps = ["Adapt", "LinearAlgebra", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "44691067188f6bd1b2289552a23e4b7572f4528d" +git-tree-sha1 = "5c9b74c973181571deb6442d41e5c902e6b9f38e" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "7.9.0" +version = "7.12.0" [deps.ArrayInterface.extensions] ArrayInterfaceBandedMatricesExt = "BandedMatrices" ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices" ArrayInterfaceCUDAExt = "CUDA" + ArrayInterfaceCUDSSExt = "CUDSS" ArrayInterfaceChainRulesExt = "ChainRules" ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore" ArrayInterfaceReverseDiffExt = "ReverseDiff" @@ -70,6 +71,7 @@ version = "7.9.0" BandedMatrices = "aae01518-5342-5314-be14-df237901396f" BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e" ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2" GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" @@ -119,38 +121,43 @@ version = "0.1.3" [[deps.CUDA]] deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "DataFrames", "ExprTools", "GPUArrays", "GPUCompiler", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "StaticArrays", "Statistics"] -git-tree-sha1 = "baa8ea7a1ea63316fa3feb454635215773c9c845" +git-tree-sha1 = "fdd9dfb67dfefd548f51000cc400bb51003de247" uuid = "052768ef-5323-5732-b1bb-66c8b64840ba" -version = "5.2.0" -weakdeps = ["ChainRulesCore", "SpecialFunctions"] +version = "5.4.3" [deps.CUDA.extensions] ChainRulesCoreExt = "ChainRulesCore" + EnzymeCoreExt = "EnzymeCore" SpecialFunctionsExt = "SpecialFunctions" + [deps.CUDA.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" + SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" + [[deps.CUDA_Driver_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] -git-tree-sha1 = "d01bfc999768f0a31ed36f5d22a76161fc63079c" +git-tree-sha1 = "97df9d4d6be8ac6270cb8fd3b8fc413690820cbd" uuid = "4ee394cb-3365-5eb0-8335-949819d2adfc" -version = "0.7.0+1" +version = "0.9.1+1" [[deps.CUDA_Runtime_Discovery]] deps = ["Libdl"] -git-tree-sha1 = "2cb12f6b2209f40a4b8967697689a47c50485490" +git-tree-sha1 = "f3b237289a5a77c759b2dd5d4c2ff641d67c4030" uuid = "1af6417a-86b4-443c-805f-a4643ffb695f" -version = "0.2.3" +version = "0.3.4" [[deps.CUDA_Runtime_jll]] deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] -git-tree-sha1 = "8e25c009d2bf16c2c31a70a6e9e8939f7325cc84" +git-tree-sha1 = "afea94249b821dc754a8ca6695d3daed851e1f5a" uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2" -version = "0.11.1+0" +version = "0.14.1+0" [[deps.ChainRulesCore]] deps = ["Compat", "LinearAlgebra"] -git-tree-sha1 = "575cd02e080939a33b6df6c5853d14924c08e35b" +git-tree-sha1 = "71acdbf594aab5bbb2cec89b208c41b4c411e49f" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.23.0" +version = "1.24.0" weakdeps = ["SparseArrays"] [deps.ChainRulesCore.extensions] @@ -158,21 +165,21 @@ weakdeps = ["SparseArrays"] [[deps.ColorTypes]] deps = ["FixedPointNumbers", "Random"] -git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4" +git-tree-sha1 = "b10d0b65641d57b8b4d5e234446582de5047050d" uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" -version = "0.11.4" +version = "0.11.5" [[deps.Colors]] deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] -git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a" +git-tree-sha1 = "362a287c3aa50601b0bc359053d5c2468f0e7ce0" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" -version = "0.12.10" +version = "0.12.11" [[deps.CommonDataModel]] deps = ["CFTime", "DataStructures", "Dates", "Preferences", "Printf", "Statistics"] -git-tree-sha1 = "d7d7b58e149f19c322840a50d1bc20e8c23addb4" +git-tree-sha1 = "d6fb5bf939a2753c74984b11434ea25d6c397a58" uuid = "1fbeeb36-5f17-413c-809b-666fb144f157" -version = "0.3.5" +version = "0.3.6" [[deps.CommonSolve]] git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c" @@ -185,11 +192,16 @@ git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" version = "0.3.0" +[[deps.CommonWorldInvalidations]] +git-tree-sha1 = "ae52d1c52048455e85a387fbee9be553ec2b68d0" +uuid = "f70d9fcc-98c5-4d4a-abd7-e4cdeebd8ca8" +version = "1.0.0" + [[deps.Compat]] deps = ["TOML", "UUIDs"] -git-tree-sha1 = "c955881e3c981181362ae4088b35995446298b80" +git-tree-sha1 = "b1c55339b7c6c350ee89f2c1604299660525b248" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.14.0" +version = "4.15.0" weakdeps = ["Dates", "LinearAlgebra"] [deps.Compat.extensions] @@ -198,7 +210,7 @@ weakdeps = ["Dates", "LinearAlgebra"] [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.1.0+0" +version = "1.1.1+0" [[deps.CompositionsBase]] git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad" @@ -247,9 +259,9 @@ version = "1.6.1" [[deps.DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "0f4b5d62a88d8f59003e43c25a8a90de9eb76317" +git-tree-sha1 = "1d0a14036acb104d9e89698bd408f63ab58cdc82" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.18" +version = "0.18.20" [[deps.DataValueInterfaces]] git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" @@ -337,9 +349,9 @@ uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" [[deps.FixedPointNumbers]] deps = ["Statistics"] -git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc" +git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172" uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" -version = "0.8.4" +version = "0.8.5" [[deps.ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"] @@ -362,9 +374,9 @@ version = "6.2.1+6" [[deps.GPUArrays]] deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"] -git-tree-sha1 = "47e4686ec18a9620850bad110b79966132f14283" +git-tree-sha1 = "04661708f5301394a1f1be86a07a89e835900db6" uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" -version = "10.0.2" +version = "10.2.3" [[deps.GPUArraysCore]] deps = ["Adapt"] @@ -373,10 +385,10 @@ uuid = "46192b85-c4d5-4398-a991-12ede77f4527" version = "0.1.6" [[deps.GPUCompiler]] -deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "TimerOutputs", "UUIDs"] -git-tree-sha1 = "a846f297ce9d09ccba02ead0cae70690e072a119" +deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Preferences", "Scratch", "Serialization", "TOML", "TimerOutputs", "UUIDs"] +git-tree-sha1 = "ab29216184312f99ff957b32cd63c2fe9c928b91" uuid = "61eb1bfa-7361-4325-ad38-22787b887f55" -version = "0.25.0" +version = "0.26.7" [[deps.Glob]] git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496" @@ -397,9 +409,9 @@ version = "1.14.2+1" [[deps.Hwloc_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "ca0f6bf568b4bfc807e7537f081c81e35ceca114" +git-tree-sha1 = "1d334207121865ac8c1c97eb7f42d0339e4635bf" uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8" -version = "2.10.0+0" +version = "2.11.0+0" [[deps.IfElse]] git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" @@ -414,15 +426,21 @@ version = "0.2.1" [[deps.InlineStrings]] deps = ["Parsers"] -git-tree-sha1 = "9cc2baf75c6d09f9da536ddf58eb2f29dedaf461" +git-tree-sha1 = "86356004f30f8e737eff143d57d41bd580e437aa" uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48" -version = "1.4.0" +version = "1.4.1" + + [deps.InlineStrings.extensions] + ArrowTypesExt = "ArrowTypes" + + [deps.InlineStrings.weakdeps] + ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd" [[deps.IntelOpenMP_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "5fdf2fe6724d8caabf43b557b84ce53f3b7e2f6b" +git-tree-sha1 = "14eb2b542e748570b56446f4c50fbfb2306ebc45" uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" -version = "2024.0.2+0" +version = "2024.2.0+0" [[deps.InteractiveUtils]] deps = ["Markdown"] @@ -430,9 +448,9 @@ uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" [[deps.InverseFunctions]] deps = ["Test"] -git-tree-sha1 = "896385798a8d49a255c398bd49162062e4a4c435" +git-tree-sha1 = "e7cbed5032c4c397a6ac23d1493f3289e01231c4" uuid = "3587e190-3f89-42d0-90ee-14403ec27112" -version = "0.1.13" +version = "0.1.14" weakdeps = ["Dates"] [deps.InverseFunctions.extensions] @@ -460,10 +478,10 @@ uuid = "82899510-4779-5014-852e-03e436cf321d" version = "1.0.0" [[deps.JLD2]] -deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "PrecompileTools", "Printf", "Reexport", "Requires", "TranscodingStreams", "UUIDs"] -git-tree-sha1 = "5ea6acdd53a51d897672edb694e3cc2912f3f8a7" +deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "PrecompileTools", "Reexport", "Requires", "TranscodingStreams", "UUIDs", "Unicode"] +git-tree-sha1 = "84642bc18a79d715b39d3724b03cbdd2e7d48c62" uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -version = "0.4.46" +version = "0.4.49" [[deps.JLLWrappers]] deps = ["Artifacts", "Preferences"] @@ -491,9 +509,9 @@ version = "0.2.1+0" [[deps.KernelAbstractions]] deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "Requires", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"] -git-tree-sha1 = "ed7167240f40e62d97c1f5f7735dea6de3cc5c49" +git-tree-sha1 = "d0448cebd5919e06ca5edc7a264631790de810ec" uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c" -version = "0.9.18" +version = "0.9.22" [deps.KernelAbstractions.extensions] EnzymeExt = "EnzymeCore" @@ -503,9 +521,9 @@ version = "0.9.18" [[deps.LLVM]] deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"] -git-tree-sha1 = "839c82932db86740ae729779e610f07a1640be9a" +git-tree-sha1 = "020abd49586480c1be84f57da0017b5d3db73f7c" uuid = "929cbde3-209d-540e-8aea-75f648917ca0" -version = "6.6.3" +version = "8.0.0" weakdeps = ["BFloat16s"] [deps.LLVM.extensions] @@ -513,9 +531,9 @@ weakdeps = ["BFloat16s"] [[deps.LLVMExtra_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"] -git-tree-sha1 = "88b916503aac4fb7f701bb625cd84ca5dd1677bc" +git-tree-sha1 = "c2636c264861edc6d305e6b4d528f09566d24c5e" uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab" -version = "0.0.29+0" +version = "0.0.30+0" [[deps.LLVMLoopInfo]] git-tree-sha1 = "2e5c102cfc41f48ae4740c7eca7743cc7e7b75ea" @@ -585,9 +603,9 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[deps.LogExpFunctions]] deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "18144f3e9cbe9b15b070288eef858f71b291ce37" +git-tree-sha1 = "a2d09619db4e765091ee5c6ffe8872849de0feea" uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.27" +version = "0.3.28" [deps.LogExpFunctions.extensions] LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" @@ -609,10 +627,10 @@ uuid = "5ced341a-0733-55b8-9ab6-a4889d929147" version = "1.9.4+0" [[deps.MKL_jll]] -deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl"] -git-tree-sha1 = "72dc3cf284559eb8f53aa593fe62cb33f83ed0c0" +deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"] +git-tree-sha1 = "f046ccd0c6db2832a9f639e2c669c6fe867e5f4f" uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" -version = "2024.0.0+0" +version = "2024.2.0+0" [[deps.MPI]] deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "PkgVersion", "PrecompileTools", "Requires", "Serialization", "Sockets"] @@ -630,21 +648,21 @@ version = "0.20.16" [[deps.MPICH_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] -git-tree-sha1 = "656036b9ed6f942d35e536e249600bc31d0f9df8" +git-tree-sha1 = "4099bb6809ac109bfc17d521dad33763bcf026b7" uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4" -version = "4.2.0+0" +version = "4.2.1+1" [[deps.MPIPreferences]] deps = ["Libdl", "Preferences"] -git-tree-sha1 = "8f6af051b9e8ec597fa09d8885ed79fd582f33c9" +git-tree-sha1 = "c105fe467859e7f6e9a852cb15cb4301126fac07" uuid = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267" -version = "0.1.10" +version = "0.1.11" [[deps.MPItrampoline_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] -git-tree-sha1 = "77c3bd69fdb024d75af38713e883d0f249ce19c2" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] +git-tree-sha1 = "8c35d5420193841b2f367e658540e8d9e0601ed0" uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748" -version = "5.3.2+0" +version = "5.4.0+0" [[deps.MacroTools]] deps = ["Markdown", "Random"] @@ -669,9 +687,9 @@ version = "10.1.4+2" [[deps.Missings]] deps = ["DataAPI"] -git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272" +git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d" uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" -version = "1.1.0" +version = "1.2.0" [[deps.Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" @@ -682,9 +700,9 @@ version = "2023.1.10" [[deps.NCDatasets]] deps = ["CFTime", "CommonDataModel", "DataStructures", "Dates", "DiskArrays", "NetCDF_jll", "NetworkOptions", "Printf"] -git-tree-sha1 = "d40d24d12f710c39d3a66be99c567ce0032f28a7" +git-tree-sha1 = "a640912695952b074672edb5f9aaee2f7f9fd59a" uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" -version = "0.14.3" +version = "0.14.4" [[deps.NVTX]] deps = ["Colors", "JuliaNVTXCallbacks_jll", "Libdl", "NVTX_jll"] @@ -722,9 +740,9 @@ version = "1.2.0" [[deps.Oceananigans]] deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"] -git-tree-sha1 = "4672af7242405313743af45168bfce3d87b84b2c" +git-tree-sha1 = "942c78ad3c95e47bcbab3eeacb0b8a846bcfb33b" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" -version = "0.90.11" +version = "0.91.4" [deps.Oceananigans.extensions] OceananigansEnzymeExt = "Enzyme" @@ -733,9 +751,9 @@ version = "0.90.11" Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" [[deps.OffsetArrays]] -git-tree-sha1 = "6a731f2b5c03157418a20c12195eb4b74c8f8621" +git-tree-sha1 = "1a27764e945a152f7ca7efa04de513d473e9542e" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.13.0" +version = "1.14.1" weakdeps = ["Adapt"] [deps.OffsetArrays.extensions] @@ -752,16 +770,16 @@ uuid = "05823500-19ac-5b8b-9628-191a04bc5112" version = "0.8.1+2" [[deps.OpenMPI_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "PMIx_jll", "TOML", "Zlib_jll", "libevent_jll", "prrte_jll"] -git-tree-sha1 = "f46caf663e069027a06942d00dced37f1eb3d8ad" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML", "Zlib_jll"] +git-tree-sha1 = "a9de2f1fc98b92f8856c640bf4aec1ac9b2a0d86" uuid = "fe0851c0-eecd-5654-98d4-656369965a5c" -version = "5.0.2+0" +version = "5.0.3+0" [[deps.OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "3da7367955dcc5c54c1ba4d402ccdc09a1a3e046" +git-tree-sha1 = "a028ee3cb5641cccc4c24e90c36b0a4f7707bdf5" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.0.13+1" +version = "3.0.14+0" [[deps.OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] @@ -780,12 +798,6 @@ git-tree-sha1 = "2cd396108e178f3ae8dedbd8e938a18726ab2fbf" uuid = "c2071276-7c44-58a7-b746-946036e04d0a" version = "0.24.1+0" -[[deps.PMIx_jll]] -deps = ["Artifacts", "Hwloc_jll", "JLLWrappers", "Libdl", "Zlib_jll", "libevent_jll"] -git-tree-sha1 = "8b3b19351fa24791f94d7ae85faf845ca1362541" -uuid = "32165bc3-0280-59bc-8c0b-c33b6203efab" -version = "4.2.7+0" - [[deps.PackageExtensionCompat]] git-tree-sha1 = "fb28e33b8a95c4cee25ce296c817d89cc2e53518" uuid = "65ce6f38-6b18-4e1d-a461-8949797d7930" @@ -800,15 +812,17 @@ version = "2.8.1" [[deps.PencilArrays]] deps = ["Adapt", "JSON3", "LinearAlgebra", "MPI", "OffsetArrays", "Random", "Reexport", "StaticArrayInterface", "StaticArrays", "StaticPermutations", "Strided", "TimerOutputs", "VersionParsing"] -git-tree-sha1 = "6510e851700a851944f7ffa5cd990cced4802ad2" +git-tree-sha1 = "fa85ac32172d96cfdb91dbc53e8e57007e5a2b5a" uuid = "0e08944d-e94e-41b1-9406-dcf66b6a9d2e" -version = "0.19.3" +version = "0.19.5" [deps.PencilArrays.extensions] + PencilArraysAMDGPUExt = ["AMDGPU"] PencilArraysDiffEqExt = ["DiffEqBase"] PencilArraysHDF5Ext = ["HDF5"] [deps.PencilArrays.weakdeps] + AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" @@ -849,9 +863,9 @@ version = "1.4.3" [[deps.PrettyTables]] deps = ["Crayons", "LaTeXStrings", "Markdown", "PrecompileTools", "Printf", "Reexport", "StringManipulation", "Tables"] -git-tree-sha1 = "88b895d13d53b5577fd53379d913b9ab9ac82660" +git-tree-sha1 = "66b20dd35966a748321d3b2537c4584cf40387c7" uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" -version = "2.3.1" +version = "2.3.2" [[deps.Printf]] deps = ["Unicode"] @@ -938,9 +952,9 @@ version = "2.1.5" [[deps.Rotations]] deps = ["LinearAlgebra", "Quaternions", "Random", "StaticArrays"] -git-tree-sha1 = "2a0a5d8569f481ff8840e3b7c84bbf188db6a3fe" +git-tree-sha1 = "5680a9276685d392c87407df00d57c9924d9f11e" uuid = "6038ab10-8711-5258-84ad-4b1120ba62dc" -version = "1.7.0" +version = "1.7.1" weakdeps = ["RecipesBase"] [deps.Rotations.extensions] @@ -963,9 +977,9 @@ version = "0.3.4" [[deps.SentinelArrays]] deps = ["Dates", "Random"] -git-tree-sha1 = "0e7508ff27ba32f26cd459474ca2ede1bc10991f" +git-tree-sha1 = "ff11acffdb082493657550959d4feb4b6149e73a" uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c" -version = "1.4.1" +version = "1.4.5" [[deps.Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -986,25 +1000,25 @@ version = "1.10.0" [[deps.SpecialFunctions]] deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d" +git-tree-sha1 = "2f5d4697f21388cbe1ff299430dd169ef97d7e14" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.3.1" +version = "2.4.0" weakdeps = ["ChainRulesCore"] [deps.SpecialFunctions.extensions] SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" [[deps.Static]] -deps = ["IfElse"] -git-tree-sha1 = "d2fdac9ff3906e27f7a618d47b676941baa6c80c" +deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools"] +git-tree-sha1 = "0bbff21027dd8a107551847528127b62a35f7594" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.8.10" +version = "1.1.0" [[deps.StaticArrayInterface]] deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"] -git-tree-sha1 = "5d66818a39bb04bf328e92bc933ec5b4ee88e436" +git-tree-sha1 = "8963e5a083c837531298fc41599182a759a87a6d" uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718" -version = "1.5.0" +version = "1.5.1" weakdeps = ["OffsetArrays", "StaticArrays"] [deps.StaticArrayInterface.extensions] @@ -1013,9 +1027,9 @@ weakdeps = ["OffsetArrays", "StaticArrays"] [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "bf074c045d3d5ffd956fa0a461da38a44685d6b2" +git-tree-sha1 = "eeafab08ae20c62c44c8399ccb9354a04b80db50" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.9.3" +version = "1.9.7" weakdeps = ["ChainRulesCore", "Statistics"] [deps.StaticArrays.extensions] @@ -1023,9 +1037,9 @@ weakdeps = ["ChainRulesCore", "Statistics"] StaticArraysStatisticsExt = "Statistics" [[deps.StaticArraysCore]] -git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d" +git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682" uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" -version = "1.4.2" +version = "1.4.3" [[deps.StaticPermutations]] git-tree-sha1 = "193c3daa18ff3e55c1dae66acb6a762c4a3bdb0b" @@ -1045,15 +1059,15 @@ version = "1.7.0" [[deps.Strided]] deps = ["LinearAlgebra", "StridedViews", "TupleTools"] -git-tree-sha1 = "40c69be0e1b72ee2f42923b7d1ff13e0b04e675c" +git-tree-sha1 = "bd9bd1c70cfc115cc3a30213fc725125a6b43652" uuid = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" -version = "2.0.4" +version = "2.1.0" [[deps.StridedViews]] deps = ["LinearAlgebra", "PackageExtensionCompat"] -git-tree-sha1 = "5b765c4e401693ab08981989f74a36a010aa1d8e" +git-tree-sha1 = "2917996ce0fa6b8a3a85240a5e9ff930e2aeaa43" uuid = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143" -version = "0.2.2" +version = "0.3.1" weakdeps = ["CUDA"] [deps.StridedViews.extensions] @@ -1133,14 +1147,14 @@ uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [[deps.TimerOutputs]] deps = ["ExprTools", "Printf"] -git-tree-sha1 = "f548a9e9c490030e545f72074a41edfd0e5bcdd7" +git-tree-sha1 = "5a13ae8a41237cff5ecf34f73eb1b8f42fff6531" uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" -version = "0.5.23" +version = "0.5.24" [[deps.TranscodingStreams]] -git-tree-sha1 = "71509f04d045ec714c4748c785a59045c3736349" +git-tree-sha1 = "60df3f8126263c0d6b357b9a1017bb94f53e3582" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" -version = "0.10.7" +version = "0.11.0" weakdeps = ["Random", "Test"] [deps.TranscodingStreams.extensions] @@ -1165,9 +1179,9 @@ version = "0.2.1" [[deps.UnsafeAtomicsLLVM]] deps = ["LLVM", "UnsafeAtomics"] -git-tree-sha1 = "323e3d0acf5e78a56dfae7bd8928c989b4f3083e" +git-tree-sha1 = "bf2c553f25e954a9b38c9c0593a59bb13113f9e5" uuid = "d80eeb9a-aca5-4d75-85e5-170c8b632249" -version = "0.1.3" +version = "0.1.5" [[deps.VersionParsing]] git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868" @@ -1176,9 +1190,9 @@ version = "1.3.0" [[deps.XML2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"] -git-tree-sha1 = "532e22cf7be8462035d092ff21fada7527e2c488" +git-tree-sha1 = "d9717ce3518dc68a99e6b96300813760d887a01d" uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" -version = "2.12.6+0" +version = "2.13.1+0" [[deps.XZ_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -1208,12 +1222,6 @@ deps = ["Artifacts", "Libdl"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" version = "5.8.0+1" -[[deps.libevent_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "OpenSSL_jll"] -git-tree-sha1 = "f04ec6d9a186115fb38f858f05c0c4e1b7fc9dcb" -uuid = "1080aeaf-3a6a-583e-a51c-c537b09f60ec" -version = "2.1.13+1" - [[deps.libzip_jll]] deps = ["Artifacts", "Bzip2_jll", "GnuTLS_jll", "JLLWrappers", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"] git-tree-sha1 = "3282b7d16ae7ac3e57ec2f3fa8fafb564d8f9f7f" @@ -1225,13 +1233,13 @@ deps = ["Artifacts", "Libdl"] uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" version = "1.52.0+1" +[[deps.oneTBB_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "7d0ea0f4895ef2f5cb83645fa689e52cb55cf493" +uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e" +version = "2021.12.0+0" + [[deps.p7zip_jll]] deps = ["Artifacts", "Libdl"] uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" version = "17.4.0+2" - -[[deps.prrte_jll]] -deps = ["Artifacts", "Hwloc_jll", "JLLWrappers", "Libdl", "PMIx_jll", "libevent_jll"] -git-tree-sha1 = "5adb2d7a18a30280feb66cad6f1a1dfdca2dc7b0" -uuid = "eb928a42-fffd-568d-ab9c-3f5d54fc65b9" -version = "3.0.2+0" diff --git a/docs/src/library/public.md b/docs/src/library/public.md index 2b9566c..d64b170 100644 --- a/docs/src/library/public.md +++ b/docs/src/library/public.md @@ -11,10 +11,23 @@ Modules = [ClimaSeaIce] Private = false ``` +## ClimaSeaIce.SeaIceThermodynamics + +```@autodocs +Modules = [ClimaSeaIce.SeaIceThermodynamics] +Private = false +``` + +## ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions + +```@autodocs +Modules = [ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions] +Private = false +``` + ## ClimaSeaIce.EnthalpyMethodSeaIceModels ```@autodocs Modules = [ClimaSeaIce.EnthalpyMethodSeaIceModels] Private = false ``` - diff --git a/examples/freezing_bucket.jl b/examples/freezing_bucket.jl index edc38cd..143d287 100644 --- a/examples/freezing_bucket.jl +++ b/examples/freezing_bucket.jl @@ -49,17 +49,22 @@ phase_transitions = PhaseTransitions(; ice_heat_capacity, ice_density) top_temperature = -10 # ᵒC top_heat_boundary_condition = PrescribedTemperature(-10) +# Construct the thermodynamics of sea ice, for this we use a simple +# slab sea ice representation of thermodynamics + +ice_thermodynamics = SlabSeaIceThermodynamics(grid; + internal_heat_flux, + phase_transitions, + top_heat_boundary_condition) + # Then we assemble it all into a model, -model = SlabSeaIceModel(grid; - internal_heat_flux, - phase_transitions, - top_heat_boundary_condition) +model = SeaIceModel(grid; ice_thermodynamics) -# Note that the default bottom heat boundary condition for `SlabSeaIceModel` is +# Note that the default bottom heat boundary condition for `SlabSeaIceThermodynamics` is # `IceWaterThermalEquilibrium` with freshwater. That's what we want! -model.heat_boundary_conditions.bottom +model.ice_thermodynamics.heat_boundary_conditions.bottom # Ok, we're ready to freeze the bucket for 10 straight days with an initial ice # thickness of 1 cm, diff --git a/src/ClimaSeaIce.jl b/src/ClimaSeaIce.jl index 0215019..e637c18 100644 --- a/src/ClimaSeaIce.jl +++ b/src/ClimaSeaIce.jl @@ -1,142 +1,37 @@ """ Ocean 🌊 Sea ice component of CliMa's Earth system model. """ module ClimaSeaIce -import Oceananigans.TimeSteppers: time_step! - -export - MeltingConstrainedFluxBalance, - PrescribedTemperature, - RadiativeEmission, - PhaseTransitions, - ConductiveFlux, - FluxFunction, - SlabSeaIceModel - -##### -##### A bit of thermodynamics to start the day -##### - -struct LinearLiquidus{FT} - freshwater_melting_temperature :: FT - slope :: FT -end - -""" - LinearLiquidus(FT=Float64, - slope = 0.054, # psu / ᵒC - freshwater_melting_temperature = 0) # ᵒC - -Return a linear model for the dependence of the melting temperature of -saltwater on salinity, - -```math -Tₘ(S) = T₀ - m S , -``` - -where ``Tₘ(S)`` is the melting temperature as a function of salinity ``S``, -``T₀`` is the melting temperature of freshwater, and ``m`` is the ratio -between the melting temperature and salinity (in other words the linear model -should be thought of as defining ``m`` and could be written ``m ≡ (T₀ - Tₘ) / S``. -The signs are arranged so that ``m > 0`` for saltwater). - -The defaults assume that salinity is given in practical salinity units `psu` and -temperature is in degrees Celsius. - -Note: the function `melting_temperature(liquidus, salinity)` returns the -melting temperature given `salinity`. -""" -function LinearLiquidus(FT::DataType=Float64; - slope = 0.054, # psu / ᵒC - freshwater_melting_temperature = 0) # ᵒC - - return LinearLiquidus(convert(FT, freshwater_melting_temperature), - convert(FT, slope)) -end - -@inline function melting_temperature(liquidus::LinearLiquidus, salinity) - return liquidus.freshwater_melting_temperature - liquidus.slope * salinity -end - -struct PhaseTransitions{FT, L} - ice_density :: FT - ice_heat_capacity :: FT - liquid_density :: FT - liquid_heat_capacity :: FT - reference_latent_heat :: FT - reference_temperature :: FT - liquidus :: L -end - -""" - PhaseTransitions(FT=Float64, - ice_density = 917, # kg m⁻³ - ice_heat_capacity = 2000, # J / (kg ᵒC) - liquid_density = 999.8, # kg m⁻³ - liquid_heat_capacity = 4186, # J / (kg ᵒC) - reference_latent_heat = 334e3 # J kg⁻³ - liquidus = LinearLiquidus(FT)) # default assumes psu, ᵒC - -Return a representation of transitions between the solid and liquid phases -of salty water: in other words, the freezing and melting of sea ice. - -The latent heat of fusion ``ℒ(T)`` (more simply just "latent heat") is -a function of temperature ``T`` via - -```math -ρᵢ ℒ(T) = ρᵢ ℒ₀ + (ρ_ℓ c_ℓ - ρᵢ cᵢ) (T - T₀) -``` - -where ``ρᵢ`` is the `ice_density`, ``ρ_ℓ`` is the liquid density, -``cᵢ`` is the heat capacity of ice, and ``c_ℓ`` is the heat capacity of -liquid, and ``T₀`` is a reference temperature, all of which are assumed constant. - -The default `liquidus` assumes that salinity has practical salinity units (psu) -and that temperature is degrees Celsius. -""" -@inline function PhaseTransitions(FT=Float64; - ice_density = 917, # kg m⁻³ - ice_heat_capacity = 2000, # J / (kg ᵒC) - liquid_density = 999.8, # kg m⁻³ - liquid_heat_capacity = 4186, # J / (kg ᵒC) - reference_latent_heat = 334e3, # J kg⁻³ - reference_temperature = 0, # ᵒC - liquidus = LinearLiquidus(FT)) - - return PhaseTransitions(convert(FT, ice_density), - convert(FT, ice_heat_capacity), - convert(FT, liquid_density), - convert(FT, liquid_heat_capacity), - convert(FT, reference_latent_heat), - convert(FT, reference_temperature), - liquidus) -end - -@inline function latent_heat(thermo::PhaseTransitions, T) - T₀ = thermo.reference_temperature - ℒ₀ = thermo.reference_latent_heat - ρᵢ = thermo.ice_density - ρℓ = thermo.liquid_density - cᵢ = thermo.ice_heat_capacity - cℓ = thermo.liquid_heat_capacity - - return ρℓ * ℒ₀ + (ρℓ * cℓ - ρᵢ * cᵢ) * (T - T₀) -end +using Oceananigans +using Oceananigans.Utils: prettysummary +using Oceananigans.TimeSteppers: Clock +using Oceananigans.Fields: field, Field, Center, ZeroField, ConstantField + +# Simulations interface +import Oceananigans: fields, prognostic_fields +import Oceananigans.Fields: set! +import Oceananigans.Models: AbstractModel +import Oceananigans.OutputWriters: default_included_properties +import Oceananigans.Simulations: reset!, initialize!, iteration +import Oceananigans.TimeSteppers: time_step!, update_state! +import Oceananigans.Utils: prettytime + +export SeaIceModel, + MeltingConstrainedFluxBalance, + PrescribedTemperature, + RadiativeEmission, + PhaseTransitions, + ConductiveFlux, + FluxFunction, + SlabSeaIceThermodynamics struct ForwardEulerTimestepper end -include("HeatBoundaryConditions/HeatBoundaryConditions.jl") - -using .HeatBoundaryConditions: - IceWaterThermalEquilibrium, - MeltingConstrainedFluxBalance, - RadiativeEmission, - FluxFunction, - PrescribedTemperature - -include("EnthalpyMethodSeaIceModels.jl") -include("SlabSeaIceModels/SlabSeaIceModels.jl") +include("SeaIceThermodynamics/SeaIceThermodynamics.jl") +include("sea_ice_model.jl") +include("tracer_tendency_kernel_functions.jl") +include("time_stepping.jl") +include("EnthalpyMethodSeaIceModel.jl") -using .EnthalpyMethodSeaIceModels: EnthalpyMethodSeaIceModel -using .SlabSeaIceModels: SlabSeaIceModel, ConductiveFlux +using .SeaIceThermodynamics end # module diff --git a/src/EnthalpyMethodSeaIceModels.jl b/src/EnthalpyMethodSeaIceModel.jl similarity index 89% rename from src/EnthalpyMethodSeaIceModels.jl rename to src/EnthalpyMethodSeaIceModel.jl index 78f6a61..da79b79 100644 --- a/src/EnthalpyMethodSeaIceModels.jl +++ b/src/EnthalpyMethodSeaIceModel.jl @@ -51,16 +51,21 @@ end const reference_density = 999.8 # kg m⁻³ """ - EnthalpyMethodSeaIceModel(; grid, kw...) + EnthalpyMethodSeaIceModel(; grid, + closure = default_closure(grid), + ice_heat_capacity = 2090.0 / reference_density, + water_heat_capacity = 3991.0 / reference_density, + fusion_enthalpy = 3.3e5 / reference_density, + boundary_conditions = NamedTuple()) Return a thermodynamic model for ice sandwiched between an atmosphere and ocean on an Eulerian grid. """ function EnthalpyMethodSeaIceModel(; grid, - closure = default_closure(grid), - ice_heat_capacity = 2090.0 / reference_density, - water_heat_capacity = 3991.0 / reference_density, - fusion_enthalpy = 3.3e5 / reference_density, - boundary_conditions = NamedTuple()) + closure = default_closure(grid), + ice_heat_capacity = 2090.0 / reference_density, + water_heat_capacity = 3991.0 / reference_density, + fusion_enthalpy = 3.3e5 / reference_density, + boundary_conditions = NamedTuple()) # Prognostic fields: temperature, enthalpy, porosity field_names = (:T, :H, :ϕ) @@ -72,7 +77,7 @@ function EnthalpyMethodSeaIceModel(; grid, state = TracerFields(field_names, grid, user_boundary_conditions) tendencies = (; H=CenterField(grid)) - clock = Clock{eltype(grid)}(0, 0, 1) + clock = Clock{eltype(grid)}(time = 0) return EnthalpyMethodSeaIceModel(grid, nothing, @@ -249,4 +254,3 @@ end @inbounds temperature_flux(i, j, k, grid, κ, T) = - κᶜᶜᶜ_∂zᶜᶜᶠT(i, j, k, grid, κ, T) end # module - diff --git a/src/HeatBoundaryConditions/HeatBoundaryConditions.jl b/src/SeaIceThermodynamics/HeatBoundaryConditions/HeatBoundaryConditions.jl similarity index 76% rename from src/HeatBoundaryConditions/HeatBoundaryConditions.jl rename to src/SeaIceThermodynamics/HeatBoundaryConditions/HeatBoundaryConditions.jl index 2f5c92a..0bc942e 100644 --- a/src/HeatBoundaryConditions/HeatBoundaryConditions.jl +++ b/src/SeaIceThermodynamics/HeatBoundaryConditions/HeatBoundaryConditions.jl @@ -1,5 +1,11 @@ module HeatBoundaryConditions +export MeltingConstrainedFluxBalance, + PrescribedTemperature, + RadiativeEmission, + ConductiveFlux, + FluxFunction + using Adapt """ @@ -19,4 +25,3 @@ include("top_heat_boundary_conditions.jl") include("boundary_fluxes.jl") end - diff --git a/src/HeatBoundaryConditions/bottom_heat_boundary_conditions.jl b/src/SeaIceThermodynamics/HeatBoundaryConditions/bottom_heat_boundary_conditions.jl similarity index 96% rename from src/HeatBoundaryConditions/bottom_heat_boundary_conditions.jl rename to src/SeaIceThermodynamics/HeatBoundaryConditions/bottom_heat_boundary_conditions.jl index 940d398..4adaacb 100644 --- a/src/HeatBoundaryConditions/bottom_heat_boundary_conditions.jl +++ b/src/SeaIceThermodynamics/HeatBoundaryConditions/bottom_heat_boundary_conditions.jl @@ -1,4 +1,4 @@ -using ClimaSeaIce: melting_temperature +using ClimaSeaIce.SeaIceThermodynamics: melting_temperature ##### ##### Bottom heat boundary conditions @@ -41,17 +41,16 @@ end @inline function bottom_flux_imbalance(i, j, grid, bottom_heat_bc, top_temperature, internal_fluxes, external_fluxes, clock, model_fields) - # + # # ice ↑ Qi ≡ internal_fluxes. Example: Qi = - k ∂z T # |⎴⎴⎴| # ----------------------- ↕ hᵇ → ℒ ∂t hᵇ = δQ _given_ T = Tₘ # |⎵⎵⎵| # water ↑ Qx ≡ external_fluxes - # + # Qi = getflux(internal_fluxes, i, j, grid, top_temperature, clock, model_fields) Qx = getflux(external_fluxes, i, j, grid, top_temperature, clock, model_fields) return Qi - Qx end - diff --git a/src/HeatBoundaryConditions/boundary_fluxes.jl b/src/SeaIceThermodynamics/HeatBoundaryConditions/boundary_fluxes.jl similarity index 99% rename from src/HeatBoundaryConditions/boundary_fluxes.jl rename to src/SeaIceThermodynamics/HeatBoundaryConditions/boundary_fluxes.jl index 56d7d5a..93d5b34 100644 --- a/src/HeatBoundaryConditions/boundary_fluxes.jl +++ b/src/SeaIceThermodynamics/HeatBoundaryConditions/boundary_fluxes.jl @@ -32,7 +32,7 @@ struct FluxFunction{P, T, F} T = SurfaceTemperatureDependence P = typeof(parameters) F = typeof(func) - new{P, T, F}(func, parameters) + return new{P, T, F}(func, parameters) end end @@ -140,4 +140,3 @@ function flux_summary(fluxes::Tuple, padchar=" ") "$padchar └── ", summary(fluxes[end])) end end - diff --git a/src/HeatBoundaryConditions/top_heat_boundary_conditions.jl b/src/SeaIceThermodynamics/HeatBoundaryConditions/top_heat_boundary_conditions.jl similarity index 98% rename from src/HeatBoundaryConditions/top_heat_boundary_conditions.jl rename to src/SeaIceThermodynamics/HeatBoundaryConditions/top_heat_boundary_conditions.jl index 9fe79d7..898f014 100644 --- a/src/HeatBoundaryConditions/top_heat_boundary_conditions.jl +++ b/src/SeaIceThermodynamics/HeatBoundaryConditions/top_heat_boundary_conditions.jl @@ -41,7 +41,7 @@ The residual flux is consumed by the cost of transforming ice into liquid water, and is related to the rate of change of ice thickness, ``h``, by ```math -d/dt hₛ = δQ / ℒ(Tᵤ) +\\frac{dhₛ}{dt} = δQ / ℒ(Tᵤ) ``` where ``ℒ(Tᵤ)`` is the latent heat, equal to the different between @@ -93,9 +93,8 @@ using RootSolvers: SecantMethod, find_zero, CompactSolution flux_balance(T) = getflux(external_fluxes, i, j, grid, T, clock, model_fields) - getflux(internal_fluxes, i, j, grid, T, clock, model_fields) - + solution = find_zero(flux_balance, method, solution_type) - + return solution.root end - diff --git a/src/SeaIceThermodynamics/SeaIceThermodynamics.jl b/src/SeaIceThermodynamics/SeaIceThermodynamics.jl new file mode 100644 index 0000000..57b6709 --- /dev/null +++ b/src/SeaIceThermodynamics/SeaIceThermodynamics.jl @@ -0,0 +1,158 @@ +module SeaIceThermodynamics + +export SlabSeaIceThermodynamics, + PhaseTransitions, + MeltingConstrainedFluxBalance, + PrescribedTemperature, + RadiativeEmission, + ConductiveFlux, + FluxFunction + +using Adapt + +##### +##### A bit of thermodynamics to start the day +##### + +struct LinearLiquidus{FT} + freshwater_melting_temperature :: FT + slope :: FT +end + +""" + LinearLiquidus(FT=Float64, + slope = 0.054, # psu / ᵒC + freshwater_melting_temperature = 0) # ᵒC + +Return a linear model for the dependence of the melting temperature of +saltwater on salinity, + +```math +Tₘ(S) = T₀ - m S , +``` + +where ``Tₘ(S)`` is the melting temperature as a function of salinity ``S``, +``T₀`` is the melting temperature of freshwater, and ``m`` is the ratio +between the melting temperature and salinity (in other words the linear model +should be thought of as defining ``m`` and could be written ``m ≡ (T₀ - Tₘ) / S``. +The signs are arranged so that ``m > 0`` for saltwater). + +The defaults assume that salinity is given in practical salinity units `psu` and +temperature is in degrees Celsius. + +Note: the function `melting_temperature(liquidus, salinity)` returns the +melting temperature given `salinity`. +""" +function LinearLiquidus(FT::DataType=Float64; + slope = 0.054, # psu / ᵒC + freshwater_melting_temperature = 0) # ᵒC + + return LinearLiquidus(convert(FT, freshwater_melting_temperature), + convert(FT, slope)) +end + +@inline function melting_temperature(liquidus::LinearLiquidus, salinity) + return liquidus.freshwater_melting_temperature - liquidus.slope * salinity +end + +struct PhaseTransitions{FT, L} + ice_density :: FT + ice_heat_capacity :: FT + liquid_density :: FT + liquid_heat_capacity :: FT + reference_latent_heat :: FT + reference_temperature :: FT + liquidus :: L +end + +""" + PhaseTransitions(FT=Float64, + ice_density = 917, # kg m⁻³ + ice_heat_capacity = 2000, # J / (kg ᵒC) + liquid_density = 999.8, # kg m⁻³ + liquid_heat_capacity = 4186, # J / (kg ᵒC) + reference_latent_heat = 334e3 # J kg⁻³ + liquidus = LinearLiquidus(FT)) # default assumes psu, ᵒC + +Return a representation of transitions between the solid and liquid phases +of salty water: in other words, the freezing and melting of sea ice. + +The latent heat of fusion ``ℒ(T)`` (more simply just "latent heat") is +a function of temperature ``T`` via + +```math +ρᵢ ℒ(T) = ρᵢ ℒ₀ + (ρ_ℓ c_ℓ - ρᵢ cᵢ) (T - T₀) +``` + +where ``ρᵢ`` is the `ice_density`, ``ρ_ℓ`` is the liquid density, +``cᵢ`` is the heat capacity of ice, and ``c_ℓ`` is the heat capacity of +liquid, and ``T₀`` is a reference temperature, all of which are assumed constant. + +The default `liquidus` assumes that salinity has practical salinity units (psu) +and that temperature is degrees Celsius. +""" +@inline function PhaseTransitions(FT=Float64; + ice_density = 917, # kg m⁻³ + ice_heat_capacity = 2000, # J / (kg ᵒC) + liquid_density = 999.8, # kg m⁻³ + liquid_heat_capacity = 4186, # J / (kg ᵒC) + reference_latent_heat = 334e3, # J kg⁻³ + reference_temperature = 0, # ᵒC + liquidus = LinearLiquidus(FT)) + + return PhaseTransitions(convert(FT, ice_density), + convert(FT, ice_heat_capacity), + convert(FT, liquid_density), + convert(FT, liquid_heat_capacity), + convert(FT, reference_latent_heat), + convert(FT, reference_temperature), + liquidus) +end + +@inline function latent_heat(thermo::PhaseTransitions, T) + T₀ = thermo.reference_temperature + ℒ₀ = thermo.reference_latent_heat + ρᵢ = thermo.ice_density + ρℓ = thermo.liquid_density + cᵢ = thermo.ice_heat_capacity + cℓ = thermo.liquid_heat_capacity + + return ρℓ * ℒ₀ + (ρℓ * cℓ - ρᵢ * cᵢ) * (T - T₀) +end + +@inline external_top_heat_flux(ice_thermodynamics, top_heat_flux) = top_heat_flux + +# Fallback for no thermodynamics +@inline thickness_thermodynamic_tendency(i, j, grid, args...) = zero(grid) + +include("HeatBoundaryConditions/HeatBoundaryConditions.jl") + +using .HeatBoundaryConditions: + IceWaterThermalEquilibrium, + MeltingConstrainedFluxBalance, + RadiativeEmission, + FluxFunction, + PrescribedTemperature, + getflux + +using Oceananigans.Utils: prettysummary +using Oceananigans.TimeSteppers: Clock +using Oceananigans.Fields: field, Field, Center, ZeroField, ConstantField + +# Simulations interface +import Oceananigans: fields, prognostic_fields +import Oceananigans.Fields: set! +import Oceananigans.Models: AbstractModel +import Oceananigans.OutputWriters: default_included_properties +import Oceananigans.Simulations: reset!, initialize!, iteration +import Oceananigans.TimeSteppers: time_step!, update_state! +import Oceananigans.Utils: prettytime + +# TODO: Fix this after this PR +# include("EnthalpyMethodThermodynamics.jl") + +include("slab_sea_ice_thermodynamics.jl") +include("slab_heat_and_tracer_fluxes.jl") +include("slab_thermodynamics_tendencies.jl") + +end diff --git a/src/SlabSeaIceModels/slab_heat_and_tracer_fluxes.jl b/src/SeaIceThermodynamics/slab_heat_and_tracer_fluxes.jl similarity index 100% rename from src/SlabSeaIceModels/slab_heat_and_tracer_fluxes.jl rename to src/SeaIceThermodynamics/slab_heat_and_tracer_fluxes.jl diff --git a/src/SeaIceThermodynamics/slab_sea_ice_thermodynamics.jl b/src/SeaIceThermodynamics/slab_sea_ice_thermodynamics.jl new file mode 100644 index 0000000..bbca290 --- /dev/null +++ b/src/SeaIceThermodynamics/slab_sea_ice_thermodynamics.jl @@ -0,0 +1,96 @@ +import Oceananigans: fields + +struct SlabSeaIceThermodynamics{ST, HBC, CF, P, MIT} + top_surface_temperature :: ST + heat_boundary_conditions :: HBC + # Internal flux + internal_heat_flux :: CF + # Melting and freezing stuff + phase_transitions :: P + ice_consolidation_thickness :: MIT +end + +Adapt.adapt_structure(to, t::SlabSeaIceThermodynamics) = + SlabSeaIceThermodynamics(Adapt.adapt(to, t.top_surface_temperature), + Adapt.adapt(to, t.heat_boundary_conditions), + Adapt.adapt(to, t.internal_heat_flux), + Adapt.adapt(to, t.phase_transitions), + Adapt.adapt(to, t.ice_consolidation_thickness)) + +const SSIT = SlabSeaIceThermodynamics + +Base.summary(therm::SSIT) = "SlabThermodynamics" + +function Base.show(io::IO, therm::SSIT) + print(io, "SlabSeaIceThermodynamics", '\n') + print(io, "├── top_surface_temperature: ", summary(therm.top_surface_temperature), '\n') + print(io, "└── minimium_ice_thickness: ", prettysummary(therm.ice_consolidation_thickness), '\n') +end + +fields(therm::SSIT) = (; Tu = therm.top_surface_temperature) + +""" + SlabSeaIceThermodynamics(grid; kw...) + +Pretty simple model for sea ice. +""" +function SlabSeaIceThermodynamics(grid; + ice_consolidation_thickness = 0.0, # m + top_surface_temperature = nothing, + top_heat_boundary_condition = MeltingConstrainedFluxBalance(), + bottom_heat_boundary_condition = IceWaterThermalEquilibrium(), + # Default internal flux: thermal conductivity of 2 kg m s⁻³ K⁻¹, appropriate for freshwater ice + internal_heat_flux = ConductiveFlux(eltype(grid), conductivity=2), + phase_transitions = PhaseTransitions(eltype(grid))) + + # TODO: pass `clock` into `field`, so functions can be time-dependent? + # Wrap ice_consolidation_thickness in a field + ice_consolidation_thickness = field((Center, Center, Nothing), ice_consolidation_thickness, grid) + + # Construct an internal heat flux function that captures the liquidus and + # bottom boundary condition. + parameters = (flux = internal_heat_flux, + liquidus = phase_transitions.liquidus, + bottom_heat_boundary_condition = bottom_heat_boundary_condition) + + internal_heat_flux_function = FluxFunction(slab_internal_heat_flux; + parameters, + top_temperature_dependent=true) + + # Construct default top temperature if one is not provided + if isnothing(top_surface_temperature) + # Check top boundary condition + if top_heat_boundary_condition isa PrescribedTemperature + top_surface_temperature = top_heat_boundary_condition.temperature + else # build the default + top_surface_temperature = Field{Center, Center, Nothing}(grid) + end + end + + # Convert to `field` (does nothing if it's already a Field) + top_surface_temperature = field((Center, Center, Nothing), top_surface_temperature, grid) + + heat_boundary_conditions = (top = top_heat_boundary_condition, + bottom = bottom_heat_boundary_condition) + + return SlabSeaIceThermodynamics(top_surface_temperature, + heat_boundary_conditions, + internal_heat_flux_function, + phase_transitions, + ice_consolidation_thickness) +end + +function external_top_heat_flux(thermodynamics::SlabSeaIceThermodynamics, + top_heat_flux) # Construct default top heat flux if one is not provided + if isnothing(top_heat_flux) + if thermodynamics.heat_boundary_conditions.top isa PrescribedTemperature + # Default: external top flux is in equilibrium with internal fluxes + top_heat_flux = thermodynamics.internal_heat_flux + else + # Default: no external top surface flux + top_heat_flux = 0 + end + end + + return top_heat_flux +end \ No newline at end of file diff --git a/src/SeaIceThermodynamics/slab_thermodynamics_tendencies.jl b/src/SeaIceThermodynamics/slab_thermodynamics_tendencies.jl new file mode 100644 index 0000000..c26cc84 --- /dev/null +++ b/src/SeaIceThermodynamics/slab_thermodynamics_tendencies.jl @@ -0,0 +1,67 @@ +using ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions: bottom_temperature, top_surface_temperature +import ClimaSeaIce.SeaIceThermodynamics: thickness_thermodynamic_tendency + +@inline function thickness_thermodynamic_tendency(i, j, k, grid, + ice_thickness, + ice_concentration, + thermodynamics::SlabSeaIceThermodynamics, + top_external_heat_flux, + bottom_external_heat_flux, + clock, model_fields) + + phase_transitions = thermodynamics.phase_transitions + + top_heat_bc = thermodynamics.heat_boundary_conditions.top + bottom_heat_bc = thermodynamics.heat_boundary_conditions.bottom + liquidus = phase_transitions.liquidus + + Qi = thermodynamics.internal_heat_flux + Qu = top_external_heat_flux + Qb = bottom_external_heat_flux + Tu = thermodynamics.top_surface_temperature + + @inbounds begin + hᶜ = thermodynamics.ice_consolidation_thickness[i, j, k] + hᵢ = ice_thickness[i, j, k] + end + + # Consolidation criteria + consolidated_ice = hᵢ >= hᶜ + + # Determine top surface temperature. + # Does this really fit here? + # This is updating the temperature inside the thermodynamics module + if !isa(top_heat_bc, PrescribedTemperature) # update surface temperature? + if consolidated_ice # slab is consolidated and has an independent surface temperature + Tu⁻ = @inbounds Tu[i, j, k] + Tuⁿ = top_surface_temperature(i, j, grid, top_heat_bc, Tu⁻, Qi, Qu, clock, model_fields) + else # slab is unconsolidated and does not have an independent surface temperature + Tuⁿ = bottom_temperature(i, j, grid, bottom_heat_bc, liquidus) + end + + @inbounds Tu[i, j, k] = Tuⁿ + end + + @inbounds Tuᵢ = Tu[i, j, k] + + Tbᵢ = bottom_temperature(i, j, grid, bottom_heat_bc, liquidus) + ℰb = latent_heat(phase_transitions, Tbᵢ) + ℰu = latent_heat(phase_transitions, Tuᵢ) + + # Retrieve fluxes + Quᵢ = getflux(Qu, i, j, grid, Tuᵢ, clock, model_fields) + Qiᵢ = getflux(Qi, i, j, grid, Tuᵢ, clock, model_fields) + Qbᵢ = getflux(Qb, i, j, grid, Tuᵢ, clock, model_fields) + + # If ice is consolidated, compute tendency for an ice slab; otherwise + # just add ocean fluxes from frazil ice formation or melting + slushy_Gh = - Qbᵢ / ℰb + + # Upper (top) and bottom interface velocities + wu = (Quᵢ - Qiᵢ) / ℰu # < 0 => melting + wb = (Qiᵢ - Qbᵢ) / ℰb # < 0 => freezing + + slabby_Gh = wu + wb + + return ifelse(consolidated_ice, slabby_Gh, slushy_Gh) +end \ No newline at end of file diff --git a/src/SlabSeaIceModels/SlabSeaIceModels.jl b/src/SlabSeaIceModels/SlabSeaIceModels.jl deleted file mode 100644 index 05d7b3c..0000000 --- a/src/SlabSeaIceModels/SlabSeaIceModels.jl +++ /dev/null @@ -1,11 +0,0 @@ -module SlabSeaIceModels - -# using RootSolvers: find_zero - -include("slab_heat_and_tracer_fluxes.jl") -include("slab_tendency_kernel_functions.jl") -include("slab_sea_ice_model.jl") -include("slab_time_stepping.jl") - -end # module - diff --git a/src/SlabSeaIceModels/slab_sea_ice_model.jl b/src/SlabSeaIceModels/slab_sea_ice_model.jl deleted file mode 100644 index 6dab911..0000000 --- a/src/SlabSeaIceModels/slab_sea_ice_model.jl +++ /dev/null @@ -1,181 +0,0 @@ -using ClimaSeaIce: - PhaseTransitions, - ForwardEulerTimestepper - -using ClimaSeaIce.HeatBoundaryConditions: - MeltingConstrainedFluxBalance, - IceWaterThermalEquilibrium, - PrescribedTemperature, - FluxFunction, - flux_summary - -using Oceananigans.Utils: prettysummary -using Oceananigans.TimeSteppers: Clock -using Oceananigans.Fields: field, Field, Center, ZeroField, ConstantField - -# Simulations interface -import Oceananigans: fields, prognostic_fields -import Oceananigans.Fields: set! -import Oceananigans.Models: AbstractModel -import Oceananigans.OutputWriters: default_included_properties -import Oceananigans.Simulations: reset!, initialize!, iteration -import Oceananigans.TimeSteppers: time_step!, update_state! -import Oceananigans.Utils: prettytime - -# TODO: move to Oceananigans -# import Oceananigans.Fields: field -# field(loc, a::Number, grid) = ConstantField(a) - -struct SlabSeaIceModel{GR, CL, TS, IT, IC, ST, IS, U, STF, TBC, CF, P, MIT, A} <: AbstractModel{TS} - grid :: GR - clock :: CL - timestepper :: TS - # State - ice_thickness :: IT - ice_concentration :: IC - top_surface_temperature :: ST - ice_salinity :: IS - velocities :: U - # Boundary conditions - external_heat_fluxes :: STF - heat_boundary_conditions :: TBC - # Internal flux - internal_heat_flux :: CF - # Melting and freezing stuff - phase_transitions :: P - ice_consolidation_thickness :: MIT - # Numerics - advection :: A -end - -const SSIM = SlabSeaIceModel - -Base.summary(model::SSIM) = "SlabSeaIceModel" -prettytime(model::SSIM) = prettytime(model.clock.time) -iteration(model::SSIM) = model.clock.iteration - -function Base.show(io::IO, model::SSIM) - grid = model.grid - arch = architecture(grid) - gridname = typeof(grid).name.wrapper - timestr = string("(time = ", prettytime(model), ", iteration = ", iteration(model), ")") - - print(io, "SlabSeaIceModel{", typeof(arch), ", ", gridname, "}", timestr, '\n') - print(io, "├── grid: ", summary(model.grid), '\n') - print(io, "├── top_surface_temperature: ", summary(model.top_surface_temperature), '\n') - print(io, "├── minimium_ice_thickness: ", prettysummary(model.ice_consolidation_thickness), '\n') - print(io, "└── external_heat_fluxes: ", '\n') - print(io, " ├── top: ", flux_summary(model.external_heat_fluxes.top, " │"), '\n') - print(io, " └── bottom: ", flux_summary(model.external_heat_fluxes.bottom, " ")) -end - -reset!(::SSIM) = nothing -initialize!(::SSIM) = nothing -default_included_properties(::SSIM) = tuple(:grid) - -fields(model::SSIM) = (h = model.ice_thickness, - ℵ = model.ice_concentration, - Tᵤ = model.top_surface_temperature, - Sᵢ = model.ice_salinity) - -# TODO: make this correct -prognostic_fields(model::SSIM) = fields(model) - -""" - SlabSeaIceModel(grid; kw...) - -Pretty simple model for sea ice. -""" -function SlabSeaIceModel(grid; - clock = Clock{eltype(grid)}(time = 0), - ice_thickness = Field{Center, Center, Nothing}(grid), - ice_consolidation_thickness = 0.0, # m - ice_concentration = Field{Center, Center, Nothing}(grid), - ice_salinity = 0, # psu - top_surface_temperature = nothing, - top_heat_flux = nothing, - bottom_heat_flux = 0, - velocities = nothing, - advection = nothing, - top_heat_boundary_condition = MeltingConstrainedFluxBalance(), - bottom_heat_boundary_condition = IceWaterThermalEquilibrium(), - # Default internal flux: thermal conductivity of 2 kg m s⁻³ K⁻¹, appropriate for freshwater ice - internal_heat_flux = ConductiveFlux(eltype(grid), conductivity=2), - phase_transitions = PhaseTransitions(eltype(grid))) - - if isnothing(velocities) - velocities = (u = ZeroField(), v=ZeroField(), w=ZeroField()) - end - - # Only one time-stepper is supported currently - timestepper = ForwardEulerTimestepper() - FT = eltype(grid) - - # TODO: pass `clock` into `field`, so functions can be time-dependent? - # Wrap ice_salinity in a field - ice_salinity = field((Center, Center, Nothing), ice_salinity, grid) - ice_consolidation_thickness = field((Center, Center, Nothing), ice_consolidation_thickness, grid) - - # Construct an internal heat flux function that captures the liquidus and - # bottom boundary condition. - parameters = (flux = internal_heat_flux, - liquidus = phase_transitions.liquidus, - bottom_heat_boundary_condition = bottom_heat_boundary_condition) - - internal_heat_flux_function = FluxFunction(slab_internal_heat_flux; - parameters, - top_temperature_dependent=true) - - # Construct default top heat flux if one is not provided - if isnothing(top_heat_flux) - if top_heat_boundary_condition isa PrescribedTemperature - # Default: external top flux is in equilibrium with internal fluxes - top_heat_flux = internal_heat_flux_function - else - # Default: no external top surface flux - top_heat_flux = 0 - end - end - - # Construct default top temperature if one is not provided - if isnothing(top_surface_temperature) - # Check top boundary condition - if top_heat_boundary_condition isa PrescribedTemperature - top_surface_temperature = top_heat_boundary_condition.temperature - else # build the default - top_surface_temperature = Field{Center, Center, Nothing}(grid) - end - end - - # Convert to `field` (does nothing if it's already a Field) - top_surface_temperature = field((Center, Center, Nothing), top_surface_temperature, grid) - - # Package the external fluxes and boundary conditions - external_heat_fluxes = (top = top_heat_flux, - bottom = bottom_heat_flux) - - heat_boundary_conditions = (top = top_heat_boundary_condition, - bottom = bottom_heat_boundary_condition) - - return SlabSeaIceModel(grid, - clock, - timestepper, - ice_thickness, - ice_concentration, - top_surface_temperature, - ice_salinity, - velocities, - external_heat_fluxes, - heat_boundary_conditions, - internal_heat_flux_function, - phase_transitions, - ice_consolidation_thickness, - advection) -end - -function set!(model::SSIM; h=nothing, α=nothing) - !isnothing(h) && set!(model.ice_thickness, h) - !isnothing(α) && set!(model.ice_conentration, α) - return nothing -end - diff --git a/src/SlabSeaIceModels/slab_tendency_kernel_functions.jl b/src/SlabSeaIceModels/slab_tendency_kernel_functions.jl deleted file mode 100644 index 75501c7..0000000 --- a/src/SlabSeaIceModels/slab_tendency_kernel_functions.jl +++ /dev/null @@ -1,77 +0,0 @@ -using ClimaSeaIce: latent_heat -using Oceananigans.Advection - -# Thickness change due to accretion and melting, restricted by minimum allowable value -function ice_thickness_tendency(i, j, grid, clock, - velocities, - advection, - ice_thickness, - ice_consolidation_thickness, - top_temperature, - bottom_heat_bc, - top_external_heat_flux, - internal_heat_flux, - bottom_external_heat_flux, - phase_transitions, - h_forcing, - model_fields) - - Gh_advection = - div_Uc(i, j, 1, grid, advection, velocities, ice_thickness) - - @inbounds begin - hᶜ = ice_consolidation_thickness[i, j, 1] - hᵢ = ice_thickness[i, j, 1] - Tuᵢ = top_temperature[i, j, 1] - end - - # Consolidation criteria - consolidated_ice = hᵢ >= hᶜ - - liquidus = phase_transitions.liquidus - Tbᵢ = bottom_temperature(i, j, grid, bottom_heat_bc, liquidus) - ℰb = latent_heat(phase_transitions, Tbᵢ) - ℰu = latent_heat(phase_transitions, Tuᵢ) - - Qi = internal_heat_flux - Qb = bottom_external_heat_flux - Qu = top_external_heat_flux - - # Retrieve fluxes - Quᵢ = getflux(Qu, i, j, grid, Tuᵢ, clock, model_fields) - Qiᵢ = getflux(Qi, i, j, grid, Tuᵢ, clock, model_fields) - Qbᵢ = getflux(Qb, i, j, grid, Tuᵢ, clock, model_fields) - - # Compute forcing - Fh = zero(grid) #h_forcing(i, j, grid, clock, model_fields) - - # If ice is consolidated, compute tendency for an ice slab; otherwise - # just add ocean fluxes from frazil ice formation or melting - slushy_Gh = - Qbᵢ / ℰb + Fh - - # Upper (top) and bottom interface velocities - wu = (Quᵢ - Qiᵢ) / ℰu # < 0 => melting - wb = (Qiᵢ - Qbᵢ) / ℰb # < 0 => freezing - - slabby_Gh = wu + wb + Fh - - return Gh_advection + ifelse(consolidated_ice, slabby_Gh, slushy_Gh) -end - -#= -function tracer_tendency(i, j, grid, clock, - c_forcing, - c, - ice_thickness, - top_tracer_concentration, - bottom_tracer_concentration, - model_fields) - - wu = top_interface_velocity(i, j, grid, Tu, Qi, Qu, phase_transitions) - wb = bottom_interface_velocity(i, j, grid, Tu, Qi, Qb, phase_transitions) - - cu = top_tracer_concentration[i, j] - cb = bottom_tracer_concentration[i, j] - - return cu * wu - cb * wb + c_forcing(i, j, grid, clock, model_fields) -end -=# diff --git a/src/SlabSeaIceModels/slab_time_stepping.jl b/src/SlabSeaIceModels/slab_time_stepping.jl deleted file mode 100644 index 5b11588..0000000 --- a/src/SlabSeaIceModels/slab_time_stepping.jl +++ /dev/null @@ -1,118 +0,0 @@ -using ClimaSeaIce.HeatBoundaryConditions: - PrescribedTemperature, - bottom_temperature, - top_surface_temperature, - getflux - -using Oceananigans.Architectures: architecture -using Oceananigans.BoundaryConditions: fill_halo_regions! -using Oceananigans.TimeSteppers: tick! -using Oceananigans.Utils: launch! - -using KernelAbstractions: @index, @kernel - -function time_step!(model::SSIM, Δt; callbacks=nothing) - grid = model.grid - arch = architecture(grid) - - launch!(arch, grid, :xyz, - slab_model_time_step!, - model.ice_thickness, - Δt, - grid, - model.clock, - model.velocities, - model.advection, - model.ice_concentration, - model.top_surface_temperature, - model.heat_boundary_conditions.top, - model.heat_boundary_conditions.bottom, - model.external_heat_fluxes.top, - model.internal_heat_flux, - model.external_heat_fluxes.bottom, - model.ice_consolidation_thickness, - model.phase_transitions, - nothing, #model.forcing.h, - fields(model)) - - tick!(model.clock, Δt) - - return nothing -end - -@kernel function slab_model_time_step!(ice_thickness, Δt, - grid, - clock, - velocities, - advection, - ice_concentration, - top_temperature, - top_heat_bc, - bottom_heat_bc, - top_external_heat_flux, - internal_heat_flux, - bottom_external_heat_flux, - ice_consolidation_thickness, - phase_transitions, - h_forcing, - model_fields) - - i, j = @index(Global, NTuple) - - ℵ = ice_concentration - h = ice_thickness - hᶜ = ice_consolidation_thickness - Qi = internal_heat_flux - Qu = top_external_heat_flux - Qb = bottom_external_heat_flux - Tu = top_temperature - liquidus = phase_transitions.liquidus - - # Determine top surface temperature - if !isa(top_heat_bc, PrescribedTemperature) # update surface temperature? - - consolidated_ice = @inbounds h[i, j, 1] >= hᶜ[i, j, 1] - - if consolidated_ice # slab is consolidated and has an independent surface temperature - Tu⁻ = @inbounds Tu[i, j, 1] - Tuⁿ = top_surface_temperature(i, j, grid, top_heat_bc, Tu⁻, Qi, Qu, clock, model_fields) - else # slab is unconsolidated and does not have an independent surface temperature - Tuⁿ = bottom_temperature(i, j, grid, bottom_heat_bc, liquidus) - end - - @inbounds Tu[i, j, 1] = Tuⁿ - end - - Gh = ice_thickness_tendency(i, j, grid, clock, - velocities, - advection, - ice_thickness, - ice_consolidation_thickness, - top_temperature, - bottom_heat_bc, - top_external_heat_flux, - internal_heat_flux, - bottom_external_heat_flux, - phase_transitions, - h_forcing, - model_fields) - - # Update ice thickness, clipping negative values - - @inbounds begin - h⁺ = h[i, j, 1] + Δt * Gh - h[i, j, 1] = max(zero(grid), h⁺) - - # Belongs in update state? - # That's certainly a simple model for ice concentration - consolidated_ice = h[i, j, 1] >= hᶜ[i, j, 1] - ℵ[i, j, 1] = ifelse(consolidated_ice, one(grid), zero(grid)) - end -end - -function update_state!(model::SSIM) - h = model.ice_thickness - fill_halo_regions!(h) - return nothing -end - diff --git a/src/sea_ice_model.jl b/src/sea_ice_model.jl new file mode 100644 index 0000000..d1c34fb --- /dev/null +++ b/src/sea_ice_model.jl @@ -0,0 +1,113 @@ +using Oceananigans.Fields: TracerFields +using ClimaSeaIce.SeaIceThermodynamics: external_top_heat_flux +using Oceananigans: tupleit +using ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions: flux_summary + +struct SeaIceModel{GR, TD, CL, TS, U, T, IT, IC, STF, SMS, A} <: AbstractModel{TS} + grid :: GR + clock :: CL + timestepper :: TS + # Prognostic State + velocities :: U + tracers :: T + ice_thickness :: IT + ice_concentration :: IC + # Thermodynamics + ice_thermodynamics :: TD + # External boundary conditions + external_heat_fluxes :: STF + external_momentum_stresses :: SMS + # Numerics + advection :: A +end + +function SeaIceModel(grid; + clock = Clock{eltype(grid)}(time = 0), + ice_thickness = Field{Center, Center, Nothing}(grid), + ice_concentration = Field{Center, Center, Nothing}(grid), + ice_salinity = 0, # psu + top_heat_flux = nothing, + bottom_heat_flux = 0, + velocities = nothing, + advection = nothing, + top_momentum_stress = nothing, # Fix when introducing dynamics + tracers = (), + boundary_conditions = NamedTuple(), + ice_thermodynamics = SlabSeaIceThermodynamics(grid)) + + if isnothing(velocities) + velocities = (u = ZeroField(), v=ZeroField(), w=ZeroField()) + end + + # Only one time-stepper is supported currently + timestepper = ForwardEulerTimestepper() + + tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example) + tracers = TracerFields(tracers, grid, boundary_conditions) + + # TODO: pass `clock` into `field`, so functions can be time-dependent? + # Wrap ice_salinity in a field + ice_salinity = field((Center, Center, Nothing), ice_salinity, grid) + tracers = merge(tracers, (; S = ice_salinity)) + + top_heat_flux = external_top_heat_flux(ice_thermodynamics, top_heat_flux) + + # Package the external fluxes and boundary conditions + external_heat_fluxes = (top = top_heat_flux, + bottom = bottom_heat_flux) + + return SeaIceModel(grid, + clock, + timestepper, + velocities, + tracers, + ice_thickness, + ice_concentration, + ice_thermodynamics, + external_heat_fluxes, + top_momentum_stress, + advection) +end + +const SIM = SeaIceModel + +function set!(model::SIM; h=nothing, ℵ=nothing) + !isnothing(h) && set!(model.ice_thickness, h) + !isnothing(ℵ) && set!(model.ice_conentration, ℵ) + return nothing +end + +Base.summary(model::SIM) = "SeaIceModel" +prettytime(model::SIM) = prettytime(model.clock.time) +iteration(model::SIM) = model.clock.iteration + +function Base.show(io::IO, model::SIM) + grid = model.grid + arch = architecture(grid) + gridname = typeof(grid).name.wrapper + timestr = string("(time = ", prettytime(model), ", iteration = ", iteration(model), ")") + + print(io, "SeaIceModel{", typeof(arch), ", ", gridname, "}", timestr, '\n') + print(io, "├── grid: ", summary(model.grid), '\n') + print(io, "├── ice thermodynamics: ", summary(model.ice_thermodynamics), '\n') + print(io, "├── advection: ", summary(model.advection), '\n') + print(io, "└── external_heat_fluxes: ", '\n') + print(io, " ├── top: ", flux_summary(model.external_heat_fluxes.top, " │"), '\n') + print(io, " └── bottom: ", flux_summary(model.external_heat_fluxes.bottom, " ")) +end + +reset!(::SIM) = nothing +initialize!(::SIM) = nothing +default_included_properties(::SIM) = tuple(:grid) + +fields(model::SIM) = merge((; h = model.ice_thickness, + ℵ = model.ice_concentration), + model.tracers, + model.velocities, + fields(model.ice_thermodynamics)) + +# TODO: make this correct +prognostic_fields(model::SIM) = merge((; h = model.ice_thickness, + ℵ = model.ice_concentration), + model.tracers, + model.velocities) diff --git a/src/time_stepping.jl b/src/time_stepping.jl new file mode 100644 index 0000000..18fd29a --- /dev/null +++ b/src/time_stepping.jl @@ -0,0 +1,75 @@ +using Oceananigans.Architectures: architecture +using Oceananigans.BoundaryConditions: fill_halo_regions! +using Oceananigans.TimeSteppers: tick! +using Oceananigans.Utils: launch! + +using KernelAbstractions: @index, @kernel + +function time_step!(model::SIM, Δt; callbacks=nothing) + grid = model.grid + arch = architecture(grid) + + launch!(arch, grid, :xyz, + sea_ice_time_step!, + model.ice_thickness, + Δt, + grid, + model.clock, + model.velocities, + model.advection, + model.ice_concentration, + model.ice_thermodynamics, + model.external_heat_fluxes.top, + model.external_heat_fluxes.bottom, + nothing, #model.forcing.h, + fields(model)) + + tick!(model.clock, Δt) + + return nothing +end + +@kernel function sea_ice_time_step!(ice_thickness, Δt, + grid, + clock, + velocities, + advection, + concentration, + thermodynamics, + top_external_heat_flux, + bottom_external_heat_flux, + h_forcing, + model_fields) + + i, j, k = @index(Global, NTuple) + h = ice_thickness + + Gh = ice_thickness_tendency(i, j, k, grid, clock, + velocities, + advection, + ice_thickness, + concentration, + thermodynamics, + top_external_heat_flux, + bottom_external_heat_flux, + h_forcing, + model_fields) + + # Update ice thickness, clipping negative values + + @inbounds begin + h⁺ = h[i, j, k] + Δt * Gh + h[i, j, k] = max(zero(grid), h⁺) + # Belongs in update state? + # That's certainly a simple model for ice concentration + # consolidated_ice = h[i, j, 1] >= hᶜ[i, j, 1] + # ℵ[i, j, 1] = ifelse(consolidated_ice, one(grid), zero(grid)) + end +end + +function update_state!(model::SIM) + h = model.ice_thickness + fill_halo_regions!(h) + return nothing +end + diff --git a/src/tracer_tendency_kernel_functions.jl b/src/tracer_tendency_kernel_functions.jl new file mode 100644 index 0000000..5c53ea0 --- /dev/null +++ b/src/tracer_tendency_kernel_functions.jl @@ -0,0 +1,31 @@ +using Oceananigans.Advection +using ClimaSeaIce.SeaIceThermodynamics: thickness_thermodynamic_tendency + +# Thickness change due to accretion and melting, restricted by minimum allowable value +function ice_thickness_tendency(i, j, k, grid, clock, + velocities, + advection, + ice_thickness, + ice_concentration, + thermodynamics, + top_external_heat_flux, + bottom_external_heat_flux, + h_forcing, + model_fields) + + Gh_advection = - div_Uc(i, j, k, grid, advection, velocities, ice_thickness) + + Gh_growth = thickness_thermodynamic_tendency(i, j, k, grid, + ice_thickness, + ice_concentration, + thermodynamics, + top_external_heat_flux, + bottom_external_heat_flux, + clock, model_fields) + + + # Compute forcing + Fh = zero(grid) #h_forcing(i, j, grid, clock, model_fields) + + return Gh_advection + Gh_growth + Fh +end diff --git a/validation/ice_ocean_model/cooling_then_warming_ocean.jl b/validation/ice_ocean_model/cooling_then_warming_ocean.jl index 53856c9..a457940 100644 --- a/validation/ice_ocean_model/cooling_then_warming_ocean.jl +++ b/validation/ice_ocean_model/cooling_then_warming_ocean.jl @@ -3,7 +3,7 @@ using Oceananigans.Utils: prettysummary using Oceananigans.Units using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity using ClimaSeaIce -using ClimaSeaIce: melting_temperature +using ClimaSeaIce.SeaIceThermodynamics: melting_temperature using SeawaterPolynomials: TEOS10EquationOfState using ClimaSeaIce.HeatBoundaryConditions: RadiativeEmission, IceWaterThermalEquilibrium using GLMakie