diff --git a/.vscode/settings.json b/.vscode/settings.json index 96873de6d..a1f38a012 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -69,6 +69,7 @@ ".pixi" ], "julia.lint.run": true, + "julia.NumThreads": "1", "julia.numTestProcesses": 6, "notebook.codeActionsOnSave": { "notebook.source.fixAll": "explicit", diff --git a/Manifest.toml b/Manifest.toml index 6629d3163..865db5c8d 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,12 +2,12 @@ julia_version = "1.12.6" manifest_format = "2.0" -project_hash = "bb9d32b9745e856908beecf2f9935e2cb46d11ff" +project_hash = "125da60ef05fc7cd4f6850a7e5233b46fb684ce5" [[deps.ADTypes]] -git-tree-sha1 = "f7304359109c768cf32dc5fa2d371565bb63b68a" +git-tree-sha1 = "bbc22a9a08a0ef6460041086d8a7b27940ed4ffd" uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b" -version = "1.21.0" +version = "1.22.0" [deps.ADTypes.extensions] ADTypesChainRulesCoreExt = "ChainRulesCore" @@ -49,10 +49,10 @@ version = "0.1.44" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [[deps.Adapt]] -deps = ["LinearAlgebra", "Requires"] -git-tree-sha1 = "35ea197a51ce46fcd01c4a44befce0578a1aaeca" +deps = ["LinearAlgebra"] +git-tree-sha1 = "28e1637322d4019ed2577cbec9268fab9b7da117" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "4.5.0" +version = "4.6.0" weakdeps = ["SparseArrays", "StaticArrays"] [deps.Adapt.extensions] @@ -88,9 +88,9 @@ version = "0.4.0" [[deps.ArrayInterface]] deps = ["Adapt", "LinearAlgebra"] -git-tree-sha1 = "78b3a7a536b4b0a747a0f296ea77091ca0a9f9a3" +git-tree-sha1 = "3d0cabd25fab32390e3bcb82cd67e700aebd9816" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "7.23.0" +version = "7.25.0" [deps.ArrayInterface.extensions] ArrayInterfaceAMDGPUExt = "AMDGPU" @@ -135,23 +135,11 @@ git-tree-sha1 = "ca36dc6dc9bbeea3dbbe3901837496ece85276e5" uuid = "59605e27-edc0-445a-b93d-c09a3a50b330" version = "0.1.1" -[[deps.BenchmarkTools]] -deps = ["Compat", "JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"] -git-tree-sha1 = "6876e30dc02dc69f0613cb6ece242144f2ca9e56" -uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -version = "1.7.0" - [[deps.BitFlags]] git-tree-sha1 = "0691e34b3bb8be9307330f88d1a3c3f25466c24d" uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35" version = "0.1.9" -[[deps.BitTwiddlingConvenienceFunctions]] -deps = ["Static"] -git-tree-sha1 = "f21cfd4950cb9f0587d5067e69405ad2acd27b87" -uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b" -version = "0.1.6" - [[deps.Blosc_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Zlib_jll", "Zstd_jll"] git-tree-sha1 = "535c80f1c0847a4c967ea945fca21becc9de1522" @@ -160,9 +148,9 @@ version = "1.21.7+0" [[deps.BracketingNonlinearSolve]] deps = ["CommonSolve", "ConcreteStructs", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase"] -git-tree-sha1 = "d9b66401c1fa982c7ca984d0566af5a9b3551420" +git-tree-sha1 = "7ad7171d693ae5552ac43862e7f6b61df4471c2b" uuid = "70df07ce-3d50-431d-a3e7-ca6ddb60ac1e" -version = "1.12.0" +version = "1.12.1" [deps.BracketingNonlinearSolve.extensions] BracketingNonlinearSolveChainRulesCoreExt = ["ChainRulesCore", "ForwardDiff"] @@ -180,15 +168,9 @@ version = "1.0.9+0" [[deps.CFTime]] deps = ["Dates", "Printf"] -git-tree-sha1 = "b2c9f2d3b8014323c0a8f2b401b68fa6bb08ed06" +git-tree-sha1 = "912c24c352c4167df4ebdacb96a432a2e3dbaf7a" uuid = "179af706-886a-5703-950a-314cd64e0468" -version = "0.2.8" - -[[deps.CPUSummary]] -deps = ["CpuId", "IfElse", "PrecompileTools", "Preferences", "Static"] -git-tree-sha1 = "f3a21d7fc84ba618a779d1ed2fcca2e682865bab" -uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9" -version = "0.2.7" +version = "0.2.9" [[deps.CSV]] deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "PrecompileTools", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings", "WorkerUtilities"] @@ -197,16 +179,10 @@ uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" version = "0.10.16" [[deps.Cairo_jll]] -deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] -git-tree-sha1 = "a21c5464519504e41e0cbc91f0188e8ca23d7440" +deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] +git-tree-sha1 = "1fa950ebc3e37eccd51c6a8fe1f92f7d86263522" uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a" -version = "1.18.5+1" - -[[deps.CloseOpenIntervals]] -deps = ["Static", "StaticArrayInterface"] -git-tree-sha1 = "05ba0d07cd4fd8b7a39541e31a7b0254704ea581" -uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9" -version = "0.1.13" +version = "1.18.7+0" [[deps.CodeTracking]] deps = ["InteractiveUtils", "REPL", "UUIDs"] @@ -275,11 +251,6 @@ git-tree-sha1 = "cda2cfaebb4be89c9084adaca7dd7333369715c5" uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" version = "0.3.1" -[[deps.CommonWorldInvalidations]] -git-tree-sha1 = "ae52d1c52048455e85a387fbee9be553ec2b68d0" -uuid = "f70d9fcc-98c5-4d4a-abd7-e4cdeebd8ca8" -version = "1.0.0" - [[deps.Compat]] deps = ["TOML", "UUIDs"] git-tree-sha1 = "9d8a54ce4b17aa5bdce0ea5c34bc5e7c340d16ad" @@ -310,9 +281,9 @@ weakdeps = ["InverseFunctions"] CompositionsBaseInverseFunctionsExt = "InverseFunctions" [[deps.ConcreteStructs]] -git-tree-sha1 = "f749037478283d372048690eb3b5f92a79432b34" +git-tree-sha1 = "ed1da4eac5ba9b3f6d061c90f3ca6ba190dd6595" uuid = "2569d6c7-a4a2-43d3-a901-331e8e4be471" -version = "0.2.3" +version = "0.2.4" [[deps.ConcurrentUtilities]] deps = ["Serialization", "Sockets"] @@ -352,12 +323,6 @@ git-tree-sha1 = "439e35b0b36e2e5881738abc8857bd92ad6ff9a8" uuid = "d38c429a-6771-53c6-b99e-75d170b6e991" version = "0.6.3" -[[deps.CpuId]] -deps = ["Markdown"] -git-tree-sha1 = "fcbb72b032692610bfbdb15018ac16a36cf2e406" -uuid = "adafc99b-e345-5852-983c-f28acb93d879" -version = "0.3.1" - [[deps.Crayons]] git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15" uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" @@ -375,19 +340,20 @@ version = "1.16.0" [[deps.DataFrames]] deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] -git-tree-sha1 = "d8928e9169ff76c6281f39a659f9bca3a573f24c" +git-tree-sha1 = "5fab31e2e01e70ad66e3e24c968c264d1cf166d6" uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -version = "1.8.1" +version = "1.8.2" [[deps.DataInterpolations]] deps = ["EnumX", "FindFirstFunctions", "ForwardDiff", "LinearAlgebra", "PrettyTables", "RecipesBase", "Reexport"] -git-tree-sha1 = "db37d8739c369b9e7212f8e61e37611bda6fa2e1" +git-tree-sha1 = "4d50e83772222ce023efb7d30632198828ea56a6" uuid = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" -version = "8.9.0" +version = "8.10.0" [deps.DataInterpolations.extensions] DataInterpolationsChainRulesCoreExt = "ChainRulesCore" DataInterpolationsMakieExt = "Makie" + DataInterpolationsMooncakeExt = ["Mooncake", "ChainRulesCore"] DataInterpolationsOptimExt = "Optim" DataInterpolationsRegularizationToolsExt = "RegularizationTools" DataInterpolationsSparseConnectivityTracerExt = ["SparseConnectivityTracer", "FillArrays"] @@ -397,6 +363,7 @@ version = "8.9.0" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" + Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6" Optim = "429524aa-4258-5aef-a3af-852621145aeb" RegularizationTools = "29dad682-9a27-4bc3-9c72-016788665182" SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" @@ -432,10 +399,10 @@ uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" version = "1.9.1" [[deps.DiffEqBase]] -deps = ["ArrayInterface", "BracketingNonlinearSolve", "ConcreteStructs", "DocStringExtensions", "FastBroadcast", "FastClosures", "FastPower", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "Setfield", "Static", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"] -git-tree-sha1 = "c5fe5125fcba8f98cdc5c7221b6c324883899c07" +deps = ["ArrayInterface", "BracketingNonlinearSolve", "ConcreteStructs", "DocStringExtensions", "FastBroadcast", "FastClosures", "FastPower", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLLogging", "SciMLOperators", "SciMLStructures", "Setfield", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"] +git-tree-sha1 = "86932ebaae0425f3d69c398dd37d57745b7cf524" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "6.213.0" +version = "7.5.0" [deps.DiffEqBase.extensions] DiffEqBaseCUDAExt = "CUDA" @@ -476,9 +443,9 @@ version = "6.213.0" [[deps.DiffEqCallbacks]] deps = ["ConcreteStructs", "DataStructures", "DiffEqBase", "DifferentiationInterface", "LinearAlgebra", "Markdown", "PrecompileTools", "RecipesBase", "RecursiveArrayTools", "SciMLBase", "StaticArraysCore"] -git-tree-sha1 = "f17b863c2d5d496363fe36c8d8535cc6a33c9952" +git-tree-sha1 = "51f445c25b80c26c83c61520cde214adf8c05227" uuid = "459566f4-90b8-5000-8ac3-15dfb0a30def" -version = "4.12.0" +version = "4.17.0" [deps.DiffEqCallbacks.extensions] DiffEqCallbacksFunctorsExt = "Functors" @@ -500,9 +467,9 @@ version = "1.15.1" [[deps.DifferentiationInterface]] deps = ["ADTypes", "LinearAlgebra"] -git-tree-sha1 = "7ae99144ea44715402c6c882bfef2adbeadbc4ce" +git-tree-sha1 = "2147a95a217cc8a78ec96ee03581adf129468e49" uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" -version = "0.7.16" +version = "0.7.18" [deps.DifferentiationInterface.extensions] DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore" @@ -512,8 +479,9 @@ version = "0.7.16" DifferentiationInterfaceFiniteDiffExt = "FiniteDiff" DifferentiationInterfaceFiniteDifferencesExt = "FiniteDifferences" DifferentiationInterfaceForwardDiffExt = ["ForwardDiff", "DiffResults"] - DifferentiationInterfaceGPUArraysCoreExt = "GPUArraysCore" + DifferentiationInterfaceGPUArraysCoreExt = ["GPUArraysCore", "Adapt"] DifferentiationInterfaceGTPSAExt = "GTPSA" + DifferentiationInterfaceHyperHessiansExt = "HyperHessians" DifferentiationInterfaceMooncakeExt = "Mooncake" DifferentiationInterfacePolyesterForwardDiffExt = ["PolyesterForwardDiff", "ForwardDiff", "DiffResults"] DifferentiationInterfaceReverseDiffExt = ["ReverseDiff", "DiffResults"] @@ -526,6 +494,7 @@ version = "0.7.16" DifferentiationInterfaceZygoteExt = ["Zygote", "ForwardDiff"] [deps.DifferentiationInterface.weakdeps] + Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" Diffractor = "9f5e2b26-1114-432f-b630-d3fe2085c51c" @@ -537,6 +506,7 @@ version = "0.7.16" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" GTPSA = "b27dd330-f138-47c5-815b-40db9dd9b6e8" + HyperHessians = "06b494a0-c8e0-40cc-ad32-d99506a00a6c" Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6" PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b" ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" @@ -576,9 +546,9 @@ version = "1.7.0" [[deps.Dualization]] deps = ["LinearAlgebra", "MathOptInterface"] -git-tree-sha1 = "060e90e82b8883742dbb3ce7b15c4efb24937fdc" +git-tree-sha1 = "d95cebd172f983fe6680a612098b6afa0f443640" uuid = "191a621a-6537-11e9-281d-650236a99e60" -version = "0.6.1" +version = "0.7.5" weakdeps = ["JuMP"] [deps.Dualization.extensions] @@ -590,9 +560,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56" version = "1.0.7" [[deps.EnzymeCore]] -git-tree-sha1 = "24bbb6fc8fb87eb71c1f8d00184a60fc22c63903" +git-tree-sha1 = "c6ee69ee502060982d12dbaaf3d8fcb4e835a0d1" uuid = "f151be2c-9106-41f4-ab19-57ee4f262869" -version = "0.8.19" +version = "0.8.20" [deps.EnzymeCore.extensions] AdaptExt = "Adapt" @@ -616,9 +586,9 @@ version = "0.1.11" [[deps.Expat_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "27af30de8b5445644e8ffe3bcb0d72049c089cf1" +git-tree-sha1 = "8f05e9a2e7c2e3eb524102bb2926c5743c07fbe1" uuid = "2e619515-83b5-522b-bb60-26c02a35a201" -version = "2.7.3+0" +version = "2.8.0+0" [[deps.ExprTools]] git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec" @@ -638,15 +608,23 @@ version = "0.4.5" [[deps.FFMPEG_jll]] deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libva_jll", "libvorbis_jll", "x264_jll", "x265_jll"] -git-tree-sha1 = "66381d7059b5f3f6162f28831854008040a4e905" +git-tree-sha1 = "cac41ca6b2d399adfc95e51240566f8a60a80806" uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5" -version = "8.0.1+1" +version = "8.1.0+0" [[deps.FastBroadcast]] -deps = ["ArrayInterface", "LinearAlgebra", "Polyester", "Static", "StaticArrayInterface", "StrideArraysCore"] -git-tree-sha1 = "ab1b34570bcdf272899062e1a56285a53ecaae08" +deps = ["ArrayInterface", "LinearAlgebra"] +git-tree-sha1 = "7feeed2c9a7fa272189a5561bebf0c4ccaedb6ec" uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898" -version = "0.3.5" +version = "1.3.2" + + [deps.FastBroadcast.extensions] + FastBroadcastPolyesterExt = "Polyester" + FastBroadcastStaticExt = "Static" + + [deps.FastBroadcast.weakdeps] + Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" + Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" [[deps.FastClosures]] git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef" @@ -717,9 +695,9 @@ version = "1.8.0" [[deps.FiniteDiff]] deps = ["ArrayInterface", "LinearAlgebra", "Setfield"] -git-tree-sha1 = "9340ca07ca27093ff68418b7558ca37b05f8aeb1" +git-tree-sha1 = "f7017a4f337f8df189fcce98e32b67a1298a2115" uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" -version = "2.29.0" +version = "2.31.0" [deps.FiniteDiff.extensions] FiniteDiffBandedMatricesExt = "BandedMatrices" @@ -778,10 +756,19 @@ uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" version = "1.1.3" [[deps.FunctionWrappersWrappers]] -deps = ["FunctionWrappers"] -git-tree-sha1 = "b104d487b34566608f8b4e1c39fb0b10aa279ff8" +deps = ["FunctionWrappers", "PrecompileTools", "TruncatedStacktraces"] +git-tree-sha1 = "c1b0c3a166a2a393257aa888787ca817532e14ce" uuid = "77dc65aa-8811-40c2-897b-53d922fa7daf" -version = "0.1.3" +version = "1.8.0" + + [deps.FunctionWrappersWrappers.extensions] + FunctionWrappersWrappersEnzymeExt = ["Enzyme", "EnzymeCore"] + FunctionWrappersWrappersMooncakeExt = "Mooncake" + + [deps.FunctionWrappersWrappers.weakdeps] + Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" + EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" + Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6" [[deps.Future]] deps = ["Random"] @@ -790,9 +777,9 @@ version = "1.11.0" [[deps.GLFW_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll", "libdecor_jll", "xkbcommon_jll"] -git-tree-sha1 = "b7bfd56fa66616138dfe5237da4dc13bbd83c67f" +git-tree-sha1 = "9e0fb9e54594c47f278d75063980e43066e26e20" uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89" -version = "3.4.1+0" +version = "3.4.1+1" [[deps.GPUArraysCore]] deps = ["Adapt"] @@ -837,9 +824,9 @@ uuid = "7746bdde-850d-59dc-9ae8-88ece973131d" version = "2.86.3+0" [[deps.Glob]] -git-tree-sha1 = "83cb0092e2792b9e3a865b6655e88f5b862607e2" +git-tree-sha1 = "246c628cec062230b7d183aab88841fa94fcabe9" uuid = "c27321d9-0574-5035-807b-f59d2c89b15c" -version = "1.4.0" +version = "1.5.0" [[deps.Graphite2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -866,10 +853,10 @@ uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe" version = "1.0.2" [[deps.HDF5_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"] -git-tree-sha1 = "e94f84da9af7ce9c6be049e9067e511e17ff89ec" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "MPIABI_jll", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "aws_c_s3_jll", "dlfcn_win32_jll", "libaec_jll", "mpif_jll"] +git-tree-sha1 = "45337643a2d97262d5fe72ce1f13e8a662d13d62" uuid = "0234f1f7-429e-5d53-9886-15a909be8d59" -version = "1.14.6+0" +version = "2.1.2+0" [[deps.HTTP]] deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "PrecompileTools", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] @@ -885,15 +872,15 @@ version = "8.5.1+0" [[deps.HiGHS]] deps = ["HiGHS_jll", "LinearAlgebra", "MathOptIIS", "MathOptInterface", "OpenBLAS32_jll", "PrecompileTools", "SparseArrays"] -git-tree-sha1 = "2b592162e8b40a5cf6435c37d86a8e2cf25cd459" +git-tree-sha1 = "bf5e946f72ebd1b4620249a6be7ff34832ba9ca0" uuid = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" -version = "1.22.2" +version = "1.23.0" [[deps.HiGHS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Zlib_jll", "libblastrampoline_jll"] -git-tree-sha1 = "621d773f277b9eadac7e049eaa6418af65c7b9d7" +git-tree-sha1 = "50ed12dc8c37ebb8d2b759f21755259d8512f2bd" uuid = "8fd58aa0-07eb-5a78-9b36-339c94fd15ea" -version = "1.13.1+0" +version = "1.14.0+0" [[deps.Hwloc_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "XML2_jll", "Xorg_libpciaccess_jll"] @@ -907,16 +894,11 @@ git-tree-sha1 = "0ee181ec08df7d7c911901ea38baf16f755114dc" uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" version = "1.0.0" -[[deps.IfElse]] -git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" -uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" -version = "0.1.1" - [[deps.Infiltrator]] deps = ["InteractiveUtils", "Markdown", "REPL", "UUIDs"] -git-tree-sha1 = "b1eccca4ebc42ac306496c5a5677e50251a19ab1" +git-tree-sha1 = "0330ef9ac27dc069bf97fa5ab1a4391be4fcc4f6" uuid = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" -version = "1.9.10" +version = "1.9.11" [[deps.Inflate]] git-tree-sha1 = "d1b1b796e47d94588b3757fe84fbf65a5ec4a80d" @@ -985,9 +967,9 @@ version = "0.1.11" [[deps.JLLWrappers]] deps = ["Artifacts", "Preferences"] -git-tree-sha1 = "0533e564aae234aff59ab625543145446d8b6ec2" +git-tree-sha1 = "7204148362dafe5fe6a273f855b8ccbe4df8173e" uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.7.1" +version = "1.8.0" [[deps.JSON]] deps = ["Dates", "Mmap", "Parsers", "Unicode"] @@ -1015,15 +997,15 @@ version = "0.2.1" [[deps.JpegTurbo_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "b6893345fd6658c8e475d40155789f4860ac3b21" +git-tree-sha1 = "c0c9b76f3520863909825cbecdef58cd63de705a" uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8" -version = "3.1.4+0" +version = "3.1.5+0" [[deps.JuMP]] deps = ["LinearAlgebra", "MacroTools", "MathOptInterface", "MutableArithmetics", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays"] -git-tree-sha1 = "4091a1338a0e32766b11b9bd3fac247d34200c77" +git-tree-sha1 = "6941586d9cf3c0af718bc6e6250dcf24057d412e" uuid = "4076af6c-e467-56ae-b986-b466b2749572" -version = "1.30.0" +version = "1.30.1" [deps.JuMP.extensions] JuMPDimensionalDataExt = "DimensionalData" @@ -1064,9 +1046,9 @@ version = "3.100.3+0" [[deps.LERC_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "aaafe88dccbd957a8d82f7d05be9b69172e0cee3" +git-tree-sha1 = "17b94ecafcfa45e8360a4fc9ca6b583b049e4e37" uuid = "88015f11-f218-50d7-93a8-a6af411a945d" -version = "4.0.1+0" +version = "4.1.0+0" [[deps.LLVMOpenMP_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -1083,12 +1065,6 @@ weakdeps = ["Serialization"] [deps.LRUCache.extensions] SerializationExt = ["Serialization"] -[[deps.LZO_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "1c602b1127f4751facb671441ca72715cc95938a" -uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac" -version = "2.10.3+0" - [[deps.LaTeXStrings]] git-tree-sha1 = "dda21b8cbd6a6c40d9d02a73230f9d70fed6918c" uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" @@ -1112,12 +1088,6 @@ version = "0.16.10" SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" tectonic_jll = "d7dd28d6-a5e6-559c-9131-7eb760cdacc5" -[[deps.LayoutPointers]] -deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface"] -git-tree-sha1 = "a9eaadb366f5493a5654e843864c13d8b107548c" -uuid = "10f19ff3-798f-405d-979b-55457f8fc047" -version = "0.1.17" - [[deps.LazilyInitializedFields]] git-tree-sha1 = "0f2da712350b020bc3957f269c9caad516383ee0" uuid = "0e77f7df-68c5-4e49-93ce-4cd80f5598bf" @@ -1183,9 +1153,9 @@ version = "1.18.0+0" [[deps.Libmount_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "97bbca976196f2a1eb9607131cb108c69ec3f8a6" +git-tree-sha1 = "cc3ad4faf30015a3e8094c9b5b7f19e85bdf2386" uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9" -version = "2.41.3+0" +version = "2.42.0+0" [[deps.Libtiff_jll]] deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"] @@ -1195,9 +1165,9 @@ version = "4.7.2+0" [[deps.Libuuid_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "d0205286d9eceadc518742860bf23f703779a3d6" +git-tree-sha1 = "d620582b1f0cbe2c72dd1d5bd195a9ce73370ab1" uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" -version = "2.41.3+0" +version = "2.42.0+0" [[deps.LicenseCheck]] deps = ["licensecheck_jll"] @@ -1207,9 +1177,9 @@ version = "0.2.3" [[deps.LineSearch]] deps = ["ADTypes", "CommonSolve", "ConcreteStructs", "FastClosures", "LinearAlgebra", "MaybeInplace", "PrecompileTools", "SciMLBase", "SciMLJacobianOperators", "StaticArraysCore"] -git-tree-sha1 = "9f7253c0574b4b585c8909232adb890930da980a" +git-tree-sha1 = "fd58a77c92e7c8f1db25c9839127d52943a49349" uuid = "87fe0de2-c867-4266-b59a-2f0a94fc965b" -version = "0.1.6" +version = "0.1.9" [deps.LineSearch.extensions] LineSearchLineSearchesExt = "LineSearches" @@ -1224,9 +1194,9 @@ version = "1.12.0" [[deps.LinearSolve]] deps = ["ArrayInterface", "ConcreteStructs", "DocStringExtensions", "EnumX", "GPUArraysCore", "InteractiveUtils", "Krylov", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "OpenBLAS_jll", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLLogging", "SciMLOperators", "Setfield", "StaticArraysCore"] -git-tree-sha1 = "1ddad5f2b0717f71f1588b3519e2f7d80fc2ce65" +git-tree-sha1 = "97a6bf19ef32518268ba15f166fc030df43b619a" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -version = "3.69.0" +version = "3.79.0" [deps.LinearSolve.extensions] LinearSolveAMDGPUExt = "AMDGPU" @@ -1239,6 +1209,7 @@ version = "3.69.0" LinearSolveCUSOLVERRFExt = ["CUSOLVERRF", "SparseArrays"] LinearSolveChainRulesCoreExt = "ChainRulesCore" LinearSolveCliqueTreesExt = ["CliqueTrees", "SparseArrays"] + LinearSolveElementalExt = "Elemental" LinearSolveEnzymeExt = ["EnzymeCore", "SparseArrays"] LinearSolveFastAlmostBandedMatricesExt = "FastAlmostBandedMatrices" LinearSolveFastLapackInterfaceExt = "FastLapackInterface" @@ -1250,10 +1221,12 @@ version = "3.69.0" LinearSolveKrylovKitExt = "KrylovKit" LinearSolveMetalExt = "Metal" LinearSolveMooncakeExt = "Mooncake" - LinearSolvePETScExt = ["PETSc", "SparseArrays"] + LinearSolvePETScExt = ["PETSc", "SparseArrays", "SparseMatricesCSR"] + LinearSolvePETScMPIExt = ["PETSc", "PartitionedArrays", "SparseArrays", "SparseMatricesCSR"] LinearSolveParUExt = ["ParU_jll", "SparseArrays"] LinearSolvePardisoExt = ["Pardiso", "SparseArrays"] LinearSolveRecursiveFactorizationExt = "RecursiveFactorization" + LinearSolveSTRUMPACKExt = ["SparseArrays", "STRUMPACK_jll"] LinearSolveSparseArraysExt = "SparseArrays" LinearSolveSparspakExt = ["SparseArrays", "Sparspak"] @@ -1267,6 +1240,7 @@ version = "3.69.0" CUSOLVERRF = "a8cc9031-bad2-4722-94f5-40deabb4245c" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" + Elemental = "902c3f28-d1ec-5e7e-8399-a24c3845ee38" EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e" FastLapackInterface = "29a986be-02c6-4525-aec4-84b980013641" @@ -1282,8 +1256,11 @@ version = "3.69.0" PETSc = "ace2c81b-2b5f-4b1e-a30d-d662738edfe0" ParU_jll = "9e0b026c-e8ce-559c-a2c4-6a3d5c955bc9" Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" + PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" + STRUMPACK_jll = "86fbd0b9-476f-557c-b766-62c724b42d8c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac" blis_jll = "6136c539-28a5-5bf0-87cc-b183200dce32" @@ -1331,11 +1308,17 @@ git-tree-sha1 = "282cadc186e7b2ae0eeadbd7a4dffed4196ae2aa" uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" version = "2025.2.0+0" +[[deps.MPIABI_jll]] +deps = ["Artifacts", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] +git-tree-sha1 = "9be143b6045719e8fb019d2b3bc2aebad1184fef" +uuid = "b5ada748-db0f-5fc0-8972-9331c762740c" +version = "0.1.5+0" + [[deps.MPICH_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] -git-tree-sha1 = "9341048b9f723f2ae2a72a5269ac2f15f80534dc" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "Libdl", "MPIPreferences", "TOML"] +git-tree-sha1 = "07dbec8aab01696edc0151a401a6cdfe95b9b885" uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4" -version = "4.3.2+0" +version = "5.0.1+0" [[deps.MPIPreferences]] deps = ["Libdl", "Preferences"] @@ -1345,20 +1328,15 @@ version = "0.1.12" [[deps.MPItrampoline_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"] -git-tree-sha1 = "36c2d142e7d45fb98b5f83925213feb3292ca348" +git-tree-sha1 = "675df097f8eeb28998b2cfe3b25655af73d5f7df" uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748" -version = "5.5.5+0" +version = "5.5.6+0" [[deps.MacroTools]] git-tree-sha1 = "1e0228a030642014fe5cfe68c2c0a818f9e3f522" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" version = "0.5.16" -[[deps.ManualMemory]] -git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd" -uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667" -version = "0.1.8" - [[deps.Markdown]] deps = ["Base64", "JuliaSyntaxHighlighting", "StyledStrings"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" @@ -1372,9 +1350,9 @@ version = "1.1.0" [[deps.MathOptAnalyzer]] deps = ["Dualization", "LinearAlgebra", "MathOptIIS", "MathOptInterface", "Printf"] -git-tree-sha1 = "d85d8fc4196d5f990dadeae5950201fbadb8cbad" +git-tree-sha1 = "5c0427e61a7f4396f513cb1e72270866140e275e" uuid = "d1179b25-476b-425c-b826-c7787f0fff83" -version = "0.1.1" +version = "0.1.2" weakdeps = ["JuMP"] [deps.MathOptAnalyzer.extensions] @@ -1382,15 +1360,23 @@ weakdeps = ["JuMP"] [[deps.MathOptIIS]] deps = ["MathOptInterface"] -git-tree-sha1 = "12f262200198606174e28c59da634e9d9da839ed" +git-tree-sha1 = "3b3d69130d8ab8c39d5fa4d30e20a8e6428c9d37" uuid = "8c4f8055-bd93-4160-a86b-a0c04941dbff" -version = "0.1.2" +version = "0.2.0" [[deps.MathOptInterface]] -deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test"] -git-tree-sha1 = "ce739e3d8a21313ea418772edfc3b7b15a1dfc16" +deps = ["CodecBzip2", "CodecZlib", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test"] +git-tree-sha1 = "73939c06e863f8d68322106fdc2464f3443b5e1a" uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" -version = "1.50.1" +version = "1.51.0" + + [deps.MathOptInterface.extensions] + MathOptInterfaceBenchmarkToolsExt = "BenchmarkTools" + MathOptInterfaceCliqueTreesExt = "CliqueTrees" + + [deps.MathOptInterface.weakdeps] + BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" + CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" [[deps.MaybeInplace]] deps = ["ArrayInterface", "LinearAlgebra", "MacroTools"] @@ -1476,15 +1462,15 @@ version = "0.2.4" [[deps.MutableArithmetics]] deps = ["LinearAlgebra", "SparseArrays", "Test"] -git-tree-sha1 = "7c25249fc13a070f5ba433c50e21e22bb33c6fb0" +git-tree-sha1 = "dc5b2c4c111c46bc79ac4405eeb563523b39c004" uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" -version = "1.7.1" +version = "1.8.0" [[deps.NCDatasets]] deps = ["CFTime", "CommonDataModel", "DataStructures", "Dates", "DiskArrays", "NetCDF_jll", "NetworkOptions", "Printf"] -git-tree-sha1 = "791095699024db9dbb95bc48c0507cb0d1938369" +git-tree-sha1 = "5eb7747d10437f5acb2675c1f865ffa0353f3f9c" uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" -version = "0.14.14" +version = "0.14.15" [deps.NCDatasets.extensions] NCDatasetsMPIExt = "MPI" @@ -1499,20 +1485,20 @@ uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" version = "1.1.3" [[deps.NetCDF_jll]] -deps = ["Artifacts", "Blosc_jll", "Bzip2_jll", "HDF5_jll", "JLLWrappers", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "TOML", "XML2_jll", "Zlib_jll", "Zstd_jll", "libaec_jll", "libzip_jll"] -git-tree-sha1 = "d574803b6055116af212434460adf654ce98e345" +deps = ["Artifacts", "Blosc_jll", "Bzip2_jll", "HDF5_jll", "JLLWrappers", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPIABI_jll", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "TOML", "XML2_jll", "Zlib_jll", "Zstd_jll", "libaec_jll", "libzip_jll"] +git-tree-sha1 = "8a36db9b934b0e72583e624abc8f3b3d60554f2c" uuid = "7243133f-43d8-5620-bbf4-c2c921802cf3" -version = "401.900.300+0" +version = "401.1000.0+0" [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" version = "1.3.0" [[deps.NonlinearSolve]] -deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolveQuasiNewton", "NonlinearSolveSpectralMethods", "PrecompileTools", "Preferences", "Reexport", "SciMLBase", "SciMLLogging", "SimpleNonlinearSolve", "StaticArraysCore", "SymbolicIndexingInterface"] -git-tree-sha1 = "d27bcf0cebf8786edcc2eaa4455c959e680334e7" +deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolveQuasiNewton", "NonlinearSolveSpectralMethods", "PrecompileTools", "Preferences", "Reexport", "SciMLBase", "Setfield", "SimpleNonlinearSolve", "StaticArraysCore", "SymbolicIndexingInterface"] +git-tree-sha1 = "a6c5719bbb42985c72f4cacbfa49e86bab850d66" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -version = "4.16.0" +version = "4.19.1" [deps.NonlinearSolve.extensions] NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt" @@ -1542,10 +1528,10 @@ version = "4.16.0" Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" [[deps.NonlinearSolveBase]] -deps = ["ADTypes", "Adapt", "ArrayInterface", "CommonSolve", "Compat", "ConcreteStructs", "DifferentiationInterface", "EnzymeCore", "FastClosures", "LinearAlgebra", "LogExpFunctions", "Markdown", "MaybeInplace", "PreallocationTools", "PrecompileTools", "Preferences", "Printf", "RecursiveArrayTools", "SciMLBase", "SciMLJacobianOperators", "SciMLLogging", "SciMLOperators", "SciMLStructures", "Setfield", "StaticArraysCore", "SymbolicIndexingInterface", "TimerOutputs"] -git-tree-sha1 = "a89529d343dbb09670a24df090787dc3475fba5d" +deps = ["ADTypes", "Adapt", "ArrayInterface", "CommonSolve", "Compat", "ConcreteStructs", "DifferentiationInterface", "EnzymeCore", "FastClosures", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "LogExpFunctions", "Markdown", "MaybeInplace", "PreallocationTools", "PrecompileTools", "Preferences", "Printf", "RecursiveArrayTools", "SciMLBase", "SciMLJacobianOperators", "SciMLLogging", "SciMLOperators", "SciMLStructures", "Setfield", "StaticArraysCore", "SymbolicIndexingInterface", "TimerOutputs"] +git-tree-sha1 = "36f063516ea7b8c355ec94819b409b89db436358" uuid = "be0214bd-f91f-a760-ac4e-3421ce2b2da0" -version = "2.19.0" +version = "2.26.0" [deps.NonlinearSolveBase.extensions] NonlinearSolveBaseBandedMatricesExt = "BandedMatrices" @@ -1575,15 +1561,15 @@ version = "2.19.0" [[deps.NonlinearSolveFirstOrder]] deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConcreteStructs", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "MaybeInplace", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase", "SciMLJacobianOperators", "Setfield", "StaticArraysCore"] -git-tree-sha1 = "eea7cbe389b168c77df7ff779fb7277019c685c8" +git-tree-sha1 = "ce68820a4f421fb5bee7ec4dcf875aff33886bfb" uuid = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d" -version = "2.0.0" +version = "2.1.1" [[deps.NonlinearSolveQuasiNewton]] deps = ["ArrayInterface", "CommonSolve", "ConcreteStructs", "LinearAlgebra", "LinearSolve", "MaybeInplace", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase", "SciMLOperators", "StaticArraysCore"] -git-tree-sha1 = "ade27e8e9566b6cec63ee62f6a6650a11cf9a2eb" +git-tree-sha1 = "538432ca1aea8bf63db02929bf870501f8a7c64c" uuid = "9a2c21bd-3a47-402d-9113-8faf9a0ee114" -version = "1.12.0" +version = "1.13.1" weakdeps = ["ForwardDiff"] [deps.NonlinearSolveQuasiNewton.extensions] @@ -1591,9 +1577,9 @@ weakdeps = ["ForwardDiff"] [[deps.NonlinearSolveSpectralMethods]] deps = ["CommonSolve", "ConcreteStructs", "LineSearch", "MaybeInplace", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase"] -git-tree-sha1 = "eafd027b5cd768f19bb5de76c0e908a9065ddd36" +git-tree-sha1 = "a3781e12becdf0ce5520bd97ec617e879bf4e9f2" uuid = "26075421-4e9a-44e1-8bd1-420ed7ad02b2" -version = "1.6.0" +version = "1.7.1" weakdeps = ["ForwardDiff"] [deps.NonlinearSolveSpectralMethods.extensions] @@ -1621,10 +1607,10 @@ uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051" version = "1.3.6+0" [[deps.OpenBLAS32_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] -git-tree-sha1 = "46cce8b42186882811da4ce1f4c7208b02deb716" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "libblastrampoline_jll"] +git-tree-sha1 = "565175ce692c065e50ad32efbb61ba69b1586593" uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2" -version = "0.3.30+0" +version = "0.3.33+1" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] @@ -1638,9 +1624,9 @@ version = "0.8.7+0" [[deps.OpenMPI_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML", "Zlib_jll"] -git-tree-sha1 = "2f3d05e419b6125ffe06e55784102e99325bdbe2" +git-tree-sha1 = "6d6c0ca4824268c1a7dca1f4721c535ac63d9074" uuid = "fe0851c0-eecd-5654-98d4-656369965a5c" -version = "5.0.10+0" +version = "5.0.11+0" [[deps.OpenSSL]] deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "NetworkOptions", "OpenSSL_jll", "Sockets"] @@ -1671,64 +1657,74 @@ uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" version = "1.8.1" [[deps.OrdinaryDiffEqBDF]] -deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqSDIRK", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "StaticArrays", "TruncatedStacktraces"] -git-tree-sha1 = "4758206b7578c8f884915ad589e6bd18548e16f3" +deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqSDIRK", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"] +git-tree-sha1 = "91e905c2224aa865ec545c5e7fd53c9a1020afa9" +repo-rev = "post-newton" +repo-subdir = "lib/OrdinaryDiffEqBDF" +repo-url = "https://github.com/visr/OrdinaryDiffEq.jl" uuid = "6ad6398a-0878-4a85-9266-38940aa047c8" -version = "1.23.0" +version = "2.1.1" [[deps.OrdinaryDiffEqCore]] -deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "ConcreteStructs", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "FastPower", "FunctionWrappersWrappers", "InteractiveUtils", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "Polyester", "PrecompileTools", "Preferences", "Random", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLLogging", "SciMLOperators", "SciMLStructures", "Static", "StaticArrayInterface", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"] -git-tree-sha1 = "14566414dda89ecb60ebd2a312042b440325d477" +deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "ConcreteStructs", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "FastPower", "FunctionWrappersWrappers", "InteractiveUtils", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "PrecompileTools", "Preferences", "Random", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLLogging", "SciMLOperators", "SciMLStructures", "SymbolicIndexingInterface", "TruncatedStacktraces"] +git-tree-sha1 = "a520345feb7ad2374786d748fad064cb3eef2a29" uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" -version = "3.25.0" +version = "4.2.1" [deps.OrdinaryDiffEqCore.extensions] OrdinaryDiffEqCoreMooncakeExt = "Mooncake" + OrdinaryDiffEqCorePolyesterExt = "Polyester" OrdinaryDiffEqCoreSparseArraysExt = "SparseArrays" [deps.OrdinaryDiffEqCore.weakdeps] Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6" + Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[deps.OrdinaryDiffEqDifferentiation]] -deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqCore", "SciMLBase", "SciMLOperators", "SparseMatrixColorings", "StaticArrayInterface", "StaticArrays"] -git-tree-sha1 = "dc6330b1155eeb609208b564331e8df9b5801bc6" +deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqCore", "SciMLBase", "SciMLOperators", "SparseMatrixColorings", "StaticArraysCore"] +git-tree-sha1 = "d9e7e8b9a6f7edf8c8ce9bba1f2cc9f8ad3a6752" uuid = "4302a76b-040a-498a-8c04-15b101fed76b" -version = "2.4.0" +version = "3.1.1" weakdeps = ["SparseArrays"] [deps.OrdinaryDiffEqDifferentiation.extensions] OrdinaryDiffEqDifferentiationSparseArraysExt = "SparseArrays" [[deps.OrdinaryDiffEqLowOrderRK]] -deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"] -git-tree-sha1 = "b01ab4362ad63d094dcde4c65672518212042dbc" +deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "SciMLBase"] +git-tree-sha1 = "f5eda8e7f4eb277d7e95566f49433f63864d2606" uuid = "1344f307-1e59-4825-a18e-ace9aa3fa4c6" -version = "1.11.0" +version = "2.1.0" [[deps.OrdinaryDiffEqNonlinearSolve]] -deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "FastClosures", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MuladdMacro", "NonlinearSolve", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "PreallocationTools", "RecursiveArrayTools", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleNonlinearSolve", "SparseArrays", "StaticArrays"] -git-tree-sha1 = "08dd3dfacb900ff27ffd2671f9b6207b7d929ab7" +deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "FastClosures", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MuladdMacro", "NonlinearSolve", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "PreallocationTools", "RecursiveArrayTools", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleNonlinearSolve", "SparseArrays", "StaticArraysCore"] +git-tree-sha1 = "c3e92bdb7bf16385e6f7b505f78601a2708eff28" uuid = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" -version = "1.24.0" +version = "2.0.0" [[deps.OrdinaryDiffEqRosenbrock]] -deps = ["ADTypes", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"] -git-tree-sha1 = "e3e1878da8f47becbc779d525717d964527ba1aa" +deps = ["ADTypes", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqRosenbrockTableaus", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase"] +git-tree-sha1 = "4b2fca93815c1b15eb9cb32cacd782b487654ea4" uuid = "43230ef6-c299-4910-a778-202eb28ce4ce" -version = "1.27.0" +version = "2.2.0" + +[[deps.OrdinaryDiffEqRosenbrockTableaus]] +git-tree-sha1 = "66acd57301ea9c79ff92b13b72135af258c141ab" +uuid = "b4bd8bb3-f80f-41d2-9b21-73a655b304b9" +version = "2.0.0" [[deps.OrdinaryDiffEqSDIRK]] deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"] -git-tree-sha1 = "ccc883881d5c74e08b28b299bb135c38d3496381" +git-tree-sha1 = "90e99f8d9970ea73fb16f7f87ae91fa30b2237d5" uuid = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" -version = "1.13.0" +version = "2.3.0" [[deps.OrdinaryDiffEqTsit5]] -deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static", "TruncatedStacktraces"] -git-tree-sha1 = "6e6173f4d0db8b8136d5fccdcdb7e8abe800a352" +deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"] +git-tree-sha1 = "d9cedcfa2cdae8859b4510b9d01ada9852cd6d9d" uuid = "b1df2697-797e-41e3-8120-5422d3b24e4a" -version = "1.10.0" +version = "2.0.0" [[deps.OteraEngine]] deps = ["Markdown", "TOML"] @@ -1749,15 +1745,15 @@ version = "2.3.0" [[deps.Pango_jll]] deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "Glib_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl"] -git-tree-sha1 = "0662b083e11420952f2e62e17eddae7fc07d5997" +git-tree-sha1 = "58e5ed5e386e156bd93e86b305ebd21ac63d2d04" uuid = "36c8627f-9965-5494-a995-c6b170f724f3" -version = "1.57.0+0" +version = "1.57.1+0" [[deps.Parsers]] deps = ["Dates", "PrecompileTools", "UUIDs"] -git-tree-sha1 = "7d2f8f21da5db6a806faf7b9b292296da42b2810" +git-tree-sha1 = "5d5e0a78e971354b1c7bff0655d11fdc1b0e12c8" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.8.3" +version = "2.8.4" [[deps.Patchelf_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] @@ -1773,9 +1769,9 @@ version = "1.3.0" [[deps.Pixman_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"] -git-tree-sha1 = "db76b1ecd5e9715f3d043cec13b2ec93ce015d53" +git-tree-sha1 = "e4a6721aa89e62e5d4217c0b21bd714263779dda" uuid = "30392449-352a-5448-841d-b1acce4e97dc" -version = "0.44.2+0" +version = "0.46.4+0" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"] @@ -1824,24 +1820,22 @@ version = "1.41.6" ImageInTerminal = "d8c32880-2388-543b-8c61-d9f865259254" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" -[[deps.Polyester]] -deps = ["ArrayInterface", "BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "ManualMemory", "PolyesterWeave", "Static", "StaticArrayInterface", "StrideArraysCore", "ThreadingUtilities"] -git-tree-sha1 = "16bbc30b5ebea91e9ce1671adc03de2832cff552" -uuid = "f517fe37-dbe3-4b94-8317-1923a5111588" -version = "0.7.19" - -[[deps.PolyesterWeave]] -deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"] -git-tree-sha1 = "645bed98cd47f72f67316fd42fc47dee771aefcd" -uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad" -version = "0.2.2" - [[deps.PooledArrays]] deps = ["DataAPI", "Future"] git-tree-sha1 = "36d8b4b899628fb92c2749eb488d884a926614d3" uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720" version = "1.4.3" +[[deps.PositiveIntegrators]] +deps = ["FastBroadcast", "LinearAlgebra", "LinearSolve", "MuladdMacro", "OrdinaryDiffEqCore", "RecipesBase", "Reexport", "SciMLBase", "SparseArrays", "StaticArrays", "Statistics", "SymbolicIndexingInterface"] +git-tree-sha1 = "a10991ac78019690081fcc7d4dcb3ca7cdbeff47" +uuid = "d1b20bf0-b083-4985-a874-dc5121669aa5" +version = "0.2.18" +weakdeps = ["JuMP"] + + [deps.PositiveIntegrators.extensions] + JuMPExt = "JuMP" + [[deps.PreallocationTools]] deps = ["Adapt", "ArrayInterface", "PrecompileTools"] git-tree-sha1 = "e16b73bf892c55d16d53c9c0dbd0fb31cb7e25da" @@ -1860,9 +1854,9 @@ version = "1.2.0" [[deps.PrecompileTools]] deps = ["Preferences"] -git-tree-sha1 = "07a921781cab75691315adc645096ed5e370cb77" +git-tree-sha1 = "edbeefc7a4889f528644251bdb5fc9ab5348bc2c" uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -version = "1.3.3" +version = "1.3.4" [[deps.Preferences]] deps = ["TOML"] @@ -1887,11 +1881,6 @@ deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" version = "1.11.0" -[[deps.Profile]] -deps = ["StyledStrings"] -uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" -version = "1.11.0" - [[deps.ProgressLogging]] deps = ["Logging", "SHA", "UUIDs"] git-tree-sha1 = "f0803bc1171e455a04124affa9c21bba5ac4db32" @@ -1919,9 +1908,9 @@ version = "0.9.31" [[deps.Qt6Base_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Vulkan_Loader_jll", "Xorg_libSM_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_cursor_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "libinput_jll", "xkbcommon_jll"] -git-tree-sha1 = "d7a4bff94f42208ce3cf6bc8e4e7d1d663e7ee8b" +git-tree-sha1 = "144895f6166994730ee7ff8113b981fc360638f1" uuid = "c0090381-4147-56d7-9ebc-da0b1113ec56" -version = "6.10.2+1" +version = "6.10.2+2" [[deps.Qt6Declarative_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Qt6Base_jll", "Qt6ShaderTools_jll", "Qt6Svg_jll"] @@ -1971,16 +1960,19 @@ version = "0.6.12" [[deps.RecursiveArrayTools]] deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "LinearAlgebra", "PrecompileTools", "RecipesBase", "StaticArraysCore", "SymbolicIndexingInterface"] -git-tree-sha1 = "76a9c4d3f7b256dd34074e714ffe0ced8d3d21de" +git-tree-sha1 = "57b6fb3932fc8d1fc911f840d2c9de5fe3ba5008" uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" -version = "3.50.0" +version = "4.3.0" [deps.RecursiveArrayTools.extensions] + RecursiveArrayToolsCUDAExt = "CUDA" RecursiveArrayToolsFastBroadcastExt = "FastBroadcast" + RecursiveArrayToolsFastBroadcastPolyesterExt = ["FastBroadcast", "Polyester"] RecursiveArrayToolsForwardDiffExt = "ForwardDiff" RecursiveArrayToolsKernelAbstractionsExt = "KernelAbstractions" RecursiveArrayToolsMeasurementsExt = "Measurements" RecursiveArrayToolsMonteCarloMeasurementsExt = "MonteCarloMeasurements" + RecursiveArrayToolsMooncakeExt = "Mooncake" RecursiveArrayToolsReverseDiffExt = ["ReverseDiff", "Zygote"] RecursiveArrayToolsSparseArraysExt = ["SparseArrays"] RecursiveArrayToolsStatisticsExt = "Statistics" @@ -1990,11 +1982,14 @@ version = "3.50.0" RecursiveArrayToolsZygoteExt = "Zygote" [deps.RecursiveArrayTools.weakdeps] + CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca" + Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6" + Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" @@ -2028,35 +2023,30 @@ version = "1.3.1" [[deps.Revise]] deps = ["CodeTracking", "FileWatching", "InteractiveUtils", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Preferences", "REPL", "UUIDs"] -git-tree-sha1 = "a44fc6ab46efc588f122ee71d6e2ca83b8645ff7" +git-tree-sha1 = "d9383b639663d8220ac9c523927e38bc21cad16a" uuid = "295af30f-e4ad-537b-8983-00126c2a3abe" -version = "3.14.1" +version = "3.14.3" weakdeps = ["Distributed"] [deps.Revise.extensions] DistributedExt = "Distributed" [[deps.Ribasim]] -deps = ["ADTypes", "Accessors", "BasicModelInterface", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "DelimitedFiles", "DiffEqCallbacks", "DifferentiationInterface", "EnumX", "FiniteDiff", "ForwardDiff", "Graphs", "HiGHS", "IterTools", "JuMP", "LinearAlgebra", "LinearSolve", "Logging", "LoggingExtras", "MathOptAnalyzer", "MetaGraphsNext", "Moshi", "NCDatasets", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqTsit5", "PrecompileTools", "Printf", "SQLite", "SciMLBase", "SciMLOperators", "SparseArrays", "SparseConnectivityTracer", "SparseMatrixColorings", "StructArrays", "Tables", "TerminalLoggers", "TranscodingStreams"] +deps = ["ADTypes", "Accessors", "BasicModelInterface", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "DelimitedFiles", "DiffEqBase", "DiffEqCallbacks", "DifferentiationInterface", "EnumX", "FiniteDiff", "ForwardDiff", "Graphs", "HiGHS", "IterTools", "JuMP", "LinearAlgebra", "LinearSolve", "Logging", "LoggingExtras", "MathOptAnalyzer", "MetaGraphsNext", "Moshi", "NCDatasets", "NaNMath", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqTsit5", "PositiveIntegrators", "PrecompileTools", "Printf", "SQLite", "SciMLBase", "SciMLOperators", "SparseArrays", "SparseConnectivityTracer", "SparseMatrixColorings", "StructArrays", "Tables", "TerminalLoggers", "TranscodingStreams"] path = "core" uuid = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635" version = "2026.1.1" [[deps.RuntimeGeneratedFunctions]] deps = ["ExprTools", "SHA", "Serialization"] -git-tree-sha1 = "7257165d5477fd1025f7cb656019dcb6b0512c38" +git-tree-sha1 = "cfcdc949c4660544ab0fdeed169561cb22f835f4" uuid = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" -version = "0.5.17" +version = "0.5.18" [[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" version = "0.7.0" -[[deps.SIMDTypes]] -git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c" -uuid = "94e857df-77ce-4151-89e5-788b33177be4" -version = "0.1.0" - [[deps.SPDX]] deps = ["DataStructures", "Dates", "JSON", "Logging", "SHA", "TimeZones", "UUIDs"] git-tree-sha1 = "b8a683ba47344db3bf9d91c848da375c1f280c3b" @@ -2077,11 +2067,16 @@ git-tree-sha1 = "0b5f220f90642566b65ba86549d1ee4118ab2579" uuid = "76ed43ae-9a5d-5a62-8c75-30186b810ce8" version = "3.51.2+0" +[[deps.SafeTestsets]] +git-tree-sha1 = "81ec49d645af090901120a1542e67ecbbe044db3" +uuid = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +version = "0.1.0" + [[deps.SciMLBase]] -deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "Moshi", "PreallocationTools", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLLogging", "SciMLOperators", "SciMLPublic", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface"] -git-tree-sha1 = "bdbc7c751d36fdb104d437361950d884d489af85" +deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PreallocationTools", "PrecompileTools", "Preferences", "Printf", "Random", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLLogging", "SciMLOperators", "SciMLPublic", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface"] +git-tree-sha1 = "91c452da69d9e349159b71cf5995973b9d6a0367" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -version = "2.153.1" +version = "3.13.0" [deps.SciMLBase.extensions] SciMLBaseChainRulesCoreExt = "ChainRulesCore" @@ -2089,7 +2084,6 @@ version = "2.153.1" SciMLBaseDistributionsExt = "Distributions" SciMLBaseEnzymeExt = "Enzyme" SciMLBaseForwardDiffExt = "ForwardDiff" - SciMLBaseMLStyleExt = "MLStyle" SciMLBaseMakieExt = "Makie" SciMLBaseMeasurementsExt = "Measurements" SciMLBaseMonteCarloMeasurementsExt = "MonteCarloMeasurements" @@ -2109,7 +2103,6 @@ version = "2.153.1" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" - MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca" @@ -2124,15 +2117,15 @@ version = "2.153.1" [[deps.SciMLJacobianOperators]] deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DifferentiationInterface", "FastClosures", "LinearAlgebra", "SciMLBase", "SciMLOperators"] -git-tree-sha1 = "e96d5e96debf7f80a50d0b976a13dea556ccfd3a" +git-tree-sha1 = "7156a5b51cba1bea33a82a036939ead4131f92bc" uuid = "19f34311-ddf3-4b8b-af20-060888a46c0e" -version = "0.1.12" +version = "0.1.13" [[deps.SciMLLogging]] deps = ["Logging", "LoggingExtras", "Preferences"] -git-tree-sha1 = "0161be062570af4042cf6f69e3d5d0b0555b6927" +git-tree-sha1 = "3f98a53703f925cbd5aac5da4924f82ca34d09ab" uuid = "a6db7da4-7206-11f0-1eab-35f2a5dbe1d1" -version = "1.9.1" +version = "2.0.0" [deps.SciMLLogging.extensions] SciMLLoggingTracyExt = "Tracy" @@ -2141,16 +2134,21 @@ version = "1.9.1" Tracy = "e689c965-62c8-4b79-b2c5-8359227902fd" [[deps.SciMLOperators]] -deps = ["Accessors", "ArrayInterface", "DocStringExtensions", "LinearAlgebra"] -git-tree-sha1 = "794c760e6aafe9f40dcd7dd30526ea33f0adc8b7" +deps = ["Accessors", "Adapt", "ArrayInterface", "DocStringExtensions", "LinearAlgebra", "SafeTestsets"] +git-tree-sha1 = "50e6ec6879eab12b039924d5a10b91c95bf9bf7f" uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" -version = "1.15.1" -weakdeps = ["SparseArrays", "StaticArraysCore"] +version = "1.21.0" [deps.SciMLOperators.extensions] + SciMLOperatorsLoopVectorizationExt = "LoopVectorization" SciMLOperatorsSparseArraysExt = "SparseArrays" SciMLOperatorsStaticArraysCoreExt = "StaticArraysCore" + [deps.SciMLOperators.weakdeps] + LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" + [[deps.SciMLPublic]] git-tree-sha1 = "0ba076dbdce87ba230fff48ca9bca62e1f345c9b" uuid = "431bcebd-1456-4ced-9d72-93c2757fff0b" @@ -2197,9 +2195,9 @@ version = "1.2.0" [[deps.SimpleNonlinearSolve]] deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "MaybeInplace", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase", "Setfield", "StaticArraysCore"] -git-tree-sha1 = "744c3f0fb186ad28376199c1e72ca39d0c614b5d" +git-tree-sha1 = "d688de789b7e643326caf9a1051376dadbcd8873" uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7" -version = "2.11.0" +version = "2.11.1" [deps.SimpleNonlinearSolve.extensions] SimpleNonlinearSolveChainRulesCoreExt = "ChainRulesCore" @@ -2254,22 +2252,25 @@ version = "1.2.1" [[deps.SparseMatrixColorings]] deps = ["ADTypes", "DocStringExtensions", "LinearAlgebra", "PrecompileTools", "Random", "SparseArrays"] -git-tree-sha1 = "1c1be8c6fdfaf9b6c9e156c509e672953b8e6af7" +git-tree-sha1 = "f63d76c7b7c329cf11badd564fd8ba877b09c3fe" uuid = "0a514795-09f3-496d-8182-132a7b665d35" -version = "0.4.26" +version = "0.4.27" [deps.SparseMatrixColorings.extensions] - SparseMatrixColoringsCUDAExt = "CUDA" + SparseMatrixColoringsCUDAExt = ["CUDA", "cuSPARSE"] SparseMatrixColoringsCliqueTreesExt = "CliqueTrees" SparseMatrixColoringsColorsExt = "Colors" + SparseMatrixColoringsGPUArraysExt = "GPUArrays" SparseMatrixColoringsJuMPExt = ["JuMP", "MathOptInterface"] [deps.SparseMatrixColorings.weakdeps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" + GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" + cuSPARSE = "b26da814-b3bc-49ef-b0ee-c816305aa060" [[deps.SpecialFunctions]] deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] @@ -2289,23 +2290,6 @@ git-tree-sha1 = "4f96c596b8c8258cc7d3b19797854d368f243ddc" uuid = "860ef19b-820b-49d6-a774-d7a799459cd3" version = "1.0.4" -[[deps.Static]] -deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools", "SciMLPublic"] -git-tree-sha1 = "49440414711eddc7227724ae6e570c7d5559a086" -uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "1.3.1" - -[[deps.StaticArrayInterface]] -deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "SciMLPublic", "Static"] -git-tree-sha1 = "aa1ea41b3d45ac449d10477f65e2b40e3197a0d2" -uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718" -version = "1.9.0" -weakdeps = ["OffsetArrays", "StaticArrays"] - - [deps.StaticArrayInterface.extensions] - StaticArrayInterfaceOffsetArraysExt = "OffsetArrays" - StaticArrayInterfaceStaticArraysExt = "StaticArrays" - [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] git-tree-sha1 = "246a8bb2e6667f832eea063c3a56aef96429a3db" @@ -2347,12 +2331,6 @@ git-tree-sha1 = "aceda6f4e598d331548e04cc6b2124a6148138e3" uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" version = "0.34.10" -[[deps.StrideArraysCore]] -deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface", "ThreadingUtilities"] -git-tree-sha1 = "83151ba8065a73f53ca2ae98bc7274d817aa30f2" -uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da" -version = "0.5.8" - [[deps.StringManipulation]] deps = ["PrecompileTools"] git-tree-sha1 = "d05693d339e37d6ab134c5ab53c29fce5ee5d7d5" @@ -2402,9 +2380,9 @@ version = "7.8.3+2" [[deps.SymbolicIndexingInterface]] deps = ["Accessors", "ArrayInterface", "RuntimeGeneratedFunctions", "StaticArraysCore"] -git-tree-sha1 = "94c58884e013efff548002e8dc2fdd1cb74dfce5" +git-tree-sha1 = "173ecfe5f7c5a36043b5f2b8cecaa30a4fc958ef" uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5" -version = "0.3.46" +version = "0.3.47" weakdeps = ["PrettyTables"] [deps.SymbolicIndexingInterface.extensions] @@ -2466,12 +2444,6 @@ git-tree-sha1 = "42fd9023fef18b9b78c8343a4e2f3813ffbcefcb" uuid = "1c621080-faea-4a02-84b6-bbd5e436b8fe" version = "1.0.0" -[[deps.ThreadingUtilities]] -deps = ["ManualMemory"] -git-tree-sha1 = "d969183d3d244b6c33796b5ed01ab97328f2db85" -uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5" -version = "0.5.5" - [[deps.TimeZones]] deps = ["Artifacts", "Dates", "Downloads", "InlineStrings", "Mocking", "Printf", "Scratch", "TZJData", "Unicode", "p7zip_jll"] git-tree-sha1 = "d422301b2a1e294e3e4214061e44f338cafe18a2" @@ -2549,9 +2521,9 @@ version = "1.24.0+0" [[deps.WeakRefStrings]] deps = ["DataAPI", "InlineStrings", "Parsers"] -git-tree-sha1 = "b1be2855ed9ed8eac54e5caff2afcdb442d52c23" +git-tree-sha1 = "0716e01c3b40413de5dedbc9c5c69f27cddfddfc" uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5" -version = "1.4.2" +version = "1.4.3" [[deps.WorkerUtilities]] git-tree-sha1 = "cd1659ba0d57b71a464a29e64dbc67cfe83d54e7" @@ -2566,9 +2538,9 @@ version = "2.13.9+0" [[deps.XZ_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "9cce64c0fdd1960b597ba7ecda2950b5ed957438" +git-tree-sha1 = "b29c22e245d092b8b4e8d3c09ad7baa586d9f573" uuid = "ffd25f8a-64ca-5728-b0f7-c24cf3aae800" -version = "5.8.2+0" +version = "5.8.3+0" [[deps.Xorg_libICE_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -2644,9 +2616,9 @@ version = "0.9.12+0" [[deps.Xorg_libpciaccess_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"] -git-tree-sha1 = "4909eb8f1cbf6bd4b1c30dd18b2ead9019ef2fad" +git-tree-sha1 = "58972370b81423fc546c56a60ed1a009450177c3" uuid = "a65dc6b1-eb27-53a1-bb3e-dea574b5389e" -version = "0.18.1+0" +version = "0.19.0+0" [[deps.Xorg_libxcb_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libXau_jll", "Xorg_libXdmcp_jll"] @@ -2704,9 +2676,9 @@ version = "1.4.7+0" [[deps.Xorg_xkeyboard_config_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_xkbcomp_jll"] -git-tree-sha1 = "00af7ebdc563c9217ecc67776d1bbf037dbcebf4" +git-tree-sha1 = "429722587208f02b1cecbddcd20133df2f1ed796" uuid = "33bec58e-1273-512f-9401-5d533626f822" -version = "2.44.0+0" +version = "2.47.0+0" [[deps.Xorg_xtrans_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -2725,6 +2697,60 @@ git-tree-sha1 = "446b23e73536f84e8037f5dce465e92275f6a308" uuid = "3161d3a3-bdf6-5164-811a-617609db77b4" version = "1.5.7+1" +[[deps.aws_c_auth_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "aws_c_cal_jll", "aws_c_http_jll", "aws_c_sdkutils_jll"] +git-tree-sha1 = "8cab83c96af80a1be968251ce1a0548a7545484d" +uuid = "2b3700d1-4306-52e2-a478-c162f0c514be" +version = "0.9.6+0" + +[[deps.aws_c_cal_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "aws_c_common_jll"] +git-tree-sha1 = "22c0f42f4a1f0dc5dcfa8fd267c4ac407c455e7a" +uuid = "70f11efc-bab2-57f1-b0f3-22aad4e67c4b" +version = "0.9.13+0" + +[[deps.aws_c_common_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "a759cb9bf456ad792cc7898a81ae333cce9ef02a" +uuid = "73048d1d-b8c4-5092-a58d-866c5e8d1e50" +version = "0.12.6+0" + +[[deps.aws_c_compression_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "aws_c_common_jll"] +git-tree-sha1 = "7910c72f45f44afd297c39fe43b99c56d5ed22ec" +uuid = "73a04cd5-f3d7-5bac-9290-e8adb709f224" +version = "0.3.2+0" + +[[deps.aws_c_http_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "aws_c_compression_jll", "aws_c_io_jll"] +git-tree-sha1 = "e358d5a001ef7afbd4f8c5225322512819cda2f2" +uuid = "3254fc65-9028-534d-aa9d-d76d128babc6" +version = "0.10.13+0" + +[[deps.aws_c_io_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "aws_c_cal_jll", "aws_c_common_jll", "s2n_tls_jll"] +git-tree-sha1 = "7e481d474b2087ee8bbf55b81bf9119f21e396d9" +uuid = "13c41daa-f319-5298-b5eb-5754e0170d52" +version = "0.26.3+0" + +[[deps.aws_c_s3_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "aws_c_auth_jll", "aws_c_common_jll", "aws_c_http_jll", "aws_checksums_jll", "s2n_tls_jll"] +git-tree-sha1 = "3e9917ab25114feba657e71be41cad068b9f6595" +uuid = "bd1f34fb-993f-5903-a121-aaf302eed6d4" +version = "0.11.5+0" + +[[deps.aws_c_sdkutils_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "aws_c_common_jll"] +git-tree-sha1 = "c43dfba2c1ab9ea9f02f2c80e86fa16f6460244e" +uuid = "1282aa60-004d-510b-9f52-12498d409daa" +version = "0.2.4+1" + +[[deps.aws_checksums_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "aws_c_common_jll"] +git-tree-sha1 = "2570c8e23f4771a087b12a47edcaaa670ac05a01" +uuid = "b2a88e68-78e7-5e94-8c20-c02986ec140e" +version = "0.2.10+0" + [[deps.dlfcn_win32_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] git-tree-sha1 = "e141d67ffe550eadfb5af1bdbdaf138031e4805f" @@ -2745,15 +2771,15 @@ version = "0.61.1+0" [[deps.libaec_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "13b760f97c6e753b47df30cb438d4dc3b50df282" +git-tree-sha1 = "1411bc34c180946d3cef591de1384012afa6edee" uuid = "477f73a3-ac25-53e9-8cc3-50b2fa2566f0" -version = "1.1.5+0" +version = "1.1.6+0" [[deps.libaom_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "371cc681c00a3ccc3fbc5c0fb91f58ba9bec1ecf" +git-tree-sha1 = "850b06095ee71f0135d644ffd8a52850699581ed" uuid = "a4ae2306-e953-59d6-aa16-d00cac43593b" -version = "3.13.1+0" +version = "3.13.3+0" [[deps.libass_jll]] deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Zlib_jll"] @@ -2798,9 +2824,9 @@ version = "1.28.1+0" [[deps.libpng_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"] -git-tree-sha1 = "e2a7072fc0cdd7949528c1455a3e5da4122e1153" +git-tree-sha1 = "e51150d5ab85cee6fc36726850f0e627ad2e4aba" uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f" -version = "1.6.56+0" +version = "1.6.58+0" [[deps.libva_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll", "Xorg_libXext_jll", "Xorg_libXfixes_jll", "libdrm_jll"] @@ -2832,6 +2858,12 @@ git-tree-sha1 = "717df6f6892af4ee13279a73aa58474e58a88667" uuid = "f8abcde7-e9b7-5caa-b8af-a437887ae8e4" version = "2.3.1+0" +[[deps.mpif_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIABI_jll", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "TOML"] +git-tree-sha1 = "a8083ee0737c243c8f40a4ba86a0956997facb73" +uuid = "9aeb927a-4695-514f-a259-621a69f20ec0" +version = "0.1.7+0" + [[deps.mtdev_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] git-tree-sha1 = "b4d631fd51f2e9cdd93724ae25b2efc198b059b1" @@ -2845,9 +2877,9 @@ version = "1.64.0+1" [[deps.oneTBB_jll]] deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"] -git-tree-sha1 = "1350188a69a6e46f799d3945beef36435ed7262f" +git-tree-sha1 = "da8c1f6eee04831f14edcfa5dae611d309807e57" uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e" -version = "2022.0.0+1" +version = "2022.3.0+0" [[deps.p7zip_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] @@ -2860,6 +2892,12 @@ git-tree-sha1 = "f349584316617063160a947a82638f7611a8ef0f" uuid = "4d7b5844-a134-5dcd-ac86-c8f19cd51bed" version = "0.41.3+0" +[[deps.s2n_tls_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "64ae051c6f03044eb7d98027d1b552b4e21e650c" +uuid = "cddc5d3d-934d-5d3a-9747-62fc12ea3f48" +version = "1.7.3+0" + [[deps.x264_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] git-tree-sha1 = "14cc7083fc6dff3cc44f2bc435ee96d06ed79aa7" diff --git a/Project.toml b/Project.toml index 7e84c56af..1a764733b 100644 --- a/Project.toml +++ b/Project.toml @@ -59,6 +59,7 @@ OteraEngine = "b2d7f28f-acd6-4007-8b26-bc27716e5513" PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" PkgToSoftwareBOM = "6254a0f9-6143-4104-aa2e-fd339a2830a6" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +PositiveIntegrators = "d1b20bf0-b083-4985-a874-dc5121669aa5" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" diff --git a/core/Project.toml b/core/Project.toml index 1dbed9c4f..18522198b 100644 --- a/core/Project.toml +++ b/core/Project.toml @@ -16,6 +16,7 @@ DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" @@ -33,6 +34,7 @@ MathOptAnalyzer = "d1179b25-476b-425c-b826-c7787f0fff83" MetaGraphsNext = "fa8bd995-216d-47f1-8a91-f3b68fbeb377" Moshi = "2e0e35c7-a2e4-4343-998d-7ef72827ed2d" NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" +NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" OrdinaryDiffEqDifferentiation = "4302a76b-040a-498a-8c04-15b101fed76b" @@ -41,6 +43,7 @@ OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" +PositiveIntegrators = "d1b20bf0-b083-4985-a874-dc5121669aa5" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9" @@ -64,11 +67,12 @@ DataInterpolations = "8.3" DataStructures = "0.18, 0.19" Dates = "1" DelimitedFiles = "1.9.1" -DiffEqCallbacks = "3.6, 4" -DifferentiationInterface = "0.6.48, 0.7" +DiffEqBase = "7" +DiffEqCallbacks = "4" +DifferentiationInterface = "0.7" EnumX = "1.0" FiniteDiff = "2.21" -ForwardDiff = "0.10.38, 1" +ForwardDiff = "1" Graphs = "1.9" HiGHS = "1.7" IterTools = "1.4" @@ -81,18 +85,19 @@ MathOptAnalyzer = "0.1.0" MetaGraphsNext = "0.6, 0.7, 0.8" Moshi = "0.3.7" NCDatasets = "0.14.8" -OrdinaryDiffEqBDF = "1.14, 2" -OrdinaryDiffEqCore = "3.2, 4" -OrdinaryDiffEqDifferentiation = "2, 3" -OrdinaryDiffEqLowOrderRK = "1.10, 2" -OrdinaryDiffEqNonlinearSolve = "1.19, 2" -OrdinaryDiffEqRosenbrock = "1.22, 2" -OrdinaryDiffEqSDIRK = "1.11, 2" -OrdinaryDiffEqTsit5 = "1.9, 2" +NaNMath = "1" +OrdinaryDiffEqBDF = "2.1.1" +OrdinaryDiffEqCore = "4" +OrdinaryDiffEqDifferentiation = "3" +OrdinaryDiffEqLowOrderRK = "2" +OrdinaryDiffEqNonlinearSolve = "2" +OrdinaryDiffEqRosenbrock = "2" +OrdinaryDiffEqSDIRK = "2" +OrdinaryDiffEqTsit5 = "2" PrecompileTools = "1.2.1" Printf = "1" SQLite = "1.5.1" -SciMLBase = "2.36, 3" +SciMLBase = "3" SciMLOperators = "1.15.1" SparseArrays = "1" SparseConnectivityTracer = "1" diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index 963fe9884..63e774cf4 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -32,13 +32,14 @@ using DifferentiationInterface: using ForwardDiff: derivative as forward_diff # Algorithms for solving ODEs. -using OrdinaryDiffEqCore: OrdinaryDiffEqCore, get_du -using OrdinaryDiffEqDifferentiation: - OrdinaryDiffEqDifferentiation, dolinsolve, jacobian2W! -using SciMLOperators: WOperator, MatrixOperator +using OrdinaryDiffEqCore: OrdinaryDiffEqCore, get_du, AbstractNLSolver +using DiffEqBase: DiffEqBase, calculate_residuals! +using OrdinaryDiffEqNonlinearSolve: OrdinaryDiffEqNonlinearSolve, relax!, _compute_rhs! import ADTypes using ADTypes: AutoForwardDiff import ForwardDiff +import NaNMath +import OrdinaryDiffEqBDF # Interface for defining and solving the ODE problem of the physical layer. using SciMLBase: @@ -51,12 +52,7 @@ using SciMLBase: ODEProblem, get_proposed_dt, DEIntegrator, - FullSpecialize, - NoSpecialize, - SciMLOperators, - AbstractSciMLOperator, - LinearProblem, - LinearSolution + FullSpecialize # Automatically detecting the sparsity pattern of the Jacobian of water_balance! # through operator overloading @@ -64,10 +60,13 @@ using SparseConnectivityTracer: GradientTracer, TracerSparsityDetector using SparseMatrixColorings: GreedyColoringAlgorithm, sparsity_pattern # For efficient sparse computations -using SparseArrays: SparseMatrixCSC, sparse, nzrange +using SparseArrays: SparseMatrixCSC, sparse, nzrange, nonzeros + +# Positivity-preserving time integration +using PositiveIntegrators: PDSProblem, MPRK22, MPRK43I, MPRK43II # Linear algebra -using LinearAlgebra: LinearAlgebra, mul!, UniformScaling +using LinearAlgebra: LinearAlgebra, mul!, UniformScaling, Symmetric, cholesky, Factorization # Interpolation functionality, used for e.g. # basin profiles and TabulatedRatingCurve. See also the node @@ -182,7 +181,7 @@ include("allocation_init.jl") include("allocation_optim.jl") include("util.jl") include("graph.jl") -include("differentiation.jl") +include("pds.jl") include("model.jl") include("read.jl") include("write.jl") diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 836c8f4f8..7a15af05e 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -14,7 +14,6 @@ function set_simulation_data!( user_demand, tabulated_rating_curve, ) = p.p_independent - du = get_du(integrator) errors = false @@ -24,7 +23,7 @@ function set_simulation_data!( set_simulation_data!(allocation_model, linear_resistance, p, t) set_simulation_data!(allocation_model, manning_resistance, p, t) set_simulation_data!(allocation_model, tabulated_rating_curve, p, t) - set_simulation_data!(allocation_model, pump, outlet, du) + set_simulation_data!(allocation_model, pump, outlet, p.state_and_time_dependent_cache) set_simulation_data!(allocation_model, user_demand, t) if errors @@ -350,7 +349,7 @@ function set_simulation_data!( allocation_model::AllocationModel, pump::Pump, outlet::Outlet, - du::CVector, + cache::StateAndTimeDependentCache, )::Nothing (; problem, scaling) = allocation_model pump_constraints = problem[:pump] @@ -360,7 +359,7 @@ function set_simulation_data!( for node_id in only(pump_constraints.axes) constraint = pump_constraints[node_id] upstream_node_id = pump.inflow_link[node_id.idx].link[1] - q = du.pump[node_id.idx] + q = cache.current_flow_rate_pump[node_id.idx] if upstream_node_id.type == NodeType.Basin low_storage_factor = get_low_storage_factor(problem, upstream_node_id) JuMP.set_normalized_coefficient( @@ -377,7 +376,7 @@ function set_simulation_data!( for node_id in only(outlet_constraints.axes) constraint = outlet_constraints[node_id] upstream_node_id = outlet.inflow_link[node_id.idx].link[1] - q = du.outlet[node_id.idx] + q = cache.current_flow_rate_outlet[node_id.idx] if upstream_node_id.type == NodeType.Basin low_storage_factor = get_low_storage_factor(problem, upstream_node_id) JuMP.set_normalized_coefficient( @@ -769,19 +768,21 @@ function warm_start!(allocation_model::AllocationModel, integrator::DEIntegrator (; link_to_state_idx) = p.p_independent # Extrapolate the current instantaneous flow rates from the physical layer + # Flow rates are now stored in cache vectors, look up by link + flow_link_lookup = p.p_independent.graph[].flow_link_lookup for link in only(flow.axes) - state_index = get_state_index(getaxes(du), link_to_state_idx, link) - if !isnothing(state_index) - JuMP.set_start_value(flow[link], du[state_index] / scaling.flow) + link_idx = get_link_index(link, flow_link_lookup) + if !isnothing(link_idx) + JuMP.set_start_value(flow[link], p.p_independent.current_flow_rate[link_idx] / scaling.flow) end end # Extrapolate the current instantaneous storage rates from the physical layer + # du.basin[idx] now directly contains dS/dt for node_id in basin_ids_subnetwork JuMP.set_start_value( storage_change[node_id], - formulate_dstorage_wrt_time(du, p.p_independent, t, node_id) * Δt_allocation / - scaling.storage, + du.basin[node_id.idx] * Δt_allocation / scaling.storage, ) end diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 603e3f8ec..370f0ff8d 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -72,7 +72,7 @@ function BMI.get_value_ptr(model::Model, name::String)::Vector{Float64} elseif name == "basin.surface_runoff" basin.vertical_flux.surface_runoff::Vector{Float64} elseif name == "basin.cumulative_infiltration" - unsafe_array(u.infiltration)::Vector{Float64} + p_independent.cumulative_infiltration_total::Vector{Float64} elseif name == "basin.cumulative_drainage" basin.cumulative_drainage::Vector{Float64} elseif name == "basin.cumulative_surface_runoff" @@ -82,7 +82,7 @@ function BMI.get_value_ptr(model::Model, name::String)::Vector{Float64} elseif name == "user_demand.demand" vec(user_demand.demand)::Vector{Float64} elseif name == "user_demand.cumulative_inflow" - unsafe_array(u.user_demand_inflow)::Vector{Float64} + state_and_time_dependent_cache.current_flow_rate_user_demand else error("Unknown variable $name") end diff --git a/core/src/callback.jl b/core/src/callback.jl index 2580ed36b..b25c17c21 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -66,11 +66,6 @@ function create_callbacks( discrete_control_cb = FunctionCallingCallback(apply_discrete_control!) push!(callbacks, discrete_control_cb) - toltimes = get_log_tstops(config.starttime, config.endtime) - decrease_tol_cb = - FunctionCallingCallback(decrease_tolerance!; funcat = toltimes, func_start = false) - push!(callbacks, decrease_tol_cb) - saved = SavedResults( saved_flow, saved_basin_states, @@ -82,69 +77,330 @@ function create_callbacks( return callback, saved end +function sync_flow_rates!(node_sym::Symbol, current_flow_cache::Vector{T}, p_independent::ParametersIndependent) where {T} + flow_link_lookup = p_independent.graph[].flow_link_lookup + node = getproperty(p_independent, node_sym) + for (node_idx, id) in enumerate(node.node_id) + q = current_flow_cache[id.idx] + inflow_link = node.inflow_link[node_idx] + outflow_link = node.outflow_link[node_idx] + idx = get_link_index(inflow_link.link, flow_link_lookup) + if idx !== nothing + p_independent.current_flow_rate[idx] = q + end + idx = get_link_index(outflow_link.link, flow_link_lookup) + if idx !== nothing + p_independent.current_flow_rate[idx] = q + end + end + return +end + """ -Decrease the relative tolerance of the integrator over time, -to compensate for the ever increasing cumulative flows. +Synchronize UserDemand flow rates into the per-link flow rate vector. +UserDemand has multiple inflow links per node, so it needs special handling. """ -function decrease_tolerance!(u, t, integrator)::Nothing - (; p, t, opts) = integrator - - for (i, state) in enumerate(u) - p.p_independent.relmask[i] || continue - - # Use the internal norm to get the magnitude of the (cumulative) states, - # as used in calculate_residuals, and compare to an estimated average magnitude - cum_magnitude = opts.internalnorm(state, t) - iszero(cum_magnitude) && continue - avg_magnitude = max(opts.internalnorm(1.0e4, t), cum_magnitude / t) # allow for 1e4 m3/s - - # Decrease the relative tolerance based on their difference - diff_norm = max(0, log10(cum_magnitude / avg_magnitude)) - # Limit new tolerance to floating point precision (~-14) - newtol = max(10.0^(log10(integrator.p.p_independent.reltol) - diff_norm), 1.0e-14) - - if opts.reltol[i] > newtol - @debug "Relative tolerance changed at t = $t, state = $i to $(newtol)" - opts.reltol[i] = newtol +function sync_user_demand_flow_rates!(current_flow_cache::Vector{T}, p_independent::ParametersIndependent) where {T} + flow_link_lookup = p_independent.graph[].flow_link_lookup + user_demand = p_independent.user_demand + for (node_idx, id) in enumerate(user_demand.node_id) + q = current_flow_cache[id.idx] + # Sync outflow link + outflow_link = user_demand.outflow_link[node_idx] + idx = get_link_index(outflow_link.link, flow_link_lookup) + if idx !== nothing + p_independent.current_flow_rate[idx] = q + end + # Sync inflow links (distribute evenly for trapezoidal accumulation) + inflow_links = user_demand.inflow_links[node_idx] + n_links = length(inflow_links) + q_per_link = n_links == 0 ? zero(T) : q / n_links + for link_meta in inflow_links + idx = get_link_index(link_meta.link, flow_link_lookup) + if idx !== nothing + p_independent.current_flow_rate[idx] = q_per_link + end end end return end """ -Update with the latest timestep: -- Cumulative flows/forcings which are integrated exactly -- Cumulative flows/forcings which are input for the allocation algorithm -- Cumulative flows/forcings which are supplied demands in the allocation context - -During these cumulative flow updates, we can also update the mass balance of the system, -as each flow carries mass, based on the concentrations of the flow source. -Specifically, we first use all the inflows to update the mass of the Basins, recalculate -the Basin concentration(s) and then remove the mass that is being lost to the outflows. +Synchronize per-node-type flow rate caches into the per-link flow rate vector +used for trapezoidal flow accumulation. +""" +function sync_flow_rates!(p::Parameters)::Nothing + (; p_independent, state_and_time_dependent_cache, time_dependent_cache) = p + (; basin) = p_independent + cache = state_and_time_dependent_cache + + sync_flow_rates!(:pump, cache.current_flow_rate_pump, p_independent) + sync_flow_rates!(:outlet, cache.current_flow_rate_outlet, p_independent) + sync_flow_rates!(:linear_resistance, cache.current_flow_rate_linear_resistance, p_independent) + sync_flow_rates!(:tabulated_rating_curve, cache.current_flow_rate_tabulated_rating_curve, p_independent) + sync_flow_rates!(:manning_resistance, cache.current_flow_rate_manning_resistance, p_independent) + sync_user_demand_flow_rates!(cache.current_flow_rate_user_demand, p_independent) + + return nothing +end + +""" +Solve a constrained least-squares balance correction when the unconstrained +projection would violate nonnegative flow physics on one-way links. + +Minimize 0.5*(||dq||^2 + ||de||^2 + ||di||^2) +subject to: A_flow * dq + de + di == residual, +and q* + dq >= 0 on constrained links, +and corrected evaporation/infiltration remaining nonnegative. + +Returns true when an optimal solution was found and writes the solution into +`correction_flow`, `evap_correction`, and `infiltration_correction`. +""" +function solve_constrained_balance_correction!( + correction_flow::Vector{Float64}, + evap_correction::Vector{Float64}, + infiltration_correction::Vector{Float64}, + residual::Vector{Float64}, + A_flow::SparseMatrixCSC{Float64, Int}, + cumulative_flow::Vector{Float64}, + cumulative_evaporation::Vector{Float64}, + cumulative_infiltration::Vector{Float64}, + nonnegative_link::BitVector, + )::Bool + n_links = length(cumulative_flow) + n_basins = length(residual) + + model = JuMP.Model(HiGHS.Optimizer) + JuMP.set_silent(model) + + JuMP.@variable(model, dq[1:n_links]) + JuMP.@variable(model, de[1:n_basins]) + JuMP.@variable(model, di[1:n_basins]) + + for j in 1:n_links + if nonnegative_link[j] + JuMP.set_lower_bound(dq[j], -cumulative_flow[j]) + end + end + for i in 1:n_basins + # corrected fluxes are e* - de and i* - di + JuMP.set_upper_bound(de[i], cumulative_evaporation[i]) + JuMP.set_upper_bound(di[i], cumulative_infiltration[i]) + end + + JuMP.@constraint(model, A_flow * dq + de + di .== residual) + JuMP.@objective( + model, + Min, + 0.5 * sum(dq[j]^2 for j in 1:n_links) + + 0.5 * sum(de[i]^2 for i in 1:n_basins) + + 0.5 * sum(di[i]^2 for i in 1:n_basins), + ) + + JuMP.optimize!(model) + if JuMP.termination_status(model) != JuMP.OPTIMAL + return false + end + + for j in 1:n_links + correction_flow[j] = JuMP.value(dq[j]) + end + for i in 1:n_basins + evap_correction[i] = JuMP.value(de[i]) + infiltration_correction[i] = JuMP.value(di[i]) + end + return true +end + +""" +Apply mass-balance-consistent correction to trapezoidal flow estimates. +Projects the trapezoidal estimates onto the subspace satisfying exact mass balance +using the minimum-norm adjustment: δq = Aᵀ(AAᵀ+2I)⁻¹(b - Aq̄*) + +The correction distributes the per-basin water balance residual back to the +internal flow links, evaporation, and infiltration in proportion to the network +topology, guaranteeing mass conservation at every step. +""" +function apply_balance_correction!(u, p_independent, time_dependent_cache)::Nothing + (; basin, flow_boundary, balance_correction) = p_independent + (; A_flow, AAt_2I_chol, nonnegative_link, storage_prev, lambda, residual, correction_flow) = + balance_correction + n_basins = length(basin.node_id) + + if n_basins == 0 + return nothing + end + + # b[i] = ΔS[i] - exact_forcings[i] + # where exact_forcings = precipitation + drainage + surface_runoff + for i in 1:n_basins + ΔS = u.basin[i] - storage_prev[i] + precip = time_dependent_cache.basin.current_cumulative_precipitation[i] - + basin.cumulative_precipitation[i] + drainage = time_dependent_cache.basin.current_cumulative_drainage[i] - + basin.cumulative_drainage[i] + surface_runoff = time_dependent_cache.basin.current_cumulative_surface_runoff[i] - + basin.cumulative_surface_runoff[i] + residual[i] = ΔS - precip - drainage - surface_runoff + end + + # Subtract exact flow boundary contributions to receiving basins + for (outflow_link, id) in zip(flow_boundary.outflow_link, flow_boundary.node_id) + dst = outflow_link.link[2] + if dst.type == NodeType.Basin + fb_vol = + time_dependent_cache.flow_boundary.current_cumulative_boundary_flow[id.idx] - + flow_boundary.cumulative_flow[id.idx] + residual[dst.idx] -= fb_vol + end + end + + # Subtract A_ext * q_ext* (the trapezoidal estimate contribution) + # A_ext * q_ext* = A_flow * cumulative_flow - cumulative_evaporation - cumulative_infiltration + mul!(lambda, A_flow, p_independent.cumulative_flow) # reuse lambda as temp + @. residual -= lambda + @. residual += p_independent.cumulative_evaporation + @. residual += p_independent.cumulative_infiltration + # Now residual = b - A_ext * q_ext* = per-basin mass balance error (in volumes) + + # Solve for Lagrange multipliers: λ = (AAᵀ + 2I)⁻¹ r + lambda .= AAt_2I_chol \ residual + + # Compute flow corrections: δq_flow = Aᵀλ + mul!(correction_flow, A_flow', lambda) + + # If unconstrained correction would violate one-way link orientation, + # solve a constrained least-squares problem instead. + # constrained_correction = false + # for j in eachindex(correction_flow) + # if nonnegative_link[j] && (p_independent.cumulative_flow[j] + correction_flow[j] < 0.0) + # constrained_correction = true + # break + # end + # end + + # if constrained_correction + # # Reuse lambda/residual vectors as work arrays for evap/infiltration corrections. + # success = solve_constrained_balance_correction!( + # correction_flow, + # lambda, + # residual, + # residual, + # A_flow, + # p_independent.cumulative_flow, + # p_independent.cumulative_evaporation, + # p_independent.cumulative_infiltration, + # nonnegative_link, + # ) + # if !success + # @warn "Constrained balance correction failed; using unconstrained projection." + # lambda .= AAt_2I_chol \ residual + # mul!(correction_flow, A_flow', lambda) + + # @. p_independent.cumulative_flow += correction_flow + # @. p_independent.cumulative_flow_saveat += correction_flow + # @. p_independent.cumulative_evaporation -= lambda + # @. p_independent.cumulative_evaporation_saveat -= lambda + # @. p_independent.cumulative_infiltration -= lambda + # @. p_independent.cumulative_infiltration_saveat -= lambda + # else + # @. p_independent.cumulative_flow += correction_flow + # @. p_independent.cumulative_flow_saveat += correction_flow + # @. p_independent.cumulative_evaporation -= lambda + # @. p_independent.cumulative_evaporation_saveat -= lambda + # @. p_independent.cumulative_infiltration -= residual + # @. p_independent.cumulative_infiltration_saveat -= residual + # end + # else + # Apply corrections to per-step and saveat accumulators + @. p_independent.cumulative_flow += correction_flow + @. p_independent.cumulative_flow_saveat += correction_flow + @. p_independent.cumulative_evaporation -= lambda + @. p_independent.cumulative_evaporation_saveat -= lambda + @. p_independent.cumulative_infiltration -= lambda + @. p_independent.cumulative_infiltration_saveat -= lambda + # end + + # Accumulate absolute correction for per-link flow convergence output + @. p_independent.flow_convergence_saveat += abs(correction_flow) + + # Update storage_prev for next step + for i in 1:n_basins + storage_prev[i] = u.basin[i] + end + + return nothing +end + +""" +Update cumulative flows. +Also updates cumulative forcings and allocation flows. """ function update_cumulative_flows!(u, t, integrator)::Nothing - (; cache, p) = integrator - (; p_independent, p_mutable, time_dependent_cache) = p - (; basin, flow_boundary, allocation, temp_convergence, convergence, ncalls) = - p_independent + (; p) = integrator + (; p_independent, p_mutable, time_dependent_cache, state_and_time_dependent_cache) = p + (; basin, flow_boundary, allocation) = p_independent + + dt = t - p_mutable.tprev # Update tprev p_mutable.tprev = t - # Update convergence measure - if hasproperty(cache, :nlsolver) - @. temp_convergence = abs(cache.nlsolver.cache.atmp / u) - @inbounds for I in eachindex(temp_convergence) - if !isfinite(temp_convergence[I]) - temp_convergence[I] = zero(eltype(temp_convergence)) + # Sync per-node-type flow rates to per-link vector + sync_flow_rates!(p) + + # Compute per-step flow volumes via trapezoidal integration. + # cumulative_flow holds the per-step volume (used by concentration and allocation), + # cumulative_flow_saveat accumulates across the saveat interval (used for output). + p_independent.cumulative_flow .= 0.0 + if dt > 0 + # Use trapezoidal integration, but fall back to backward Euler for links where + # flow_rate_prev was set to NaN as a sentinel at a control state transition. + # This avoids a spurious spike from averaging the pre-transition and post-transition + for i in eachindex(p_independent.cumulative_flow) + curr = p_independent.current_flow_rate[i] + prev = p_independent.flow_rate_prev[i] + p_independent.cumulative_flow[i] = + isnan(prev) ? dt * curr : 0.5 * dt * (curr + prev) + p_independent.flow_rate_prev[i] = curr + end + @. p_independent.cumulative_flow_saveat += p_independent.cumulative_flow + + @. p_independent.cumulative_evaporation = 0.5 * dt * (state_and_time_dependent_cache.current_evaporation + p_independent.evaporation_prev) + @. p_independent.cumulative_evaporation_saveat += p_independent.cumulative_evaporation + @. p_independent.evaporation_prev = state_and_time_dependent_cache.current_evaporation + + @. p_independent.cumulative_infiltration = 0.5 * dt * (state_and_time_dependent_cache.current_infiltration + p_independent.infiltration_prev) + @. p_independent.cumulative_infiltration_saveat += p_independent.cumulative_infiltration + @. p_independent.infiltration_prev = state_and_time_dependent_cache.current_infiltration + + # Accumulate running total for BMI using current rate * dt (like an ODE state) + @. p_independent.cumulative_infiltration_total += dt * state_and_time_dependent_cache.current_infiltration + + # Project trapezoidal estimates onto mass-balance-consistent subspace + apply_balance_correction!(u, p_independent, time_dependent_cache) + + # Accumulate normalized Newton residual for convergence output + cache = integrator.cache + if hasproperty(cache, :nlsolver) + atmp = cache.nlsolver.cache.atmp + abs_atmp = abs.(atmp) + max_atmp = finitemaximum(abs_atmp; init = one(eltype(abs_atmp))) + for i in eachindex(p_independent.convergence) + v = abs_atmp[i] + p_independent.convergence[i] += isfinite(v) ? v / max_atmp : 0.0 end + p_independent.convergence_ncalls[1] += 1 end - convergence .+= - temp_convergence / - finitemaximum(temp_convergence; init = one(eltype(temp_convergence))) - ncalls[1] += 1 end + # Update cumulative boundary flow which is integrated exactly + @. flow_boundary.cumulative_flow_saveat += + time_dependent_cache.flow_boundary.current_cumulative_boundary_flow - + flow_boundary.cumulative_flow + @. flow_boundary.cumulative_flow = + time_dependent_cache.flow_boundary.current_cumulative_boundary_flow + # Update cumulative forcings which are integrated exactly @. basin.cumulative_drainage_saveat += time_dependent_cache.basin.current_cumulative_drainage - basin.cumulative_drainage @@ -162,18 +418,10 @@ function update_cumulative_flows!(u, t, integrator)::Nothing @. basin.cumulative_surface_runoff = time_dependent_cache.basin.current_cumulative_surface_runoff - # Update cumulative boundary flow which is integrated exactly - @. flow_boundary.cumulative_flow_saveat += - time_dependent_cache.flow_boundary.current_cumulative_boundary_flow - - flow_boundary.cumulative_flow - @. flow_boundary.cumulative_flow = - time_dependent_cache.flow_boundary.current_cumulative_boundary_flow - # Update supplied flows for allocation input and output for allocation_model in allocation.allocation_models (; cumulative_supplied_volume) = allocation_model - # Update supplied flows for allocation output for link in keys(cumulative_supplied_volume) cumulative_supplied_volume[link] += flow_update_on_link(integrator, link) end @@ -182,7 +430,7 @@ function update_cumulative_flows!(u, t, integrator)::Nothing end function update_concentrations!(u, t, integrator)::Nothing - (; uprev, p, tprev, dt) = integrator + (; p, tprev, dt) = integrator (; p_independent, state_and_time_dependent_cache) = p (; current_storage, current_level) = state_and_time_dependent_cache (; basin, flow_boundary, do_concentration) = p_independent @@ -236,7 +484,7 @@ function update_concentrations!(u, t, integrator)::Nothing add_substance_mass!( mass_node, loads_itp[node_id.idx], - dt, # loads are per second, not volume, so the flow is just the time step + dt, t, ) end @@ -281,13 +529,13 @@ function update_concentrations!(u, t, integrator)::Nothing for node_id in basin.node_id mass_node = mass[node_id.idx] - # Evaporate mass to keep the mass balance, if enabled in model config + # Evaporate mass using accumulated evaporation if evaporate_mass - evaporated_volume = u.evaporation[node_id.idx] - uprev.evaporation[node_id.idx] + evaporated_volume = p.state_and_time_dependent_cache.current_evaporation[node_id.idx] * dt mass_node .-= concentration_state[node_id.idx, :] .* evaporated_volume end - infiltrated_volume = u.infiltration[node_id.idx] - uprev.infiltration[node_id.idx] + infiltrated_volume = p.state_and_time_dependent_cache.current_infiltration[node_id.idx] * dt mass_node .-= concentration_state[node_id.idx, :] .* infiltrated_volume # Take care of infinitely small masses, possibly becoming negative due to truncation. @@ -328,8 +576,9 @@ end Compute the forcing volume entering and leaving the Basin over the last time step """ function forcing_update(integrator::DEIntegrator, node_id::NodeID)::Tuple{Float64, Float64} - (; u, uprev, p, dt) = integrator - (; basin) = p.p_independent + (; p, dt) = integrator + (; basin, p_independent) = (p.p_independent, p) + basin = p.p_independent.basin (; vertical_flux) = basin @assert node_id.type == NodeType.Basin @@ -343,23 +592,24 @@ function forcing_update(integrator::DEIntegrator, node_id::NodeID)::Tuple{Float6 vertical_flux.surface_runoff[node_id.idx] ) * dt + # Use trapezoidal accumulated evaporation/infiltration outflow_update = - (u.evaporation[node_id.idx] - uprev.evaporation[node_id.idx]) + - (u.infiltration[node_id.idx] - uprev.infiltration[node_id.idx]) + p.state_and_time_dependent_cache.current_evaporation[node_id.idx] * dt + + p.state_and_time_dependent_cache.current_infiltration[node_id.idx] * dt return inflow_update, outflow_update end """ -Given an link (from_id, to_id), compute the cumulative flow over that -link over the latest time step. +Given a link (from_id, to_id), compute the cumulative flow over that +link over the latest time step. Uses the trapezoidal-accumulated cumulative flows. """ function flow_update_on_link( integrator::DEIntegrator, link_src::Tuple{NodeID, NodeID}, )::Float64 - (; u, uprev, p, t, tprev) = integrator - (; flow_boundary, state_ranges, link_to_state_idx) = p.p_independent + (; p, t, tprev) = integrator + (; flow_boundary) = p.p_independent from_id, to_id = link_src return if from_id == to_id @@ -369,8 +619,16 @@ function flow_update_on_link( elseif from_id.type == NodeType.FlowBoundary integral(flow_boundary.flow_rate[from_id.idx], tprev, t) else - flow_idx = get_state_index(state_ranges, link_to_state_idx, link_src) - u[flow_idx] - uprev[flow_idx] + # Look up the internal flow link index and use cumulative flow + flow_link_lookup = p.p_independent.graph[].flow_link_lookup + link_idx = get_link_index(link_src, flow_link_lookup) + if isnothing(link_idx) + 0.0 + else + # Return cumulative flow since last update + # This gets reset in update_cumulative_flows! or save_flow + p.p_independent.cumulative_flow[link_idx] + end end end @@ -383,38 +641,33 @@ function save_basin_state(u, t, integrator) end """ -Save all cumulative forcings and flows over links over the latest timestep, -Both computed by the solver and integrated exactly. Also computes the total horizontal -inflow and outflow per Basin. +Save all flow rates (averaged over the saveat interval) and vertical fluxes. """ function save_flow(u, t, integrator) - (; cache, p) = integrator - (; - basin, - state_inflow_link, - state_outflow_link, - flow_boundary, - u_prev_saveat, - convergence, - ncalls, - node_id, - ) = p.p_independent + (; p) = integrator + (; p_independent) = p + (; basin, flow_boundary, graph) = p_independent + internal_flow_links = graph[].internal_flow_links + Δt = get_Δt(integrator) - flow_mean = (u - u_prev_saveat) / Δt - # Current u is previous u in next computation - u_prev_saveat .= u + # Compute mean flow rate per internal link from cumulative flows + n_links = length(internal_flow_links) + flow_mean = p_independent.cumulative_flow_saveat ./ Δt + + # Reset saveat accumulators + p_independent.cumulative_flow_saveat .= 0.0 n_basin = length(basin.node_id) inflow_mean = zeros(n_basin) outflow_mean = zeros(n_basin) - flow_convergence = fill(missing, length(u)) |> Vector{Union{Missing, Float64}} - basin_convergence = fill(missing, n_basin) |> Vector{Union{Missing, Float64}} - # Flow contributions from horizontal flow states - for (flow, inflow_link, outflow_link) in - zip(flow_mean, state_inflow_link, state_outflow_link) - inflow_id = inflow_link.link[1] + # Flow contributions from horizontal flow links + for (fi, link_meta) in enumerate(internal_flow_links) + flow = flow_mean[fi] + inflow_id = link_meta.link[1] + outflow_id = link_meta.link[2] + if inflow_id.type == NodeType.Basin if flow > 0 outflow_mean[inflow_id.idx] += flow @@ -423,7 +676,6 @@ function save_flow(u, t, integrator) end end - outflow_id = outflow_link.link[2] if outflow_id.type == NodeType.Basin if flow > 0 inflow_mean[outflow_id.idx] += flow @@ -445,30 +697,34 @@ function save_flow(u, t, integrator) end end + # Vertical fluxes from exact integration via time_dependent_cache precipitation = copy(basin.cumulative_precipitation_saveat) ./ Δt surface_runoff = copy(basin.cumulative_surface_runoff_saveat) ./ Δt drainage = copy(basin.cumulative_drainage_saveat) ./ Δt @. basin.cumulative_precipitation_saveat = 0.0 @. basin.cumulative_surface_runoff_saveat = 0.0 @. basin.cumulative_drainage_saveat = 0.0 + evaporation = p_independent.cumulative_evaporation_saveat ./ Δt + infiltration = p_independent.cumulative_infiltration_saveat ./ Δt - if hasproperty(cache, :nlsolver) - flow_convergence = convergence ./ ncalls[1] - for (i, (evap, infil)) in - enumerate(zip(flow_convergence.evaporation, flow_convergence.infiltration)) - if isnan(evap) - basin_convergence[i] = infil - elseif isnan(infil) - basin_convergence[i] = evap - else - basin_convergence[i] = max(evap, infil) - end + # Reset saveat accumulators for evaporation/infiltration + p_independent.cumulative_evaporation_saveat .= 0.0 + p_independent.cumulative_infiltration_saveat .= 0.0 + + concentration = copy(basin.concentration_data.concentration_state) + + # Compute mean convergence over the saveat interval (missing if no nlsolver calls) + n_basin = length(basin.node_id) + convergence = fill(missing, n_basin) |> Vector{Union{Missing, Float64}} + ncalls = p_independent.convergence_ncalls[1] + if ncalls > 0 + for i in 1:n_basin + convergence[i] = p_independent.convergence[i] / ncalls end - fill!(convergence, 0) - ncalls[1] = 0 + fill!(p_independent.convergence, 0.0) + p_independent.convergence_ncalls[1] = 0 end - concentration = copy(basin.concentration_data.concentration_state) saved_flow = SavedFlow(; flow = flow_mean, inflow = inflow_mean, @@ -477,11 +733,14 @@ function save_flow(u, t, integrator) precipitation, surface_runoff, drainage, + evaporation, + infiltration, concentration, - flow_convergence, - basin_convergence, + convergence, + flow_convergence = copy(p_independent.flow_convergence_saveat) ./ Δt, t, ) + p_independent.flow_convergence_saveat .= 0.0 check_water_balance_error!(saved_flow, integrator, Δt) return saved_flow end @@ -493,20 +752,11 @@ function check_water_balance_error!( )::Nothing (; u, p, t) = integrator (; p_independent, state_and_time_dependent_cache) = p - (; u_reduced, state_ranges) = p_independent (; current_storage) = state_and_time_dependent_cache (; basin, water_balance_abstol, water_balance_reltol, starttime) = p_independent errors = false - # The initial storage is irrelevant for the storage rate and can only cause - # floating point truncation errors - reduce_state!(u_reduced, u, p_independent) - formulate_storages!(u_reduced, p, t; add_initial_storage = false) - - evaporation = view(saved_flow.flow, state_ranges.evaporation) - infiltration = view(saved_flow.flow, state_ranges.infiltration) - for ( inflow_rate, outflow_rate, @@ -524,8 +774,8 @@ function check_water_balance_error!( saved_flow.precipitation, saved_flow.surface_runoff, saved_flow.drainage, - evaporation, - infiltration, + saved_flow.evaporation, + saved_flow.infiltration, current_storage, basin.Δstorage_prev_saveat, basin.node_id, @@ -553,7 +803,6 @@ function check_water_balance_error!( end @. basin.Δstorage_prev_saveat = current_storage - current_storage .+= basin.storage0 return nothing end @@ -573,14 +822,14 @@ end function check_negative_storage(u, t, integrator)::Nothing (; p) = integrator - (; p_independent, state_and_time_dependent_cache) = p + (; p_independent) = p (; basin) = p_independent du = get_du(integrator) water_balance!(du, u, p, t) errors = false for id in basin.node_id - if state_and_time_dependent_cache.current_storage[id.idx] < 0 + if u.basin[id.idx] < 0 @error "Negative storage detected in $id" errors = true end @@ -609,7 +858,6 @@ function apply_discrete_control!(u, t, integrator)::Nothing (; p) = integrator (; discrete_control) = p.p_independent (; node_id, truth_state, compound_variables) = discrete_control - du = get_du(integrator) # Loop over the discrete control nodes to determine their truth state # and detect possible control state changes @@ -625,7 +873,7 @@ function apply_discrete_control!(u, t, integrator)::Nothing # Loop over the compound variables listened to by this discrete control node for compound_variable in compound_variables_node - value = compound_variable_value(compound_variable, p, du, t) + value = compound_variable_value(compound_variable, p, t) # Loop over the threshold interpolations associated with the current compound variable for (threshold_low, threshold_high) in @@ -703,6 +951,26 @@ function set_new_control_state!( tabulated_rating_curve.allocation_controlled[target_node_id.idx] = control_state_update.allocation_controlled end + + # Mark the links of this node so that the next trapezoidal integration step + # uses backward Euler (avoids a spike from averaging pre/post-switch rates). + flow_link_lookup = p_independent.graph[].flow_link_lookup + idx = target_node_id.idx + node_links = if target_node_id.type == NodeType.Pump + (pump.inflow_link[idx].link, pump.outflow_link[idx].link) + elseif target_node_id.type == NodeType.Outlet + (outlet.inflow_link[idx].link, outlet.outflow_link[idx].link) + elseif target_node_id.type == NodeType.TabulatedRatingCurve + (tabulated_rating_curve.inflow_link[idx].link, tabulated_rating_curve.outflow_link[idx].link) + else + nothing + end + if node_links !== nothing + for link in node_links + link_idx = get_link_index(link, flow_link_lookup) + link_idx !== nothing && (p_independent.flow_rate_prev[link_idx] = NaN) + end + end end discrete_control.control_state[discrete_control_id.idx] = control_state_new @@ -715,12 +983,12 @@ end Get a value for a condition. Currently supports getting levels from Basins and flows from FlowBoundaries. """ -function get_value(subvariable::SubVariable, p::Parameters, du::CVector, t::Float64) +function get_value(subvariable::SubVariable, p::Parameters, t::Real) (; flow_boundary, level_boundary, basin) = p.p_independent (; listen_node_id, look_ahead, variable, cache_ref) = subvariable if !iszero(cache_ref.idx) - return get_value(cache_ref, p, du) + return get_value(cache_ref, p) end if variable == "level" @@ -754,10 +1022,10 @@ function get_value(subvariable::SubVariable, p::Parameters, du::CVector, t::Floa return value end -function compound_variable_value(compound_variable::CompoundVariable, p, du, t) - value = zero(eltype(du)) +function compound_variable_value(compound_variable::CompoundVariable, p, t) + value = zero(typeof(t)) for subvariable in compound_variable.subvariables - value += subvariable.weight * get_value(subvariable, p, du, t) + value += subvariable.weight * get_value(subvariable, p, t) end return value end diff --git a/core/src/concentration.jl b/core/src/concentration.jl index 5c91d370d..a30b86b9f 100644 --- a/core/src/concentration.jl +++ b/core/src/concentration.jl @@ -3,16 +3,18 @@ Process mass inflows from UserDemand separately as the inflow and outflow are decoupled in the states """ function mass_inflows_from_user_demand!(integrator::DEIntegrator)::Nothing - (; p, t) = integrator + (; p, t, dt) = integrator (; basin, user_demand) = p.p_independent (; concentration_state, mass) = basin.concentration_data + (; state_and_time_dependent_cache) = p for (node_idx, outflow_link) in enumerate(user_demand.outflow_link) to_node = outflow_link.link[2] user_demand_idx = outflow_link.link[1].idx inflow_links = user_demand.inflow_links[node_idx] - cumulative_user_demand_outflow = flow_update_on_link(integrator, outflow_link.link) + # Use current flow rate * dt as approximate cumulative flow + cumulative_user_demand_outflow = state_and_time_dependent_cache.current_flow_rate_user_demand[user_demand_idx] * dt if to_node.type == NodeType.Basin # Mix concentrations of all inflow links weighted by each link's cumulative @@ -38,7 +40,6 @@ function mass_inflows_from_user_demand!(integrator::DEIntegrator)::Nothing end mass[to_node.idx][Substance.UserDemand] = ud_mass_before - # Add fresh UserDemand tracer (= 1.0) and any user-defined substances add_substance_mass!( mass[to_node.idx], user_demand.concentration_itp[user_demand_idx], @@ -54,67 +55,52 @@ end Process all mass inflows to basins """ function mass_inflows_basin!(integrator::DEIntegrator)::Nothing - (; p, t) = integrator - (; basin, state_inflow_link, state_outflow_link, level_boundary) = p.p_independent + (; p, t, dt) = integrator + (; basin, level_boundary, graph) = p.p_independent (; cumulative_in, concentration_state, mass) = basin.concentration_data + internal_flow_links = graph[].internal_flow_links - # Loop over connections that have state - for (inflow_link, outflow_link) in zip(state_inflow_link, state_outflow_link) - from_node = inflow_link.link[1] - state_node = inflow_link.link[2] - to_node = outflow_link.link[2] + # Loop over internal flow links + for (fi, link_meta) in enumerate(internal_flow_links) + from_node = link_meta.link[1] + to_node = link_meta.link[2] - if state_node.type == NodeType.UserDemand - # UserDemand is handled separately in mass_inflows_from_user_demand + # Skip UserDemand - handled separately + if from_node.type == NodeType.UserDemand || to_node.type == NodeType.UserDemand continue end - if from_node.type == NodeType.Basin - cumulative_flow = flow_update_on_link(integrator, inflow_link.link) - # Negative flow over the inflow link means flow into the from_node - if cumulative_flow < 0 - cumulative_in[from_node.idx] -= cumulative_flow - if to_node.type == NodeType.Basin - mass[from_node.idx] .-= - concentration_state[to_node.idx, :] .* cumulative_flow - elseif to_node.type == NodeType.LevelBoundary - add_substance_mass!( - mass[from_node.idx], - level_boundary.concentration_itp[to_node.idx], - -cumulative_flow, - t, - ) - elseif (to_node.type == NodeType.Terminal && to_node.value == 0) - # UserDemand inflow is discoupled from its outflow - # The unset flow link defaults to Terminal #0 - nothing - else - @warn "Unsupported outflow from $to_node to $from_node with cumulative flow $cumulative_flow m³" - end + # Use cumulative flow from trapezoidal integration + cumulative_flow = p.p_independent.cumulative_flow[fi] + + if from_node.type == NodeType.Basin && cumulative_flow < 0 + # Negative flow means flow into from_node + cumulative_in[from_node.idx] -= cumulative_flow + if to_node.type == NodeType.Basin + mass[from_node.idx] .-= + concentration_state[to_node.idx, :] .* cumulative_flow + elseif to_node.type == NodeType.LevelBoundary + add_substance_mass!( + mass[from_node.idx], + level_boundary.concentration_itp[to_node.idx], + -cumulative_flow, + t, + ) end end - if to_node.type == NodeType.Basin - cumulative_flow = flow_update_on_link(integrator, outflow_link.link) - if cumulative_flow > 0 - cumulative_in[to_node.idx] += cumulative_flow - if from_node.type == NodeType.Basin - mass[to_node.idx] .+= - concentration_state[from_node.idx, :] .* cumulative_flow - - elseif from_node.type == NodeType.LevelBoundary - add_substance_mass!( - mass[to_node.idx], - level_boundary.concentration_itp[from_node.idx], - cumulative_flow, - t, - ) - elseif from_node.type == NodeType.Terminal && from_node.value == 0 - # The unset flow link defaults to Terminal #0 - nothing - else - @warn "Unsupported outflow from $from_node to $to_node with flow $cumulative_flow m³" - end + if to_node.type == NodeType.Basin && cumulative_flow > 0 + cumulative_in[to_node.idx] += cumulative_flow + if from_node.type == NodeType.Basin + mass[to_node.idx] .+= + concentration_state[from_node.idx, :] .* cumulative_flow + elseif from_node.type == NodeType.LevelBoundary + add_substance_mass!( + mass[to_node.idx], + level_boundary.concentration_itp[from_node.idx], + cumulative_flow, + t, + ) end end end @@ -125,24 +111,21 @@ end Process all mass outflows from Basins """ function mass_outflows_basin!(integrator::DEIntegrator)::Nothing - (; state_inflow_link, state_outflow_link, basin) = integrator.p.p_independent + (; graph, basin) = integrator.p.p_independent (; mass, concentration_state) = basin.concentration_data + internal_flow_links = graph[].internal_flow_links - @views for (inflow_link, outflow_link) in zip(state_inflow_link, state_outflow_link) - from_node = inflow_link.link[1] - to_node = outflow_link.link[2] + @views for (fi, link_meta) in enumerate(internal_flow_links) + from_node = link_meta.link[1] + to_node = link_meta.link[2] - if from_node.type == NodeType.Basin - flow = flow_update_on_link(integrator, inflow_link.link) - if flow > 0 - mass[from_node.idx] .-= concentration_state[from_node.idx, :] .* flow - end + cumulative_flow = integrator.p.p_independent.cumulative_flow[fi] + + if from_node.type == NodeType.Basin && cumulative_flow > 0 + mass[from_node.idx] .-= concentration_state[from_node.idx, :] .* cumulative_flow end - if to_node.type == NodeType.Basin - flow = flow_update_on_link(integrator, outflow_link.link) - if flow < 0 - mass[to_node.idx] .+= concentration_state[to_node.idx, :] .* flow - end + if to_node.type == NodeType.Basin && cumulative_flow < 0 + mass[to_node.idx] .+= concentration_state[to_node.idx, :] .* cumulative_flow end end return nothing diff --git a/core/src/config.jl b/core/src/config.jl index 502c173a7..c188c7b5f 100644 --- a/core/src/config.jl +++ b/core/src/config.jl @@ -22,6 +22,7 @@ using OrdinaryDiffEqSDIRK: ImplicitEuler, KenCarp4, TRBDF2 using OrdinaryDiffEqBDF: FBDF, QNDF using OrdinaryDiffEqRosenbrock: Rosenbrock23, Rodas4P, Rodas5P import OrdinaryDiffEqDifferentiation +using PositiveIntegrators: MPRK22, MPRK43I, MPRK43II using LinearSolve: KLUFactorization, SciMLLinearSolveAlgorithm, LinearSolve, SciMLLinearSolveAlgorithm @@ -29,6 +30,7 @@ export Config, Solver, Results, Logging, Toml export algorithm, camel_case, get_ad_type, + is_mprk_algorithm, snake_case, input_path, database_path, @@ -180,7 +182,7 @@ end dtmin::Float64 = 0.0 dtmax::Union{Float64, Nothing} = nothing force_dtmin::Bool = false - abstol::Float64 = 1.0e-5 + abstol::Float64 = 1.0e-4 reltol::Float64 = 1.0e-5 water_balance_abstol::Float64 = 1.0e-3 water_balance_reltol::Float64 = 1.0e-2 @@ -288,8 +290,8 @@ to avoid runtime errors, especially when writing results. function validate_config(toml::Toml)::Nothing is_valid = true - if !haskey(algorithms, toml.solver.algorithm) - options = join(keys(algorithms), ", ") + if !haskey(algorithms, toml.solver.algorithm) && toml.solver.algorithm ∉ mprk_algorithms + options = join(union(keys(algorithms), mprk_algorithms), ", ") @error("Given solver algorithm $(toml.solver.algorithm) not supported.\n\ Available options are: ($(options)).") is_valid = false @@ -375,6 +377,9 @@ const algorithms = Dict{String, Type}( "Euler" => Euler, ) +"MPRK algorithms requiring PDSProblem formulation" +const mprk_algorithms = Set(["MPRK22", "MPRK43I", "MPRK43II"]) + """ Check whether the given function has a method that accepts the given kwarg. Note that it is possible that methods exist that accept :a and :b individually, @@ -388,41 +393,32 @@ function function_accepts_kwarg(f, kwarg)::Bool end function get_ad_type(solver::Solver) - chunksize = solver.specialize ? nothing : 1 return if solver.autodiff - AutoForwardDiff(; chunksize, tag = :Ribasim) + AutoForwardDiff(; tag = :Ribasim) else AutoFiniteDiff() end end -""" -A wrapper of a SciMLLinearSolveAlgorithm to dispatch on for the specialized Jacobian -matrix of Ribasim. -""" -struct RibasimLinearSolve{AType <: SciMLLinearSolveAlgorithm} <: SciMLLinearSolveAlgorithm - algorithm::AType -end +"Whether the configured algorithm is a positivity-preserving MPRK method" +is_mprk_algorithm(solver::Solver)::Bool = solver.algorithm in mprk_algorithms -LinearSolve.needs_concrete_A(::RibasimLinearSolve) = false +"Create an algorithm from solver config" +function algorithm(solver::Solver) + if is_mprk_algorithm(solver) + return mprk_algorithm(solver) + end -"Create an OrdinaryDiffEqAlgorithm from solver config" -function algorithm(solver::Solver)::OrdinaryDiffEqAlgorithm kwargs = Dict{Symbol, Any}() algotype = algorithms[solver.algorithm] if algotype <: OrdinaryDiffEqNewtonAdaptiveAlgorithm - kwargs[:nlsolve] = NLNewton() + kwargs[:nlsolve] = NLNewton(; relax = Ribasim.BackTracking()) if solver.sparse - kwargs[:linsolve] = - RibasimLinearSolve(KLUFactorization(; check_pattern = false)) + kwargs[:linsolve] = KLUFactorization(; check_pattern = false) end end - if function_accepts_kwarg(algotype, :step_limiter!) - kwargs[:step_limiter!] = Ribasim.limit_flow! - end - if function_accepts_kwarg(algotype, :autodiff) kwargs[:autodiff] = get_ad_type(solver) end @@ -430,6 +426,22 @@ function algorithm(solver::Solver)::OrdinaryDiffEqAlgorithm return algotype(; kwargs...) end +"Create an MPRK algorithm from solver config" +function mprk_algorithm(solver::Solver) + linsolve = solver.sparse ? KLUFactorization(; check_pattern = false) : nothing + kwargs = isnothing(linsolve) ? (;) : (; linsolve) + + return if solver.algorithm == "MPRK22" + MPRK22(1.0; kwargs...) + elseif solver.algorithm == "MPRK43I" + MPRK43I(1.0, 0.5; kwargs...) + elseif solver.algorithm == "MPRK43II" + MPRK43II(2.0 / 3.0; kwargs...) + else + error("Unknown MPRK algorithm: $(solver.algorithm)") + end +end + "Convert the saveat Float64 from our Config to SciML's saveat" function convert_saveat(saveat::Float64, t_end::Float64)::Union{Float64, Vector{Float64}} errors = false diff --git a/core/src/differentiation.jl b/core/src/differentiation.jl deleted file mode 100644 index afa9df477..000000000 --- a/core/src/differentiation.jl +++ /dev/null @@ -1,496 +0,0 @@ -#= -Theoretical background ----------------------- - -The ODE system that is solved by Ribasim is formulated in terms of cumulative flow states for -reasons of water balance accuracy, see https://github.com/Deltares/Ribasim/pull/1819. The state vector -`u` thus has its various components for the connector node flows, and also the PID control integral terms. - -In the RHS of the ODE systems these states are summed to obtain their contribution to -the Basin storages, and this is the only way these states are used in the RHS. The RHS is thus of the form -`f(u) = g(A*u)` where `A` is a highly sparse matrix with the following structure: - - ⎡ max one 1 and ⎢ ⎢ ⎢ ⠀⎤ - ⎢ one -1 per ⎢ -I ⎢ -I ⎢ 0 ⎢ Basin rows -A = ⎢ col by graph ⎢ evap. ⎢ infl. ⎢ ⎢ - ⎢---------------⎢---------⎢---------⎢---⎢ - ⎣ 0 ⎢ 0 ⎢ 0 ⎢ I⠀⎦ PID integral rows - -The upper left section of `A` depends on the graph of the model. The evaporation and infiltration sections are simply -identity matrices multiplied by -1, subtracting these fluxes from the corresponding Basin Storages. -We call `u_reduced = A*u`. The matrix `A` is never explicitly constructed in the code (not even as a `SparseMatrixCSC`). -Instead, the matrix-vector multiplication `A*u` is defined in the function `reduce_state!`. - -This structure of the RHS can be taken advantage of in several ways. First of all, The Jacobian of the RHS can be expressed -as a matrix-matrix product using the chain rule: - -`Jf(u) = Jg(A*u)*A.` - -Note that the Jacobian of `g` has `length(u)` rows and `length(u_reduced)` columns. It turns out that taking this multiplication -with `A` out of the AD Jacobian computation is very advantageous in terms of computation time. In the code `f` and `g` are both methods -of `water_balance!`. We denote the Jacobian of `g` as `J_intermediate`. The custom Jacobian object that wraps `J_intermediate` is called -`HalfLazyJacobian`, which contains the parameters to implicitly define `A` and other parameters needed for AD Jacobian evaluation. - -It turns out that we can make use of this structure of the Jacobian in the linear solve as well. The linear system is of the form -`Wa = b`, where `W = -γ⁻¹I - J_intermediate*A`, as explained in https://book.sciml.ai/notes/09-Solving_Stiff_Ordinary_Differential_Equations/. -It turns out this special form of the Jacobian can be utilized to solve the linear system in 2 -steps: -- Solve `(-γ⁻¹I - A * J_intermediate)*c = A*b` for `c`, -- Compute `a = -γ * (b + J_intermediate * c)`. - -For more details on the derivation of this see https://github.com/Deltares/Ribasim/pull/2624#issuecomment-3382550431. -The crucial detail here is that the linear system to be solved has much fewer equations. We denote `J_inner = A * J_intermediate`, -which is explicitly computed as a sparse matrix in `calc_J_inner`. The above computation is incorporated in the ODE solve as follows: - -- A custom linear solve algorithm `RibasimLinearSolve` wraps a default linear solve algorithm, used for sparse solves -- `SciMLBase.init` is overloaded for `RibasimLinearSolve`, initializing the cache for the inner solve (i.e. the solve in 'storage space') - and returns this in the `RibasimLinearSolveCache` -- `OrdinaryDiffEqDifferentiation.dolinsolve` is overloaded for `RibasimLinearSolveCache`, which - transforms the problem to 'storage space', solves the problem, and translates the result back to 'flow space' -=# - -import SparseArrays - -""" -Helper to construct a WOperator from an ODEFunction, replicating the behavior -of the removed constructor from OrdinaryDiffEqDifferentiation, -see https://github.com/SciML/OrdinaryDiffEq.jl/pull/3017 -""" -function make_woperator(f, u, gamma) - mass_matrix = f.mass_matrix - if !isa(mass_matrix, Union{AbstractMatrix, UniformScaling}) - mass_matrix = convert(AbstractMatrix, mass_matrix) - end - J = deepcopy(f.jac_prototype) - if J isa AbstractMatrix - @assert SciMLBase.has_jac(f) "f needs to have an associated jacobian" - J = MatrixOperator(J; update_func! = f.jac) - end - return WOperator{true}(mass_matrix, gamma, J, u) -end - -""" -The HalfLazyJacobian represents the Ribasim Jacobian in the form `J = J_intermediate * A` -(see also the theoretical background in differentiation.jl). -`J_intermediate` is explicitly (AD) computed, and `A` is implicit in: -- `reduce_state!`, which defines the matrix-vector product `u_reduced = A * u`; -- `calc_J_inner!`, which defined the matrix-matrix product `J_inner = A * J_intermediate`. -""" -struct HalfLazyJacobian <: AbstractSciMLOperator{Float64} - J_intermediate::SparseMatrixCSC{Float64, Int64} - p_independent::ParametersIndependent - du::CVector - prep::Any - backend::Any -end - -# Used in the default GMRES linear solve for -# dense Jacobians -function LinearAlgebra.mul!( - _u::RibasimCVectorType, - J::HalfLazyJacobian, - _v::RibasimCVectorType, - ) - (; J_intermediate, p_independent) = J - (; u_reduced, state_ranges) = p_independent - # The input vectors are rewrapped because somewhere - # they obtain wrong axes - u = CVector(getdata(_u), state_ranges) - v = CVector(getdata(_v), state_ranges) - reduce_state!(u_reduced, v, p_independent) - mul!(u, J_intermediate, u_reduced) - return nothing -end - -# SciMLOperators interface -SciMLOperators.isconstant(::HalfLazyJacobian) = false -SciMLOperators.issquare(::HalfLazyJacobian) = true -SciMLOperators.islinear(::HalfLazyJacobian) = true -SciMLOperators.isconvertible(::HalfLazyJacobian) = false -SciMLOperators.has_mul!(::HalfLazyJacobian) = true - -SciMLBase.update_coefficients!(J::HalfLazyJacobian, u, p, t) = - get_jacobian!(J, J.du, u, p, t, J.prep, J.backend) - -# Overloads to make OrdinaryDiffEq happy -Base.size(J::HalfLazyJacobian) = (length(J.du), length(J.du)) -ADTypes.KnownJacobianSparsityDetector(J::HalfLazyJacobian) = - ADTypes.KnownJacobianSparsityDetector(J.J_intermediate) - -""" -The cache associated with the custom linear solve algorithm `config.RibasimLinearSolve`. -""" -struct RibasimLinearSolveCache{C, WType} - cache_inner::C - W::WType -end - -""" -Compute the product `J_inner = A * J_intermediate`, where `A` is implicitly defined -by the structure of the Ribasim model. -""" -function calc_J_inner!( - J_inner::AbstractMatrix, - J::HalfLazyJacobian - )::Nothing - J_inner .= 0 - n_states_reduced = size(J_inner)[1] - - for col_reduced in 1:n_states_reduced - update_J_inner!(J_inner, J, col_reduced) - end - return nothing -end - -function update_J_inner!( - J_inner::SparseMatrixCSC, - J::HalfLazyJacobian, - col_reduced::Int, - )::Nothing - (; J_intermediate, p_independent) = J - for nz_idx in nzrange(J_intermediate, col_reduced) - row = J_intermediate.rowval[nz_idx] - val = J_intermediate.nzval[nz_idx] - update_J_inner!(J_inner, p_independent, row, col_reduced, val) - end - return -end - -function update_J_inner!(J_inner::Matrix, J::HalfLazyJacobian, col_reduced::Int)::Nothing - (; J_intermediate, p_independent) = J - for row in 1:size(J_intermediate)[2] - val = J_inner[row, col_reduced] - !iszero(val) && update_J_inner!(J_inner, p_independent, row, col_reduced, val) - end - return -end - -function update_J_inner!( - J_inner::AbstractMatrix, - p_independent::ParametersIndependent, - row::Int, - col_reduced::Int, - val::Float64, - ) - (; - tabulated_rating_curve, - pump, - outlet, - user_demand, - linear_resistance, - manning_resistance, - basin, - state_ranges, - ) = p_independent - node_id = p_independent.node_id[row] - - return if row in state_ranges.tabulated_rating_curve - update_J_inner!(J_inner, val, node_id, col_reduced, tabulated_rating_curve) - elseif row in state_ranges.pump - update_J_inner!(J_inner, val, node_id, col_reduced, pump) - elseif row in state_ranges.outlet - update_J_inner!(J_inner, val, node_id, col_reduced, outlet) - elseif row in state_ranges.user_demand_inflow - inflow_link_meta = p_independent.state_inflow_link[row] - inflow_id = inflow_link_meta.link[1] - if inflow_id.type == NodeType.Basin - J_inner[inflow_id.idx, col_reduced] -= val - end - elseif row in state_ranges.user_demand_outflow - update_J_inner!(J_inner, val, node_id, col_reduced, user_demand; do_inflow = false) - elseif row in state_ranges.linear_resistance - update_J_inner!(J_inner, val, node_id, col_reduced, linear_resistance) - elseif row in state_ranges.manning_resistance - update_J_inner!(J_inner, val, node_id, col_reduced, manning_resistance) - elseif row in state_ranges.evaporation - @assert node_id.type == NodeType.Basin - @assert node_id.idx == col_reduced - J_inner[node_id.idx, col_reduced] -= val - elseif row in state_ranges.infiltration - @assert node_id.type == NodeType.Basin - @assert node_id.idx == col_reduced - J_inner[node_id.idx, col_reduced] -= val - else - @assert node_id.type == NodeType.PidControl - @assert row in state_ranges.integral - row_reduced = length(basin.node_id) + (row - state_ranges.integral.start + 1) - J_inner[row_reduced, col_reduced] += val - end -end - -function update_J_inner!( - J_inner::AbstractMatrix, - val::Float64, - node_id::NodeID, - col_reduced::Int, - node::AbstractParameterNode; - do_inflow::Bool = true, - )::Nothing - if do_inflow - inflow_id = node.inflow_link[node_id.idx].link[1] - if inflow_id.type == NodeType.Basin - J_inner[inflow_id.idx, col_reduced] -= val - end - end - - outflow_id = node.outflow_link[node_id.idx].link[2] - if outflow_id.type == NodeType.Basin - J_inner[outflow_id.idx, col_reduced] += val - end - return nothing -end - -""" -Calculate `W`, which is the matrix in the linear solve of the ODE solve algorithm. -""" -function calc_W_inner!(W_inner, J)::Nothing - calc_J_inner!(W_inner.J.A, J) - jacobian2W!(W_inner._concrete_form, W_inner.mass_matrix, W_inner.gamma, W_inner.J.A) - return nothing -end - -""" -Initialize the `RibasimLinearSolveCache` for the `config.RibasimLinearSolve` algorithm. -This cache contains `cache_inner` for the actual linear solve in the reduced state space, and `W` -from the original state space which directly interacts with `OrdinaryDiffEqNonlinearSolve.jl`. -""" -function SciMLBase.init( - prob::LinearProblem, - alg::config.RibasimLinearSolve, - args...; - kwargs..., - ) - W = prob.A - (; J) = W - (; u_reduced) = J.p_independent - n_states_reduced = length(u_reduced) - J_inner = similar(J.J_intermediate, (n_states_reduced, n_states_reduced)) - - # In this first call memory is allocated for the non zeros in the sparse case - calc_J_inner!(J_inner, J) - - W_inner = make_woperator( - ODEFunction(Returns(nothing); jac_prototype = J_inner, jac = Returns(nothing)), - u_reduced, - 1.0, - ) - - b_inner = copy(u_reduced) - prob_inner = LinearProblem(W_inner, b_inner) - cache_inner = init(prob_inner, alg.algorithm, args...; kwargs...) - - return RibasimLinearSolveCache(cache_inner, W) -end - -""" -This is a wrapper of the standard method of `dolinsolve`. It performs the transformations -between the original state space and the reduced state space and the linear solve in the -reduced state space. -""" -function OrdinaryDiffEqDifferentiation.dolinsolve( - integrator, - linsolve::RibasimLinearSolveCache; - b = nothing, - linu = nothing, - kwargs..., - ) - @assert !isnothing(b) - @assert !isnothing(linu) - (; cache_inner, W) = linsolve - (; J) = W - (; J_intermediate, p_independent) = J - γ = W.gamma - - W_inner = cache_inner.A - W_inner.gamma = γ - - # Translate the problem to the reduced state space - reduce_state!(cache_inner.b, b, p_independent) - calc_W_inner!(cache_inner.A, J) - cache_inner.isfresh = true - - # Solve the problem in the reduced state space - linres = dolinsolve( - integrator, - cache_inner; - kwargs..., - A = nothing, - linu = nothing, - b = nothing, - ) - - # Translate the solution back to the full state space - mul!(linu, J_intermediate, linres.u) - linu .-= b - linu .*= γ - - # Build new solution object - return LinearSolution{ - Float64, - 1, - RibasimCVectorType{Float64}, - typeof(linres.resid), - typeof(linres.alg), - typeof(linsolve), - typeof(linres.stats), - }( - linu, - linres.resid, - linres.alg, - linres.retcode, - linres.iters, - linsolve, - linres.stats, - ) -end - -function get_jacobian!(J::HalfLazyJacobian, du, u, p, t, prep, backend) - (; J_intermediate) = J - (; u_reduced) = p.p_independent - reduce_state!(u_reduced, u, p.p_independent) - - saved_td_t_prev = p.time_dependent_cache.t_prev_call[1] - # Invalidate t_prev_call so the first AD call's check_new_input! always sees t != -1, - # Otherwise, it would read garbage values from p.state_and_time_dependent_cache, - # since it is marked as Cache(), which means it starts as uninitialised dual number arrays. - p.time_dependent_cache.t_prev_call[1] = -1 - - jacobian!( - water_balance!, - du, - J_intermediate, - prep, - backend, - u_reduced, - Constant(p.p_independent), - Cache(p.state_and_time_dependent_cache), - Constant(p.time_dependent_cache), - Constant(p.p_mutable), - Constant(t), - ) - - # Restore shared state so next real RHS call works correctly - p.time_dependent_cache.t_prev_call[1] = saved_td_t_prev - - return J -end - -""" -Get the Jacobian evaluation function via DifferentiationInterface.jl. -The time derivative is also supplied in case a Rosenbrock method is used. -""" -function get_diff_eval(du::CVector, p::Parameters, solver::Solver) - (; p_independent, state_and_time_dependent_cache, time_dependent_cache, p_mutable) = p - (; u_reduced) = p_independent - backend = get_ad_type(solver) - sparsity_detector = TracerSparsityDetector() - # Use non-zero u to avoid missing connections in the sparsity - u_reduced_ = copy(u_reduced) - u_reduced_ .= 1 - - backend_jac = if solver.sparse - AutoSparse(backend; sparsity_detector, coloring_algorithm = GreedyColoringAlgorithm()) - else - backend - end - - t = 0.0 - - jac_prep = prepare_jacobian( - water_balance!, - du, - backend_jac, - u_reduced_, - Constant(p_independent), - Cache(state_and_time_dependent_cache), - Constant(time_dependent_cache), - Constant(p_mutable), - Constant(t); - strict = Val(true), - ) - - J_intermediate = - solver.sparse ? Float64.(sparsity_pattern(jac_prep)) : - zeros(length(du), length(u_reduced)) - jac_prototype = - HalfLazyJacobian(J_intermediate, p_independent, copy(du), jac_prep, backend_jac) - W_prototype = make_woperator( - ODEFunction(water_balance!; jac_prototype, jac = Returns(nothing)), - copy(du), - 0.0, - ) - - jac(J, u, p, t) = get_jacobian!(J, du, u, p, t, jac_prep, backend_jac) - - tgrad_prep = prepare_derivative( - water_balance!, - du, - backend, - t, - Constant(copy(du)), - Constant(p_independent), - Cache(state_and_time_dependent_cache), - Cache(time_dependent_cache), - Constant(p_mutable); - strict = Val(true), - ) - tgrad(dT, u, p, t) = derivative!( - water_balance!, - du, - dT, - tgrad_prep, - backend, - t, - Constant(u), - Constant(p.p_independent), - Cache(state_and_time_dependent_cache), - Cache(time_dependent_cache), - Constant(p.p_mutable), - ) - - time_dependent_cache.t_prev_call[1] = -1.0 - - return (; jac_prototype, W_prototype, jac, tgrad) -end - -# Method with `t` as second argument parsable by DifferentiationInterface.jl for time derivative computation -water_balance!( - du::RibasimCVectorType, - t::Number, - u::RibasimCVectorType, - p_independent::ParametersIndependent, - state_and_time_dependent_cache::StateAndTimeDependentCache, - time_dependent_cache::TimeDependentCache, - p_mutable::ParametersMutable, -) = water_balance!( - du, - u, - p_independent, - state_and_time_dependent_cache, - time_dependent_cache, - p_mutable, - t, -) - -# Method with `u` as second argument parsable by DifferentiationInterface.jl for Jacobian computation -function water_balance!( - du::RibasimCVectorType, - u::RibasimCVectorType, - p_independent::ParametersIndependent, - state_and_time_dependent_cache::StateAndTimeDependentCache, - time_dependent_cache::TimeDependentCache, - p_mutable::ParametersMutable, - t::Number, - )::Nothing - (; u_reduced) = p_independent - reduce_state!(u_reduced, u, p_independent) - return water_balance!( - du, - u_reduced, - p_independent, - state_and_time_dependent_cache, - time_dependent_cache, - p_mutable, - t, - ) -end diff --git a/core/src/graph.jl b/core/src/graph.jl index f92409e38..f72d3bf0e 100644 --- a/core/src/graph.jl +++ b/core/src/graph.jl @@ -131,6 +131,9 @@ function create_graph(db::DB, config::Config)::MetaGraph internal_flow_links, external_flow_links, flow_link_map, + flow_link_lookup = Dict{Tuple{NodeID, NodeID}, Int}( + link_meta.link => i for (i, link_meta) in enumerate(internal_flow_links) + ), ) @reset graph.graph_data = graph_data @@ -342,19 +345,18 @@ function inflow_id(graph::MetaGraph, id::NodeID)::NodeID end """ -Get the specific q from the input vector `flow` which has the same components as -the state vector, given an link (inflow_id, outflow_id). -`flow` can be either instantaneous or integrated/averaged. Instantaneous FlowBoundary flows can be obtained -from the parameters, but integrated/averaged FlowBoundary flows must be provided via `boundary_flow`. +Get the flow rate for a given link. +For internal flow links, looks up in the flow vector by link index. +For flow boundaries, returns the interpolated boundary flow rate. """ function get_flow( - flow::CVector, + flow::Vector{Float64}, p_independent::ParametersIndependent, t::Number, link::Tuple{NodeID, NodeID}; boundary_flow = nothing, ) - (; flow_boundary, state_ranges, link_to_state_idx) = p_independent + (; flow_boundary, graph) = p_independent from_id = link[1] return if from_id.type == NodeType.FlowBoundary if boundary_flow === nothing @@ -363,25 +365,9 @@ function get_flow( boundary_flow[from_id.idx] end else - flow[get_state_index(state_ranges, link_to_state_idx, link)] - end -end - -""" -Like `get_flow` and `get_state_index`, but for convergence, so without the boundary flow. -""" -function get_convergence( - convergence::CVector, - link::Tuple{NodeID, NodeID}, - )::Union{Missing, Float64} - a = get_state_index(getaxes(convergence), link[1]; inflow = false) - b = get_state_index(getaxes(convergence), link[2]) - return if isnothing(a) && isnothing(b) - missing - elseif isnothing(a) - convergence[b] - elseif isnothing(b) - convergence[a] + internal_flow_links = graph[].internal_flow_links + link_idx = get_link_index(link, internal_flow_links) + isnothing(link_idx) ? 0.0 : flow[link_idx] end end diff --git a/core/src/logging.jl b/core/src/logging.jl index cbe634e87..e1d503a49 100644 --- a/core/src/logging.jl +++ b/core/src/logging.jl @@ -61,32 +61,34 @@ end "Log the convergence bottlenecks." function log_bottlenecks(model; interrupt::Bool) (; cache, p, u) = model.integrator - (; p_independent) = p level = LoggingExtras.Warn # Indicate convergence bottlenecks if possible with the current algorithm return if hasproperty(cache, :nlsolver) - flow_error = if interrupt && p.p_independent.ncalls[1] > 0 - flow_error = p.p_independent.convergence ./ p.p_independent.ncalls[1] - else - temp_convergence = @. abs(cache.nlsolver.cache.atmp / u) - temp_convergence / finitemaximum(temp_convergence) - end + flow_error = @. abs(cache.nlsolver.cache.atmp / u) + flow_error ./= finitemaximum(flow_error) + (; state_ranges, basin, pid_control) = p.p_independent errors = Pair{Symbol, String}[] error_count = 0 max_errors = 5 # Iterate over the errors in descending order for i in sortperm(flow_error; rev = true) - node_id = Symbol(p_independent.node_id[i]) error = flow_error[i] - isnan(error) && continue # NaN are sorted as largest - # Stop reporting errors if they are too small or too many + isnan(error) && continue if error < 1 / length(flow_error) || error_count >= max_errors break end - push!(errors, node_id => @sprintf("%.2f", error * 100) * "%") + # Map state index to node ID + label = if i in state_ranges.basin + string(basin.node_id[i - first(state_ranges.basin) + 1]) + elseif i in state_ranges.integral + string(pid_control.node_id[i - first(state_ranges.integral) + 1]) + else + "state_$i" + end + push!(errors, Symbol(label) => @sprintf("%.2f", error * 100) * "%") error_count += 1 end if !isempty(errors) diff --git a/core/src/model.jl b/core/src/model.jl index 9e2764b4e..9eb805ec0 100644 --- a/core/src/model.jl +++ b/core/src/model.jl @@ -24,83 +24,6 @@ struct Model end end -""" -Get the Jacobian evaluation function via DifferentiationInterface.jl. -The time derivative is also supplied in case a Rosenbrock method is used. -""" -function get_diff_eval(du::CVector, u::CVector, p::Parameters, solver::Solver) - (; p_independent, state_and_time_dependent_cache, time_dependent_cache, p_mutable) = p - backend = get_ad_type(solver) - sparsity_detector = TracerSparsityDetector() - - backend_jac = if solver.sparse - AutoSparse(backend; sparsity_detector, coloring_algorithm = GreedyColoringAlgorithm()) - else - backend - end - - t = 0.0 - - jac_prep = prepare_jacobian( - water_balance!, - du, - backend_jac, - u, - Constant(p_independent), - Cache(state_and_time_dependent_cache), - Constant(time_dependent_cache), - Constant(p_mutable), - Constant(t); - strict = Val(true), - ) - - jac_prototype = solver.sparse ? sparsity_pattern(jac_prep) : nothing - - jac(J, u, p, t) = jacobian!( - water_balance!, - du, - J, - jac_prep, - backend_jac, - u, - Constant(p.p_independent), - Cache(state_and_time_dependent_cache), - Constant(time_dependent_cache), - Constant(p.p_mutable), - Constant(t), - ) - - tgrad_prep = prepare_derivative( - water_balance!, - du, - backend, - t, - Constant(u), - Constant(p_independent), - Cache(state_and_time_dependent_cache), - Cache(time_dependent_cache), - Constant(p_mutable); - strict = Val(true), - ) - tgrad(dT, u, p, t) = derivative!( - water_balance!, - du, - dT, - tgrad_prep, - backend, - t, - Constant(u), - Constant(p.p_independent), - Cache(state_and_time_dependent_cache), - Cache(time_dependent_cache), - Constant(p.p_mutable), - ) - - time_dependent_cache.t_prev_call[1] = -1.0 - - return jac_prototype, jac, tgrad -end - function Model(config_path::AbstractString)::Model return Model(Config(config_path)) end @@ -177,18 +100,18 @@ function Model(config::Config)::Model error("Model has no state.") end - reltol, relmask = build_reltol_vector(u0, config.solver.reltol) - parameters.p_independent.relmask .= relmask + # Initialize basin storages in the state vector + initialize_state_vector!(u0, parameters.p_independent) + du0 = zero(u0) # The Solver algorithm alg = algorithm(config.solver) # Synchronize level with storage - set_current_basin_properties!(p_independent.u_reduced, parameters, t0) + set_current_basin_properties!(u0, parameters, t0) - # Previous level is used to estimate the minimum level that was attained during a time step - # in limit_flow! + # Previous level is used to track level changes p_independent.basin.level_prev .= state_and_time_dependent_cache.current_level saveat = convert_saveat(config.solver.saveat, t_end) @@ -196,12 +119,16 @@ function Model(config::Config)::Model tstops = sort(unique(reduce(vcat, tstops))) adaptive = is_adaptive(config.solver.dt) - specialize = config.solver.specialize ? FullSpecialize : NoSpecialize - RHS = ODEFunction{true, specialize}( - water_balance!; - get_diff_eval(du0, parameters, config.solver)..., - ) - prob = ODEProblem{true, specialize}(RHS, u0, timespan, parameters) + prob = if is_mprk_algorithm(config.solver) + # Use PDSProblem for positivity-preserving MPRK methods + build_pds_problem(u0, timespan, parameters, p_independent) + else + RHS = ODEFunction{true, FullSpecialize}( + water_balance!; + get_diff_eval(du0, u0, parameters, config.solver)..., + ) + ODEProblem{true, FullSpecialize}(RHS, u0, timespan, parameters) + end @debug "Setup ODEProblem." callback, saved = create_callbacks(p_independent, config, saveat) @@ -212,29 +139,40 @@ function Model(config::Config)::Model # time which depends on the flows formulated in water_balance! water_balance!(du0, u0, parameters, t0) + # Initialize flow_rate_prev from the initial water_balance! so that + # the first trapezoidal integration step uses correct initial flow rates + sync_flow_rates!(parameters) + p_independent.flow_rate_prev .= p_independent.current_flow_rate + p_independent.evaporation_prev .= state_and_time_dependent_cache.current_evaporation + p_independent.infiltration_prev .= state_and_time_dependent_cache.current_infiltration + # Initialize the integrator, providing all solver options as described in # https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/ # Not all keyword arguments (e.g. `dtmax`) support `nothing`, in which case we follow # https://github.com/SciML/OrdinaryDiffEq.jl/blob/v6.57.0/src/solve.jl#L10 - integrator = init( - prob, - alg; - progress = true, - progress_name = "Simulating", - progress_steps = 100, - save_everystep = false, - callback, - tstops, - isoutofdomain, - adaptive, - config.solver.dt, - config.solver.dtmin, - dtmax = something(config.solver.dtmax, t_end), - config.solver.force_dtmin, - config.solver.abstol, - reltol, - config.solver.maxiters, + kwargs = Dict{Symbol, Any}( + :progress => true, + :progress_name => "Simulating", + :progress_steps => 100, + :save_everystep => false, + :callback => callback, + :tstops => tstops, + :adaptive => adaptive, + :dt => config.solver.dt, + :dtmin => config.solver.dtmin, + :dtmax => something(config.solver.dtmax, t_end), + :force_dtmin => config.solver.force_dtmin, + :abstol => config.solver.abstol, + :reltol => config.solver.reltol, + :maxiters => config.solver.maxiters, ) + + # MPRK methods guarantee positivity; skip isoutofdomain to avoid interference + if !is_mprk_algorithm(config.solver) + kwargs[:isoutofdomain] = isoutofdomain + end + + integrator = init(prob, alg; kwargs...) @debug "Setup integrator." model = Model(integrator, config, saved) diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 04523b3d3..26ec84c2e 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -8,28 +8,11 @@ const SolverStats = @NamedTuple{ dt::Float64, } -const state_components = ( - :tabulated_rating_curve, - :pump, - :outlet, - :user_demand_inflow, - :user_demand_outflow, - :linear_resistance, - :manning_resistance, - :evaporation, - :infiltration, - :integral, -) +const state_components = (:basin, :integral) const n_components = length(state_components) const StateTuple{V} = NamedTuple{state_components, NTuple{n_components, V}} const RibasimCVectorType{T} = Ribasim.CArrays.CArray{T, 1, Vector{T}, StateTuple{UnitRange{Int}}} -const RibasimReducedCVectorType{T} = Ribasim.CArrays.CArray{ - T, - 1, - Vector{T}, - @NamedTuple{combined_cumulative_flows::UnitRange{Int}, integral::UnitRange{Int}} -} # LinkType.flow and NodeType.FlowBoundary @enumx LinkType flow control listen observation none @@ -377,6 +360,80 @@ The parameter update associated with a certain control state for discrete contro allocation_controlled::Bool = false end +""" +Cache for mass-balance-consistent correction of trapezoidal flow estimates. +Projects trapezoidal estimates onto the subspace satisfying exact mass balance +using the minimum-norm adjustment: q̄ = q̄* + Aᵀ(AAᵀ+2I)⁻¹(b - Aq̄*) + +- `A_flow`: Signed incidence matrix (n_basins × n_links), +1 for inflow, -1 for outflow +- `AAt_2I_inv`: Precomputed (AAᵀ+2I)⁻¹, always positive definite +- `storage_prev`: Basin storage at previous accepted step +- `lambda`, `residual`, `correction_flow`: Preallocated work vectors +""" +struct BalanceCorrectionCache + A_flow::SparseMatrixCSC{Float64, Int} + AAt_2I_chol::Factorization{Float64} + nonnegative_link::BitVector + storage_prev::Vector{Float64} + lambda::Vector{Float64} + residual::Vector{Float64} + correction_flow::Vector{Float64} +end + +function BalanceCorrectionCache( + internal_flow_links::Vector{LinkMetadata}, + basin_node_ids::Vector{NodeID}, + n_links::Int, + ) + n_basins = length(basin_node_ids) + + # Build sparse signed incidence matrix: n_basins × n_links + I_idx = Int[] + J_idx = Int[] + V_val = Float64[] + for (j, link_metadata) in enumerate(internal_flow_links) + src, dst = link_metadata.link + if src.type == NodeType.Basin + push!(I_idx, src.idx) + push!(J_idx, j) + push!(V_val, -1.0) # outflow from source basin + end + if dst.type == NodeType.Basin + push!(I_idx, dst.idx) + push!(J_idx, j) + push!(V_val, +1.0) # inflow to destination basin + end + end + A_flow = sparse(I_idx, J_idx, V_val, n_basins, n_links) + + # (AAᵀ + 2I) is always positive definite; precompute sparse Cholesky factorization + AAt_sparse = A_flow * A_flow' + for i in 1:n_basins + AAt_sparse[i, i] += 2.0 + end + AAt_2I_chol = cholesky(Symmetric(AAt_sparse)) + + # Links connected to one-way controlled components should remain nonnegative + # in link orientation after correction. + nonnegative_link = BitVector(undef, n_links) + for (j, link_metadata) in enumerate(internal_flow_links) + src, dst = link_metadata.link + nonnegative_link[j] = + src.type in (NodeType.Pump, NodeType.Outlet, NodeType.TabulatedRatingCurve, NodeType.UserDemand) || + dst.type in (NodeType.Pump, NodeType.Outlet, NodeType.TabulatedRatingCurve, NodeType.UserDemand) + end + + return BalanceCorrectionCache( + A_flow, + AAt_2I_chol, + nonnegative_link, + zeros(n_basins), # storage_prev + zeros(n_basins), # lambda + zeros(n_basins), # residual + zeros(n_links), # correction_flow + ) +end + """ In-memory storage of saved mean flows for writing to results. @@ -393,6 +450,8 @@ In-memory storage of saved mean flows for writing to results. - `t`: Endtime of the interval over which is averaged """ @kwdef struct SavedFlow + + # Mean flow rates per internal flow link flow::Vector{Float64} inflow::Vector{Float64} outflow::Vector{Float64} @@ -400,12 +459,15 @@ In-memory storage of saved mean flows for writing to results. precipitation::Vector{Float64} surface_runoff::Vector{Float64} drainage::Vector{Float64} + # Mean evaporation and infiltration per basin + evaporation::Vector{Float64} + infiltration::Vector{Float64} concentration::Matrix{Float64} storage_rate::Vector{Float64} = zero(precipitation) balance_error::Vector{Float64} = zero(precipitation) relative_error::Vector{Float64} = zero(precipitation) - basin_convergence::Vector{Union{Missing, Float64}} - flow_convergence::Vector{Union{Missing, Float64}} + convergence::Vector{Union{Missing, Float64}} + flow_convergence::Vector{Float64} t::Float64 end @@ -787,12 +849,18 @@ const StateAndTimeDependentCache{T} = @NamedTuple{ current_area::Vector{T}, current_flow_rate_pump::Vector{T}, current_flow_rate_outlet::Vector{T}, + current_flow_rate_tabulated_rating_curve::Vector{T}, + current_flow_rate_linear_resistance::Vector{T}, + current_flow_rate_manning_resistance::Vector{T}, + current_flow_rate_user_demand::Vector{T}, + current_evaporation::Vector{T}, + current_infiltration::Vector{T}, current_error_pid_control::Vector{T}, - u_reduced_prev_call::Vector{T}, + u_prev_call::Vector{T}, t_prev_call::Vector{T}, } where {T} -@enumx CacheType flow_rate_pump flow_rate_outlet basin_level basin_storage +@enumx CacheType flow_rate_pump flow_rate_outlet flow_rate_tabulated_rating_curve flow_rate_linear_resistance flow_rate_manning_resistance flow_rate_user_demand basin_level basin_storage """ A cache for intermediate results in `water_balance!` which depend only on the time `t`. A second version of this @@ -839,7 +907,6 @@ of `water_balance!` (AD versus 'normal') which version of these objects is passe @kwdef struct CacheRef type::CacheType.T = CacheType.flow_rate_pump idx::Int = 0 - from_du::Bool = false end """ @@ -853,6 +920,14 @@ function get_cache_vector( state_and_time_dependent_cache.current_flow_rate_pump elseif type == CacheType.flow_rate_outlet state_and_time_dependent_cache.current_flow_rate_outlet + elseif type == CacheType.flow_rate_tabulated_rating_curve + state_and_time_dependent_cache.current_flow_rate_tabulated_rating_curve + elseif type == CacheType.flow_rate_linear_resistance + state_and_time_dependent_cache.current_flow_rate_linear_resistance + elseif type == CacheType.flow_rate_manning_resistance + state_and_time_dependent_cache.current_flow_rate_manning_resistance + elseif type == CacheType.flow_rate_user_demand + state_and_time_dependent_cache.current_flow_rate_user_demand elseif type == CacheType.basin_level state_and_time_dependent_cache.current_level elseif type == CacheType.basin_storage @@ -1087,6 +1162,7 @@ saveat: The time interval between saves of output data (storage, flow, ...) internal_flow_links: The metadata of the flow links used in the core without any Junctions. external_flow_links: The metadata of all flow links including those with Junctions. flow_link_map: A sparse matrix mapping internal_flow_ids to external_flow_ids. +flow_link_lookup: A dictionary mapping link tuples to their index in internal_flow_links for O(1) lookup. """ const ModelGraph = MetaGraph{ Int64, @@ -1100,6 +1176,7 @@ const ModelGraph = MetaGraph{ internal_flow_links::Vector{LinkMetadata}, external_flow_links::Vector{LinkMetadata}, flow_link_map::SparseMatrixCSC{Bool, Int}, + flow_link_lookup::Dict{Tuple{NodeID, NodeID}, Int}, }, Returns{Float64}, Float64, @@ -1127,8 +1204,6 @@ the object itself is not. """ @kwdef struct ParametersIndependent{C1} starttime::DateTime - reltol::Float64 - relmask::Vector{Bool} graph::ModelGraph allocation::Allocation basin::Basin @@ -1148,32 +1223,42 @@ the object itself is not. level_demand::LevelDemand flow_demand::FlowDemand subgrid::Subgrid - # Per state the in- and outflow links associated with that state (if they exist) - state_inflow_link::Vector{LinkMetadata} = LinkMetadata[] - state_outflow_link::Vector{LinkMetadata} = LinkMetadata[] - # Map each flow link to its state index. Used for link→state lookups where the - # destination node can have multiple inflow-link states (currently only UserDemand). - link_to_state_idx::Dict{Tuple{NodeID, NodeID}, Int} = - Dict{Tuple{NodeID, NodeID}, Int}() # Water balance tolerances water_balance_abstol::Float64 water_balance_reltol::Float64 - # State at previous saveat - u_prev_saveat::Vector{Float64} = Float64[] - # Node ID associated with each state - node_id::Vector{NodeID} = NodeID[] + # State ranges: basin and integral state_ranges::StateTuple{UnitRange{Int}} # Callback configurations do_concentration::Bool do_subgrid::Bool - temp_convergence::RibasimCVectorType{Float64} - convergence::RibasimCVectorType{Float64} ncalls::Vector{Int} = [0] - # Reduced state where the cumulative flows are combined into Basin - # storages (without non-state cumulative_flows) - u_reduced::RibasimReducedCVectorType{Float64} # Solver constants level_difference_threshold::Float64 + # Flow accumulation for trapezoidal callback (per internal flow link) + # current_flow_rate[i] holds the instantaneous flow for internal_flow_links[i] + current_flow_rate::Vector{Float64} = Float64[] + # Previous step's flow rates for trapezoidal integration + flow_rate_prev::Vector{Float64} = Float64[] + # Cumulative flow volumes per link (reset each timestep) + cumulative_flow::Vector{Float64} = Float64[] + # Cumulative flow volumes per link since previous saveat + cumulative_flow_saveat::Vector{Float64} = Float64[] + # Previous step evaporation/infiltration for trapezoidal integration + evaporation_prev::Vector{Float64} = Float64[] + infiltration_prev::Vector{Float64} = Float64[] + cumulative_evaporation::Vector{Float64} = Float64[] + cumulative_infiltration::Vector{Float64} = Float64[] + cumulative_evaporation_saveat::Vector{Float64} = Float64[] + cumulative_infiltration_saveat::Vector{Float64} = Float64[] + # Running totals for BMI (never reset) + cumulative_infiltration_total::Vector{Float64} = Float64[] + # Mass-balance-consistent correction cache + balance_correction::BalanceCorrectionCache + # Convergence tracking: accumulated normalized Newton residual per basin + convergence::Vector{Float64} = Float64[] + convergence_ncalls::Vector{Int} = [0] + # Accumulated absolute flow correction from balance correction (per internal flow link) + flow_convergence_saveat::Vector{Float64} = Float64[] end """ @@ -1185,7 +1270,12 @@ function StateAndTimeDependentCache( n_basin = length(p_independent.basin.node_id) n_pump = length(p_independent.pump.node_id) n_outlet = length(p_independent.outlet.node_id) + n_trc = length(p_independent.tabulated_rating_curve.node_id) + n_lr = length(p_independent.linear_resistance.node_id) + n_mr = length(p_independent.manning_resistance.node_id) + n_ud = length(p_independent.user_demand.node_id) n_pid_control = length(p_independent.pid_control.node_id) + n_states = n_basin + n_pid_control return (; current_storage = zeros(n_basin), @@ -1194,8 +1284,14 @@ function StateAndTimeDependentCache( current_area = zeros(n_basin), current_flow_rate_pump = zeros(n_pump), current_flow_rate_outlet = zeros(n_outlet), + current_flow_rate_tabulated_rating_curve = zeros(n_trc), + current_flow_rate_linear_resistance = zeros(n_lr), + current_flow_rate_manning_resistance = zeros(n_mr), + current_flow_rate_user_demand = zeros(n_ud), + current_evaporation = zeros(n_basin), + current_infiltration = zeros(n_basin), current_error_pid_control = zeros(n_pid_control), - u_reduced_prev_call = getdata(p_independent.u_reduced) .- 1.0, + u_prev_call = zeros(n_states) .- 1.0, t_prev_call = [-1.0], ) end @@ -1275,15 +1371,10 @@ end Base.show(io::IO, ::Parameters) = print(io, "Ribasim Parameters") Base.show(io::IO, ::MIME"text/plain", ::Parameters) = print(io, "Ribasim Parameters") -function get_value(ref::CacheRef, p::Parameters, du::CVector) - return if ref.from_du - du[ref.idx] - else - get_cache_vector(p.state_and_time_dependent_cache, ref.type)[ref.idx] - end +function get_value(ref::CacheRef, p::Parameters) + return get_cache_vector(p.state_and_time_dependent_cache, ref.type)[ref.idx] end function set_value!(ref::CacheRef, p::Parameters, value) - @assert !ref.from_du return get_cache_vector(p.state_and_time_dependent_cache, ref.type)[ref.idx] = value end diff --git a/core/src/pds.jl b/core/src/pds.jl new file mode 100644 index 000000000..562416099 --- /dev/null +++ b/core/src/pds.jl @@ -0,0 +1,577 @@ +""" +Production-Destruction System (PDS) formulation for PositiveIntegrators.jl. + +This module provides the PDS representation of the Ribasim water balance, +enabling the use of Modified Patankar-Runge-Kutta (MPRK) schemes that +unconditionally guarantee positivity of basin storages. + +The production matrix P[i,j] (off-diagonal) represents flow from basin j to basin i. +The diagonal P[i,i] represents external sources into basin i. +The destruction vector D[i] represents external sinks from basin i. +""" + +""" + build_pds_prototype(p_independent::ParametersIndependent) -> SparseMatrixCSC + +Build the sparsity pattern for the production matrix P from the network topology. +Returns a sparse matrix with structural nonzeros at positions where flows between +basins can occur. The diagonal has entries for every basin (external sources). +""" +function build_pds_prototype(p_independent::ParametersIndependent) + (; + basin, linear_resistance, manning_resistance, tabulated_rating_curve, + pump, outlet, user_demand, + ) = p_independent + + n_basins = length(basin.node_id) + + # Collect (row, col) pairs for structural nonzeros + I_idx = Int[] + J_idx = Int[] + + # Every basin has a diagonal entry for external sources (precipitation, drainage, etc.) + for i in 1:n_basins + push!(I_idx, i) + push!(J_idx, i) + end + + # Helper: add entry for a flow link between two nodes. + # If both are basins, add P[dst, src]. For bidirectional links, add both directions. + function add_link!(src_id::NodeID, dst_id::NodeID; bidirectional::Bool = false) + return if src_id.type == NodeType.Basin && dst_id.type == NodeType.Basin + push!(I_idx, dst_id.idx) + push!(J_idx, src_id.idx) + if bidirectional + push!(I_idx, src_id.idx) + push!(J_idx, dst_id.idx) + end + end + end + + # LinearResistance: bidirectional + for node_idx in eachindex(linear_resistance.node_id) + inflow_id = linear_resistance.inflow_link[node_idx].link[1] + outflow_id = linear_resistance.outflow_link[node_idx].link[2] + add_link!(inflow_id, outflow_id; bidirectional = true) + end + + # ManningResistance: bidirectional + for node_idx in eachindex(manning_resistance.node_id) + inflow_id = manning_resistance.inflow_link[node_idx].link[1] + outflow_id = manning_resistance.outflow_link[node_idx].link[2] + add_link!(inflow_id, outflow_id; bidirectional = true) + end + + # TabulatedRatingCurve: unidirectional + for node_idx in eachindex(tabulated_rating_curve.node_id) + inflow_id = tabulated_rating_curve.inflow_link[node_idx].link[1] + outflow_id = tabulated_rating_curve.outflow_link[node_idx].link[2] + add_link!(inflow_id, outflow_id; bidirectional = false) + end + + # Pump: unidirectional + for node_idx in eachindex(pump.node_id) + inflow_id = pump.inflow_link[node_idx].link[1] + outflow_id = pump.outflow_link[node_idx].link[2] + add_link!(inflow_id, outflow_id; bidirectional = false) + end + + # Outlet: unidirectional + for node_idx in eachindex(outlet.node_id) + inflow_id = outlet.inflow_link[node_idx].link[1] + outflow_id = outlet.outflow_link[node_idx].link[2] + add_link!(inflow_id, outflow_id; bidirectional = false) + end + + # Build sparse matrix; duplicate entries are summed by `sparse`, giving correct structure + V_val = ones(length(I_idx)) + p_prototype = sparse(I_idx, J_idx, V_val, n_basins, n_basins) + + return p_prototype +end + +""" + ribasim_P!(P::AbstractMatrix, u::CVector, p::Parameters, t) + +Fill the production matrix P for the PDS formulation. +- Off-diagonal P[i,j]: flow rate from basin j to basin i (must be ≥ 0) +- Diagonal P[i,i]: external sources into basin i (precipitation, drainage, flow boundaries) +""" +function ribasim_P!(P::AbstractMatrix, u::CVector, p::Parameters, t) + fill!(nonzeros(P), 0.0) + + # Ensure caches are up to date + check_new_input!(p, u, t) + set_current_basin_properties!(u, p, t) + + # External sources on the diagonal + fill_P_external_sources!(P, p, t) + + # Basin-to-basin flows (off-diagonal) and flows from non-basin sources + fill_P_linear_resistance!(P, p, t) + fill_P_manning_resistance!(P, p, t) + fill_P_tabulated_rating_curve!(P, p, t) + fill_P_pump_or_outlet!(P, p.p_independent.pump, p, t, false) + fill_P_pump_or_outlet!(P, p.p_independent.outlet, p, t, true) + fill_P_user_demand!(P, p, t) + + return nothing +end + +""" + ribasim_D!(D::AbstractVector, u::CVector, p::Parameters, t) + +Fill the destruction vector D for the PDS formulation. +D[i] contains the total rate of external sinks from basin i (evaporation, infiltration, +flows leaving the system through Terminal/LevelBoundary/UserDemand). +""" +function ribasim_D!(D::AbstractVector, u::CVector, p::Parameters, t) + fill!(D, 0.0) + + # Ensure caches are up to date (may already be current from P! call) + check_new_input!(p, u, t) + set_current_basin_properties!(u, p, t) + + (; p_independent, state_and_time_dependent_cache) = p + (; basin) = p_independent + (; vertical_flux) = basin + (; current_area, current_low_storage_factor) = state_and_time_dependent_cache + + # Evaporation and infiltration + for id in basin.node_id + i = id.idx + area = current_area[i] + factor = current_low_storage_factor[i] + evaporation = area * factor * vertical_flux.potential_evaporation[i] + infiltration = factor * vertical_flux.infiltration[i] + D[i] += evaporation + infiltration + end + + # Flows from basins to non-basin nodes (Terminal, LevelBoundary, etc.) + fill_D_outflows!(D, p, t) + + # UserDemand abstractions from source basins + fill_D_user_demand!(D, p, t) + + return nothing +end + +# --- Production matrix fill helpers --- + +function fill_P_external_sources!(P::AbstractMatrix, p::Parameters, t) + (; p_independent) = p + (; basin, flow_boundary) = p_independent + (; vertical_flux) = basin + + for id in basin.node_id + i = id.idx + fixed_area = basin_areas(basin, i)[end] + P[i, i] += fixed_area * vertical_flux.precipitation[i] + P[i, i] += vertical_flux.drainage[i] + P[i, i] += vertical_flux.surface_runoff[i] + end + + # Flow boundary contributions (positive = inflow to basin) + formulate_flow_boundary!(p, t) + for outflow_link in flow_boundary.outflow_link + from_id = outflow_link.link[1] + to_id = outflow_link.link[2] + if to_id.type == NodeType.Basin + q = flow_boundary.flow_rate[from_id.idx](t) + if q > 0 + P[to_id.idx, to_id.idx] += q + end + end + end + + return nothing +end + +function fill_P_linear_resistance!(P::AbstractMatrix, p::Parameters, t) + (; p_independent, state_and_time_dependent_cache) = p + (; linear_resistance) = p_independent + + for node_idx in eachindex(linear_resistance.node_id) + id = linear_resistance.node_id[node_idx] + inflow_id = linear_resistance.inflow_link[node_idx].link[1] + outflow_id = linear_resistance.outflow_link[node_idx].link[2] + + h_a = get_level(p, inflow_id, t) + h_b = get_level(p, outflow_id, t) + q = linear_resistance_flow(linear_resistance, id, h_a, h_b, p) + state_and_time_dependent_cache.current_flow_rate_linear_resistance[node_idx] = q + + _assign_flow_to_P!(P, q, inflow_id, outflow_id) + end + return nothing +end + +function fill_P_manning_resistance!(P::AbstractMatrix, p::Parameters, t) + (; p_independent, state_and_time_dependent_cache) = p + (; manning_resistance) = p_independent + + for node_idx in eachindex(manning_resistance.node_id) + id = manning_resistance.node_id[node_idx] + inflow_id = manning_resistance.inflow_link[node_idx].link[1] + outflow_id = manning_resistance.outflow_link[node_idx].link[2] + + h_a = get_level(p, inflow_id, t) + h_b = get_level(p, outflow_id, t) + q = manning_resistance_flow(manning_resistance, id, h_a, h_b, p) + state_and_time_dependent_cache.current_flow_rate_manning_resistance[node_idx] = q + + _assign_flow_to_P!(P, q, inflow_id, outflow_id) + end + return nothing +end + +function fill_P_tabulated_rating_curve!(P::AbstractMatrix, p::Parameters, t) + (; p_independent, state_and_time_dependent_cache) = p + (; tabulated_rating_curve) = p_independent + + for node_idx in eachindex(tabulated_rating_curve.node_id) + id = tabulated_rating_curve.node_id[node_idx] + inflow_id = tabulated_rating_curve.inflow_link[node_idx].link[1] + outflow_id = tabulated_rating_curve.outflow_link[node_idx].link[2] + + h_a = get_level(p, inflow_id, t) + h_b = get_level(p, outflow_id, t) + q_h = tabulated_rating_curve_flow(tabulated_rating_curve, id, h_a, h_b, p, t) + q = if tabulated_rating_curve.allocation_controlled[node_idx] + q_alloc = allocated_rating_curve_flow(tabulated_rating_curve, id, h_a, h_b, p) + min(q_alloc, q_h) + else + q_h + end + state_and_time_dependent_cache.current_flow_rate_tabulated_rating_curve[node_idx] = q + + _assign_flow_to_P!(P, q, inflow_id, outflow_id) + end + return nothing +end + +function fill_P_pump_or_outlet!( + P::AbstractMatrix, + node::Union{Pump, Outlet}, + p::Parameters, + t::Number, + reduce_Δlevel::Bool, + ) + (; p_independent, state_and_time_dependent_cache, time_dependent_cache) = p + (; allocation, flow_demand, level_difference_threshold) = p_independent + + current_flow_rate = node isa Pump ? + state_and_time_dependent_cache.current_flow_rate_pump : + state_and_time_dependent_cache.current_flow_rate_outlet + component_cache = node isa Pump ? time_dependent_cache.pump : time_dependent_cache.outlet + + for node_idx in eachindex(node.node_id) + # Phase 1: only handle non-controlled nodes + if node.control_type[node_idx] != ContinuousControlType.None + continue + end + + id = node.node_id[node_idx] + inflow_id = node.inflow_link[node_idx].link[1] + outflow_id = node.outflow_link[node_idx].link[2] + + q = _compute_pump_or_outlet_flow( + node, node_idx, p, t, + current_flow_rate, component_cache, reduce_Δlevel + ) + + _assign_flow_to_P!(P, q, inflow_id, outflow_id) + end + return nothing +end + +function fill_P_user_demand!(P::AbstractMatrix, p::Parameters, t) + (; p_independent, time_dependent_cache, state_and_time_dependent_cache) = p + (; user_demand, allocation, level_difference_threshold) = p_independent + (; current_return_factor) = time_dependent_cache.user_demand + + for node_idx in eachindex(user_demand.node_id) + id = user_demand.node_id[node_idx] + outflow_link = user_demand.outflow_link[node_idx] + outflow_id = outflow_link.link[2] + inflow_links = user_demand.inflow_links[node_idx] + has_demand_priority = view(user_demand.has_demand_priority, node_idx, :) + allocated = view(user_demand.allocated, node_idx, :) + return_factor = user_demand.return_factor[node_idx] + min_level = user_demand.min_level[node_idx] + link_alloc = user_demand.inflow_link_allocated[node_idx] + n_links = length(inflow_links) + + q_total_demand = 0.0 + for demand_priority_idx in eachindex(allocation.demand_priorities_all) + !has_demand_priority[demand_priority_idx] && continue + q_total_demand += min( + allocated[demand_priority_idx], + get_demand(user_demand, id, demand_priority_idx, t), + ) + end + equal_split = n_links == 0 ? 0.0 : q_total_demand / n_links + + q_total_actual = 0.0 + for (k, link_meta) in enumerate(inflow_links) + src_id = link_meta.link[1] + f_low_storage = get_low_storage_factor(p, src_id) + source_level = get_level(p, src_id, t) + f_reduction = reduction_factor(source_level - min_level, level_difference_threshold) + q_k_target = isinf(link_alloc[k]) ? equal_split : link_alloc[k] + q_total_actual += q_k_target * f_low_storage * f_reduction + end + + q_return = q_total_actual * + eval_time_interpolation(return_factor, current_return_factor, id.idx, p, t) + + state_and_time_dependent_cache.current_flow_rate_user_demand[id.idx] = q_total_actual + + # Return flow is an external source into the outflow basin + if outflow_id.type == NodeType.Basin && q_return > 0 + P[outflow_id.idx, outflow_id.idx] += q_return + end + end + return nothing +end + +# --- Destruction vector fill helpers --- + +function fill_D_outflows!(D::AbstractVector, p::Parameters, t) + (; p_independent, state_and_time_dependent_cache) = p + (; linear_resistance, manning_resistance, tabulated_rating_curve, pump, outlet, flow_boundary) = p_independent + + # LinearResistance: basin → non-basin + for node_idx in eachindex(linear_resistance.node_id) + inflow_id = linear_resistance.inflow_link[node_idx].link[1] + outflow_id = linear_resistance.outflow_link[node_idx].link[2] + q = state_and_time_dependent_cache.current_flow_rate_linear_resistance[node_idx] + _assign_flow_to_D!(D, q, inflow_id, outflow_id) + end + + # ManningResistance: basin → non-basin + for node_idx in eachindex(manning_resistance.node_id) + inflow_id = manning_resistance.inflow_link[node_idx].link[1] + outflow_id = manning_resistance.outflow_link[node_idx].link[2] + q = state_and_time_dependent_cache.current_flow_rate_manning_resistance[node_idx] + _assign_flow_to_D!(D, q, inflow_id, outflow_id) + end + + # TabulatedRatingCurve: basin → non-basin + for node_idx in eachindex(tabulated_rating_curve.node_id) + inflow_id = tabulated_rating_curve.inflow_link[node_idx].link[1] + outflow_id = tabulated_rating_curve.outflow_link[node_idx].link[2] + q = state_and_time_dependent_cache.current_flow_rate_tabulated_rating_curve[node_idx] + _assign_flow_to_D!(D, q, inflow_id, outflow_id) + end + + # Pump: basin → non-basin + for node_idx in eachindex(pump.node_id) + pump.control_type[node_idx] != ContinuousControlType.None && continue + inflow_id = pump.inflow_link[node_idx].link[1] + outflow_id = pump.outflow_link[node_idx].link[2] + q = state_and_time_dependent_cache.current_flow_rate_pump[pump.node_id[node_idx].idx] + _assign_flow_to_D!(D, q, inflow_id, outflow_id) + end + + # Outlet: basin → non-basin + for node_idx in eachindex(outlet.node_id) + outlet.control_type[node_idx] != ContinuousControlType.None && continue + inflow_id = outlet.inflow_link[node_idx].link[1] + outflow_id = outlet.outflow_link[node_idx].link[2] + q = state_and_time_dependent_cache.current_flow_rate_outlet[outlet.node_id[node_idx].idx] + _assign_flow_to_D!(D, q, inflow_id, outflow_id) + end + + # Negative flow boundary (outflow from basin) + for outflow_link in flow_boundary.outflow_link + from_id = outflow_link.link[1] + to_id = outflow_link.link[2] + if to_id.type == NodeType.Basin + q = flow_boundary.flow_rate[from_id.idx](t) + if q < 0 + D[to_id.idx] += -q + end + end + end + + return nothing +end + +function fill_D_user_demand!(D::AbstractVector, p::Parameters, t) + (; p_independent, state_and_time_dependent_cache) = p + (; user_demand, allocation, level_difference_threshold) = p_independent + + for node_idx in eachindex(user_demand.node_id) + id = user_demand.node_id[node_idx] + inflow_links = user_demand.inflow_links[node_idx] + has_demand_priority = view(user_demand.has_demand_priority, node_idx, :) + allocated = view(user_demand.allocated, node_idx, :) + min_level = user_demand.min_level[node_idx] + link_alloc = user_demand.inflow_link_allocated[node_idx] + n_links = length(inflow_links) + + q_total_demand = 0.0 + for demand_priority_idx in eachindex(allocation.demand_priorities_all) + !has_demand_priority[demand_priority_idx] && continue + q_total_demand += min( + allocated[demand_priority_idx], + get_demand(user_demand, id, demand_priority_idx, t), + ) + end + equal_split = n_links == 0 ? 0.0 : q_total_demand / n_links + + for (k, link_meta) in enumerate(inflow_links) + src_id = link_meta.link[1] + if src_id.type == NodeType.Basin + f_low_storage = get_low_storage_factor(p, src_id) + source_level = get_level(p, src_id, t) + f_reduction = reduction_factor(source_level - min_level, level_difference_threshold) + q_k_target = isinf(link_alloc[k]) ? equal_split : link_alloc[k] + q_k = q_k_target * f_low_storage * f_reduction + D[src_id.idx] += q_k + end + end + end + return nothing +end + +# --- Shared helpers --- + +""" +Assign a computed flow q to the production matrix P. +- If both nodes are basins: off-diagonal entry (sign-split for bidirectional) +- If only the destination is a basin: diagonal entry (external source) +- If only the source is a basin: skip (handled in D) +""" +function _assign_flow_to_P!(P::AbstractMatrix, q::Number, inflow_id::NodeID, outflow_id::NodeID) + if inflow_id.type == NodeType.Basin && outflow_id.type == NodeType.Basin + # Basin-to-basin: off-diagonal, sign-split + if q >= 0 + P[outflow_id.idx, inflow_id.idx] += q + else + P[inflow_id.idx, outflow_id.idx] += -q + end + elseif outflow_id.type == NodeType.Basin + # Non-basin (e.g. LevelBoundary) → Basin: external source on diagonal + if q >= 0 + P[outflow_id.idx, outflow_id.idx] += q + end + # If q < 0 it means flow FROM basin TO non-basin → handled in D + end + # If only inflow is basin → destruction, handled in D + return nothing +end + +""" +Assign a flow to the destruction vector D for flows that leave a basin to a non-basin node. +Only adds to D when a basin is losing water to the outside (non-basin destination). +""" +function _assign_flow_to_D!(D::AbstractVector, q::Number, inflow_id::NodeID, outflow_id::NodeID) + # Only handle basin → non-basin flows (off-system outflow) + if inflow_id.type == NodeType.Basin && outflow_id.type != NodeType.Basin + if q > 0 + D[inflow_id.idx] += q + end + elseif outflow_id.type == NodeType.Basin && inflow_id.type != NodeType.Basin + # Reverse flow: non-basin → basin, but q < 0 means basin is losing water + if q < 0 + D[outflow_id.idx] += -q + end + end + return nothing +end + +""" +Compute pump/outlet flow rate without modifying du. +""" +function _compute_pump_or_outlet_flow( + node::Union{Pump, Outlet}, + node_idx::Int, + p::Parameters, + t::Number, + current_flow_rate::Vector{<:Number}, + component_cache::NamedTuple, + reduce_Δlevel::Bool, + )::Float64 + (; p_independent) = p + (; allocation, flow_demand, level_difference_threshold) = p_independent + (; + current_min_flow_rate, current_max_flow_rate, + current_min_upstream_level, current_max_downstream_level, + ) = component_cache + + id = node.node_id[node_idx] + inflow_id = node.inflow_link[node_idx].link[1] + outflow_id = node.outflow_link[node_idx].link[2] + min_flow_rate_itp = node.min_flow_rate[node_idx] + max_flow_rate_itp = node.max_flow_rate[node_idx] + min_upstream_level = node.min_upstream_level[node_idx] + max_downstream_level = node.max_downstream_level[node_idx] + + flow_rate = if node.control_type[node_idx] != ContinuousControlType.None + current_flow_rate[id.idx] + elseif isassigned(node.time_dependent_flow_rate, node_idx) + eval_time_interpolation(node.time_dependent_flow_rate[node_idx], current_flow_rate, id.idx, p, t) + else + node.flow_rate[id.idx] + end + + src_level = get_level(p, inflow_id, t) + dst_level = get_level(p, outflow_id, t) + q = flow_rate * get_low_storage_factor(p, inflow_id) + + lower_bound = eval_time_interpolation(min_flow_rate_itp, current_min_flow_rate, node_idx, p, t) + upper_bound = eval_time_interpolation(max_flow_rate_itp, current_max_flow_rate, node_idx, p, t) + + if !is_active(allocation) + has_demand, flow_demand_id = has_external_demand(node, id) + if has_demand + total_demand = 0.0 + has_any_demand_priority = false + for (dpi, di) in enumerate(flow_demand.demand_interpolation[flow_demand_id.idx]) + if flow_demand.has_demand_priority[flow_demand_id.idx, dpi] + has_any_demand_priority = true + total_demand += di(t) + end + end + if has_any_demand_priority + lower_bound = clamp(total_demand, lower_bound, upper_bound) + end + end + end + q = clamp(q, lower_bound, upper_bound) + + if reduce_Δlevel + q *= reduction_factor(src_level - dst_level, level_difference_threshold) + end + + min_upstream_level_ = eval_time_interpolation(min_upstream_level, current_min_upstream_level, node_idx, p, t) + q *= reduction_factor(src_level - min_upstream_level_, level_difference_threshold) + + max_downstream_level_ = eval_time_interpolation(max_downstream_level, current_max_downstream_level, node_idx, p, t) + q *= reduction_factor(max_downstream_level_ - dst_level, level_difference_threshold) + + current_flow_rate[id.idx] = q + return q +end + +""" + build_pds_problem(u0, timespan, parameters, p_independent) -> PDSProblem + +Construct the PDSProblem for use with MPRK solvers. +The existing water_balance! is provided as std_rhs so that non-MPRK solvers can also be used. +""" +function build_pds_problem(u0::CVector, timespan, parameters::Parameters, p_independent::ParametersIndependent) + p_prototype = build_pds_prototype(p_independent) + + return PDSProblem( + ribasim_P!, + ribasim_D!, + u0, + timespan, + parameters; + p_prototype, + std_rhs = water_balance!, + ) +end diff --git a/core/src/read.jl b/core/src/read.jl index 4b8f54cbd..aa29b4f2c 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -1648,69 +1648,59 @@ function Parameters(db::DB, config::Config)::Parameters u_ids = state_node_ids( (; - nodes.tabulated_rating_curve, - nodes.pump, - nodes.outlet, - nodes.user_demand, - nodes.linear_resistance, - nodes.manning_resistance, nodes.basin, nodes.pid_control, ) ) - node_id = reduce(vcat, u_ids) - n_states = length(node_id) state_ranges = count_state_ranges(u_ids) - state_inflow_link, state_outflow_link = get_state_flow_links(graph, nodes) - link_to_state_idx = build_link_to_state_idx(state_inflow_link) set_target_ref!( nodes.pid_control.target_ref, nodes.pid_control.node_id, - fill("flow_rate", length(node_id)), - state_ranges, + fill("flow_rate", length(nodes.pid_control.node_id)), graph, ) set_target_ref!( nodes.continuous_control.target_ref, nodes.continuous_control.node_id, nodes.continuous_control.controlled_variable, - state_ranges, graph, ) n_basin = length(nodes.basin.node_id) - n_pid_control = length(nodes.pid_control.node_id) - u_reduced = CVector( - zeros(n_basin + n_pid_control), - (; - combined_cumulative_flows = 1:n_basin, - integral = (n_basin + 1):(n_basin + n_pid_control), - ), + n_internal_flow_links = length(graph[].internal_flow_links) + + balance_correction = BalanceCorrectionCache( + graph[].internal_flow_links, nodes.basin.node_id, n_internal_flow_links, ) p_independent = ParametersIndependent(; config.starttime, - config.solver.reltol, - relmask = collect(trues(n_states)), graph, allocation, nodes..., subgrid, - state_inflow_link, - state_outflow_link, - link_to_state_idx, config.solver.water_balance_abstol, config.solver.water_balance_reltol, - u_prev_saveat = zeros(n_states), - node_id, state_ranges, do_concentration = config.experimental.concentration, do_subgrid = config.results.subgrid, - temp_convergence = CVector(zeros(n_states), state_ranges), - convergence = CVector(zeros(n_states), state_ranges), - u_reduced, config.solver.level_difference_threshold, + # Flow accumulation vectors + current_flow_rate = zeros(n_internal_flow_links), + flow_rate_prev = zeros(n_internal_flow_links), + cumulative_flow = zeros(n_internal_flow_links), + cumulative_flow_saveat = zeros(n_internal_flow_links), + evaporation_prev = zeros(n_basin), + infiltration_prev = zeros(n_basin), + cumulative_evaporation = zeros(n_basin), + cumulative_infiltration = zeros(n_basin), + cumulative_evaporation_saveat = zeros(n_basin), + cumulative_infiltration_saveat = zeros(n_basin), + cumulative_infiltration_total = zeros(n_basin), + balance_correction, + convergence = zeros(n_basin), + flow_convergence_saveat = zeros(n_internal_flow_links), ) collect_control_mappings!(p_independent) diff --git a/core/src/solve.jl b/core/src/solve.jl index d984d5f1b..e89f83711 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -1,5 +1,7 @@ """ The right hand side function of the system of ODEs set up by Ribasim. +State vector u contains basin storages and PID integral terms. +du contains dS/dt per basin and d(integral)/dt per PID control. """ water_balance!(du::CVector, u::CVector, p::Parameters, t::Number)::Nothing = water_balance!( du::RibasimCVectorType, @@ -30,11 +32,9 @@ water_balance!( t, ) -# Method where u is already parsed to u_reduced so this part is skipped in AD Jacobian computation - function water_balance!( du::RibasimCVectorType, - u_reduced::RibasimReducedCVectorType, + u::RibasimCVectorType, p_independent::ParametersIndependent, state_and_time_dependent_cache::StateAndTimeDependentCache, time_dependent_cache::TimeDependentCache, @@ -49,12 +49,12 @@ function water_balance!( ) # Check whether t or u is different from the last water_balance! call - check_new_input!(p, u_reduced, t) + check_new_input!(p, u, t) du .= 0.0 - # Ensures current_* vectors are current - set_current_basin_properties!(u_reduced, p, t) + # Ensures current_* vectors are current (storage, level, area from u.basin) + set_current_basin_properties!(u, p, t) # Notes on the ordering of these formulations: # - Continuous control can depend on flows (which are not continuously controlled themselves), @@ -62,8 +62,8 @@ function water_balance!( # - Pid control can depend on the du of basins and subsequently change them # because of the error derivative term. - # Basin forcings - update_vertical_flux!(du, p) + # Basin forcings (precipitation, evaporation, infiltration, drainage, surface_runoff) + update_vertical_flux!(du, p, t) # Formulate intermediate flows (non continuously controlled) formulate_flows!(du, p, t) @@ -75,7 +75,7 @@ function water_balance!( formulate_flows!(du, p, t; control_type = ContinuousControlType.Continuous) # Compute PID control - formulate_pid_control!(du, u_reduced, p, t) + formulate_pid_control!(du, u, p, t) # Formulate intermediate flow (controlled by PID control) formulate_flows!(du, p, t; control_type = ContinuousControlType.PID) @@ -105,7 +105,7 @@ function formulate_continuous_control!(du::CVector, p::Parameters, t::Number)::N cvar = compound_variable[i] ref = target_ref[i] func_ = func[i] - value = compound_variable_value(cvar, p, du, t) + value = compound_variable_value(cvar, p, t) set_value!(ref, p, func_(value)) end @@ -113,12 +113,10 @@ function formulate_continuous_control!(du::CVector, p::Parameters, t::Number)::N end """ -Compute the storages, levels and areas of all Basins given the -state u and the time t. +Compute the levels and areas of all Basins given the state u (which contains storage directly). """ - function set_current_basin_properties!( - u_reduced::RibasimReducedCVectorType, + u::RibasimCVectorType, p::Parameters, t::Number, )::Nothing @@ -126,18 +124,14 @@ function set_current_basin_properties!( (; basin) = p_independent (; - node_id, - cumulative_precipitation, - cumulative_surface_runoff, - cumulative_drainage, - vertical_flux, - low_storage_threshold, + low_storage_threshold, vertical_flux, + cumulative_precipitation, cumulative_surface_runoff, cumulative_drainage, ) = basin - # The exact cumulative precipitation and drainage up to the t of this water_balance call + # The exact cumulative precipitation, drainage, surface_runoff up to the t of this water_balance! call if p_mutable.new_time_dependent_cache dt = t - p_mutable.tprev - for id in node_id + for id in basin.node_id fixed_area = basin_areas(basin, id.idx)[end] time_dependent_cache.basin.current_cumulative_precipitation[id.idx] = cumulative_precipitation[id.idx] + @@ -150,11 +144,11 @@ function set_current_basin_properties!( end return if p_mutable.new_state_and_time_dependent_cache - formulate_storages!(u_reduced, p, t) + # Storage is directly in u.basin + state_and_time_dependent_cache.current_storage .= u.basin for i in eachindex(basin.node_id) - id = basin.node_id[i] - s = state_and_time_dependent_cache.current_storage[i] - i = id.idx + # Force storage positive, since Rosenbrock methods can overshoot to negative storages + s = ifelse(u.basin[i] > 0.0, u.basin[i], 0.0) state_and_time_dependent_cache.current_low_storage_factor[i] = reduction_factor(s, low_storage_threshold[i]) @inbounds state_and_time_dependent_cache.current_level[i] = @@ -165,62 +159,45 @@ function set_current_basin_properties!( end end -function formulate_storages!( - u_reduced::RibasimReducedCVectorType, - p::Parameters, - t::Number; - add_initial_storage::Bool = true, - )::Nothing - (; p_independent, state_and_time_dependent_cache, time_dependent_cache, p_mutable) = p - (; basin, flow_boundary) = p_independent - (; current_storage) = state_and_time_dependent_cache - - # Current storage: initial condition + - # total inflows and outflows since the start - # of the simulation - if add_initial_storage - current_storage .= basin.storage0 - else - current_storage .= 0.0 - end - - current_storage .+= u_reduced.combined_cumulative_flows - current_storage .+= time_dependent_cache.basin.current_cumulative_precipitation - current_storage .+= time_dependent_cache.basin.current_cumulative_surface_runoff - current_storage .+= time_dependent_cache.basin.current_cumulative_drainage - - # Formulate storage contributions of flow boundaries - formulate_flow_boundary!(p, t) - for (outflow_link, cumulative_flow) in zip( - flow_boundary.outflow_link, - time_dependent_cache.flow_boundary.current_cumulative_boundary_flow, - ) - outflow_id = outflow_link.link[2] - if outflow_id.type == NodeType.Basin - current_storage[outflow_id.idx] += cumulative_flow - end - end - return nothing -end - """ -Smoothly let the evaporation and infiltration flux go to 0 when the storage is less than 10 m^3 +Add vertical flux contributions (precipitation, evaporation, infiltration, drainage, surface_runoff) +directly to basin storage derivatives. """ -function update_vertical_flux!(du::CVector, p::Parameters)::Nothing +function update_vertical_flux!(du::CVector, p::Parameters, t::Number)::Nothing (; p_independent, state_and_time_dependent_cache) = p (; basin) = p_independent (; vertical_flux) = basin (; current_area, current_low_storage_factor) = state_and_time_dependent_cache for id in basin.node_id - area = current_area[id.idx] - factor = current_low_storage_factor[id.idx] - - evaporation = area * factor * vertical_flux.potential_evaporation[id.idx] - infiltration = factor * vertical_flux.infiltration[id.idx] + i = id.idx + area = current_area[i] + factor = current_low_storage_factor[i] + fixed_area = basin_areas(basin, i)[end] + + evaporation = area * factor * vertical_flux.potential_evaporation[i] + infiltration = factor * vertical_flux.infiltration[i] + + # Store current evaporation/infiltration rates in cache for flow output + state_and_time_dependent_cache.current_evaporation[i] = evaporation + state_and_time_dependent_cache.current_infiltration[i] = infiltration + + # Add vertical fluxes to basin storage derivative + du.basin[i] += fixed_area * vertical_flux.precipitation[i] + du.basin[i] += vertical_flux.surface_runoff[i] + du.basin[i] += vertical_flux.drainage[i] + du.basin[i] -= evaporation + du.basin[i] -= infiltration + end - du.evaporation[id.idx] = evaporation - du.infiltration[id.idx] = infiltration + # Flow boundary contributions + formulate_flow_boundary!(p, t) + for outflow_link in p_independent.flow_boundary.outflow_link + from_id = outflow_link.link[1] + to_id = outflow_link.link[2] + if to_id.type == NodeType.Basin + du.basin[to_id.idx] += p_independent.flow_boundary.flow_rate[from_id.idx](t) + end end return nothing @@ -245,7 +222,7 @@ end function formulate_pid_control!( du::CVector, - u_reduced::CVector, + u::CVector, p::Parameters, t::Number, )::Nothing @@ -291,7 +268,7 @@ function formulate_pid_control!( end if !iszero(K_i) - flow_rate += K_i * u_reduced.integral[i] / D + flow_rate += K_i * u.integral[i] / D end if !iszero(K_d) @@ -301,8 +278,8 @@ function formulate_pid_control!( else dtarget = derivative(target[i], t) end - dstorage_listened_basin_old = - formulate_dstorage_wrt_time(du, p_independent, t, listened_node_id) + # With basin storage states, du.basin[idx] already contains the current dstorage + dstorage_listened_basin_old = du.basin[listened_node_id.idx] # The expression below is the solution to an implicit equation for # dstorage_listened_basin. This equation results from the fact that if the derivative # term in the PID controller is used, the controlled pump flow rate depends on itself. @@ -316,33 +293,22 @@ function formulate_pid_control!( end """ -Formulate the time derivative of the storage in a single Basin. +Apply a flow rate q from inflow_id → outflow_id to basin storage derivatives. +Positive q flows from inflow_id to outflow_id. """ -function formulate_dstorage_wrt_time( +function apply_flow_to_basins!( du::CVector, - p_independent::ParametersIndependent, - t::Number, - node_id::NodeID, - ) - (; basin) = p_independent - (; inflow_ids, outflow_ids, vertical_flux) = basin - @assert node_id.type == NodeType.Basin - dstorage = 0.0 - for inflow_id in inflow_ids[node_id.idx] - dstorage += get_flow(du, p_independent, t, (inflow_id, node_id)) + q::Number, + inflow_id::NodeID, + outflow_id::NodeID, + )::Nothing + if inflow_id.type == NodeType.Basin + du.basin[inflow_id.idx] -= q end - for outflow_id in outflow_ids[node_id.idx] - dstorage -= get_flow(du, p_independent, t, (node_id, outflow_id)) + if outflow_id.type == NodeType.Basin + du.basin[outflow_id.idx] += q end - - fixed_area = basin_areas(basin, node_id.idx)[end] - dstorage += fixed_area * vertical_flux.precipitation[node_id.idx] - dstorage += vertical_flux.surface_runoff[node_id.idx] - dstorage += vertical_flux.drainage[node_id.idx] - dstorage -= du.evaporation[node_id.idx] - dstorage -= du.infiltration[node_id.idx] - - return dstorage + return nothing end function formulate_flow!( @@ -351,7 +317,7 @@ function formulate_flow!( p::Parameters, t::Number, )::Nothing - (; p_independent, time_dependent_cache) = p + (; p_independent, time_dependent_cache, state_and_time_dependent_cache) = p (; current_return_factor) = time_dependent_cache.user_demand (; allocation, level_difference_threshold) = p_independent @@ -392,13 +358,22 @@ function formulate_flow!( ) q_k_target = isinf(link_alloc[k]) ? equal_split : link_alloc[k] q_k = q_k_target * f_low_storage * f_reduction - du.user_demand_inflow[link_offset + k] = q_k + # Apply each inflow link's abstraction to the source basin + apply_flow_to_basins!(du, q_k, src_id, id) q_total_actual += q_k end - du.user_demand_outflow[id.idx] = + q_return = q_total_actual * eval_time_interpolation(return_factor, current_return_factor, id.idx, p, t) + + # Store total flow rate in cache for control listening + state_and_time_dependent_cache.current_flow_rate_user_demand[id.idx] = q_total_actual + + # Apply return flow (from UserDemand to downstream basin) + outflow_link = user_demand.outflow_link[node_idx] + outflow_id = outflow_link.link[2] + apply_flow_to_basins!(du, q_return, id, outflow_id) end return nothing end @@ -409,7 +384,7 @@ function formulate_flow!( p::Parameters, t::Number, )::Nothing - (; p_mutable) = p + (; state_and_time_dependent_cache) = p (; node_id) = linear_resistance for node_idx in eachindex(linear_resistance.node_id) @@ -423,7 +398,8 @@ function formulate_flow!( h_a = get_level(p, inflow_id, t) h_b = get_level(p, outflow_id, t) q = linear_resistance_flow(linear_resistance, id, h_a, h_b, p) - du.linear_resistance[node_idx] = q + state_and_time_dependent_cache.current_flow_rate_linear_resistance[node_idx] = q + apply_flow_to_basins!(du, q, inflow_id, outflow_id) end return nothing end @@ -500,7 +476,7 @@ function formulate_flow!( p::Parameters, t::Number, )::Nothing - (; p_mutable) = p + (; state_and_time_dependent_cache) = p for node_idx in eachindex(tabulated_rating_curve.node_id) id = tabulated_rating_curve.node_id[node_idx] inflow_link = tabulated_rating_curve.inflow_link[node_idx] @@ -512,15 +488,14 @@ function formulate_flow!( q_h = tabulated_rating_curve_flow(tabulated_rating_curve, id, h_a, h_b, p, t) q = if tabulated_rating_curve.allocation_controlled[node_idx] - # Ensure q is always >= to the Q(h) relationship, since errors in the linear approximations in allocation could lead to - # a higher q at the current h than the user defined q(h) would allow q_alloc = allocated_rating_curve_flow(tabulated_rating_curve, id, h_a, h_b, p) min(q_alloc, q_h) else q_h end - du.tabulated_rating_curve[node_idx] = q + state_and_time_dependent_cache.current_flow_rate_tabulated_rating_curve[node_idx] = q + apply_flow_to_basins!(du, q, inflow_id, outflow_id) end return nothing end @@ -633,7 +608,7 @@ function formulate_flow!( p::Parameters, t::Number, )::Nothing - (; p_mutable) = p + (; state_and_time_dependent_cache) = p (; node_id) = manning_resistance for node_idx in eachindex(manning_resistance.node_id) @@ -649,13 +624,14 @@ function formulate_flow!( q = manning_resistance_flow(manning_resistance, id, h_a, h_b, p) - du.manning_resistance[node_idx] = q + state_and_time_dependent_cache.current_flow_rate_manning_resistance[node_idx] = q + apply_flow_to_basins!(du, q, inflow_id, outflow_id) end return nothing end function formulate_pump_or_outlet_flow!( - du_component::SubArray{<:Number}, + du::CVector, node::Union{Pump, Outlet}, p::Parameters, t::Number, @@ -689,7 +665,6 @@ function formulate_pump_or_outlet_flow!( flow_rate = if control_type != ContinuousControlType.None current_flow_rate[id.idx] elseif isassigned(node.time_dependent_flow_rate, node_idx) - # get the time dependent flow rate from interpolation or cached value eval_time_interpolation( node.time_dependent_flow_rate[node_idx], current_flow_rate, @@ -698,7 +673,6 @@ function formulate_pump_or_outlet_flow!( t, ) else - # get the scalar flow rate from (for DiscreteControl, Control by allocation or flows from the Static table) node.flow_rate[id.idx] end @@ -764,7 +738,10 @@ function formulate_pump_or_outlet_flow!( ) q *= reduction_factor(max_downstream_level_ - dst_level, level_difference_threshold) - du_component[node_idx] = q + # Store in cache for control listening + current_flow_rate[id.idx] = q + # Apply to basin derivatives + apply_flow_to_basins!(du, q, inflow_id, outflow_id) end return nothing end @@ -778,7 +755,7 @@ function formulate_flow!( )::Nothing (; time_dependent_cache, state_and_time_dependent_cache) = p return formulate_pump_or_outlet_flow!( - du.pump, + du, pump, p, t, @@ -797,7 +774,7 @@ function formulate_flow!( )::Nothing (; time_dependent_cache, state_and_time_dependent_cache) = p return formulate_pump_or_outlet_flow!( - du.outlet, + du, outlet, p, t, @@ -834,168 +811,78 @@ function formulate_flows!( end """ -Clamp the cumulative flow states within the minimum and maximum -flow rates for the last time step if these flow rate bounds are known. +Get the Jacobian evaluation function via DifferentiationInterface.jl. +The time derivative is also supplied in case a Rosenbrock method is used. """ -function limit_flow!( - u::CVector, - integrator::DEIntegrator, - p::Parameters, - t::Number, - )::Nothing - (; uprev, dt) = integrator - (; p_independent, state_and_time_dependent_cache) = p - (; - pump, - outlet, - linear_resistance, - user_demand, - tabulated_rating_curve, - basin, - allocation, - u_reduced, - level_difference_threshold, - ) = p_independent - (; current_storage, current_level) = state_and_time_dependent_cache - - # The current storage and level based on the proposed u are used to estimate the lowest - # storage and level attained in the last time step to estimate whether there was an effect - # of reduction factors - - reduce_state!(u_reduced, u, p_independent) - set_current_basin_properties!(u_reduced, p, t) - - # TabulatedRatingCurve flow is in [0, ∞) - for id in tabulated_rating_curve.node_id - limit_flow!( - u.tabulated_rating_curve, - uprev.tabulated_rating_curve, - id, - 0.0, - Inf, - dt, - ) - end - - # Pump flow is in [min_flow_rate, max_flow_rate] - for (id, min_flow_rate, max_flow_rate) in - zip(pump.node_id, pump.min_flow_rate, pump.max_flow_rate) - limit_flow!(u.pump, uprev.pump, id, min_flow_rate(t), max_flow_rate(t), dt) - end - - # Outlet flow is in [min_flow_rate, max_flow_rate] - for (id, min_flow_rate, max_flow_rate) in - zip(outlet.node_id, outlet.min_flow_rate, outlet.max_flow_rate) - limit_flow!( - u.outlet, - uprev.outlet, - id, - min_flow_rate(t), - max_flow_rate(t), - dt, - ) - end +function get_diff_eval(du::CVector, u::CVector, p::Parameters, solver::Solver) + (; p_independent, state_and_time_dependent_cache, time_dependent_cache, p_mutable) = p + backend = get_ad_type(solver) + sparsity_detector = TracerSparsityDetector() - # LinearResistance flow is in [-max_flow_rate, max_flow_rate] - for (id, max_flow_rate) in zip( - linear_resistance.node_id, - linear_resistance.max_flow_rate, - ) - limit_flow!( - u.linear_resistance, - uprev.linear_resistance, - id, - -max_flow_rate, - max_flow_rate, - dt, - ) + backend_jac = if solver.sparse + AutoSparse(backend; sparsity_detector, coloring_algorithm = GreedyColoringAlgorithm()) + else + backend end - # UserDemand per inflow link bounds - for node_idx in eachindex(user_demand.node_id) - id = user_demand.node_id[node_idx] - inflow_links = user_demand.inflow_links[node_idx] - link_offset = user_demand.inflow_link_offsets[node_idx] - n_links = length(inflow_links) - demand_from_timeseries = user_demand.demand_from_timeseries[node_idx] - link_alloc = user_demand.inflow_link_allocated[node_idx] - - allocated_total = if demand_from_timeseries - 0.0 - else - sum( - min( - user_demand.demand[id.idx, demand_priority_idx], - user_demand.allocated[id.idx, demand_priority_idx], - ) for demand_priority_idx in eachindex(allocation.demand_priorities_all) - ) - end - equal_split = n_links == 0 ? 0.0 : allocated_total / n_links + t = 0.0 + + jac_prep = prepare_jacobian( + water_balance!, + du, + backend_jac, + u, + Constant(p_independent), + Cache(state_and_time_dependent_cache), + Constant(time_dependent_cache), + Constant(p_mutable), + Constant(t); + strict = Val(true), + ) - for (k, link_meta) in enumerate(inflow_links) - state_idx = link_offset + k - q_k_max = isinf(link_alloc[k]) ? equal_split : link_alloc[k] - min_flow_rate, max_flow_rate = if demand_from_timeseries - 0.0, Inf - else - src_id = link_meta.link[1] - factor_basin_min = min_low_storage_factor( - current_storage, - basin.storage_prev, - basin, - src_id, - ) - factor_level_min = min_low_user_demand_level_factor( - current_level, - basin.level_prev, - user_demand.min_level, - id, - src_id, - level_difference_threshold, - ) - factor_basin_min * factor_level_min * q_k_max, q_k_max - end - u_prev = uprev.user_demand_inflow[state_idx] - u.user_demand_inflow[state_idx] = clamp( - u.user_demand_inflow[state_idx], - u_prev + min_flow_rate * dt, - u_prev + max_flow_rate * dt, - ) - end - end + jac_prototype = solver.sparse ? Float64.(sparsity_pattern(jac_prep)) : nothing + + jac(J, u, p, t) = jacobian!( + water_balance!, + du, + J, + jac_prep, + backend_jac, + u, + Constant(p.p_independent), + Cache(state_and_time_dependent_cache), + Constant(time_dependent_cache), + Constant(p.p_mutable), + Constant(t), + ) - # Evaporation is in [0, ∞) (stricter bounds would require also estimating the area) - # Infiltration is in [f * infiltration, infiltration] where f is a rough estimate of the smallest low storage factor - # reduction factor value that was attained over the last timestep - for (id, infiltration) in zip(basin.node_id, basin.vertical_flux.infiltration) - factor_min = min_low_storage_factor(current_storage, basin.storage_prev, basin, id) - limit_flow!(u.evaporation, uprev.evaporation, id, 0.0, Inf, dt) - limit_flow!( - u.infiltration, - uprev.infiltration, - id, - factor_min * infiltration, - infiltration, - dt, - ) - end + tgrad_prep = prepare_derivative( + water_balance!, + du, + backend, + t, + Constant(u), + Constant(p_independent), + Cache(state_and_time_dependent_cache), + Cache(time_dependent_cache), + Constant(p_mutable); + strict = Val(true), + ) + tgrad(dT, u, p, t) = derivative!( + water_balance!, + du, + dT, + tgrad_prep, + backend, + t, + Constant(u), + Constant(p.p_independent), + Cache(state_and_time_dependent_cache), + Cache(time_dependent_cache), + Constant(p.p_mutable), + ) - return nothing -end + time_dependent_cache.t_prev_call[1] = -1.0 -function limit_flow!( - u_component, - uprev_component, - id::NodeID, - min_flow_rate::Number, - max_flow_rate::Number, - dt::Number, - )::Nothing - u_prev = uprev_component[id.idx] - u_component[id.idx] = clamp( - u_component[id.idx], - u_prev + min_flow_rate * dt, - u_prev + max_flow_rate * dt, - ) - return nothing + return (; jac_prototype, jac, tgrad) end diff --git a/core/src/util.jl b/core/src/util.jl index b1d8ef4e3..44e5815dc 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -40,7 +40,9 @@ end Compute the level of a basin given its storage. """ function get_level_from_storage(basin::Basin, state_idx::Int, storage::T)::T where {T} - return basin.storage_to_level[state_idx](storage) + # Clamp storage positive, since Rosenbrock methods can overshoot to negative + s = ifelse(storage > zero(T), storage, zero(T)) + return basin.storage_to_level[state_idx](s) end """ @@ -485,13 +487,31 @@ function NodeID(type::Symbol, value::Integer, p_independent::ParametersIndepende return NodeID(node_type, value, idx) end +"Map a NodeType to its corresponding CacheType for flow rate listening." +function cache_type_for_node(node_type::NodeType.T)::CacheType.T + return if node_type == NodeType.Pump + CacheType.flow_rate_pump + elseif node_type == NodeType.Outlet + CacheType.flow_rate_outlet + elseif node_type == NodeType.TabulatedRatingCurve + CacheType.flow_rate_tabulated_rating_curve + elseif node_type == NodeType.LinearResistance + CacheType.flow_rate_linear_resistance + elseif node_type == NodeType.ManningResistance + CacheType.flow_rate_manning_resistance + elseif node_type == NodeType.UserDemand + CacheType.flow_rate_user_demand + else + error("No cache type for node type $node_type") + end +end + """ Get the reference to a parameter """ function get_cache_ref( node_id::NodeID, - variable::String, - state_ranges::StateTuple{UnitRange{Int}}; + variable::String; listen::Bool = true, )::Tuple{CacheRef, Bool} errors = false @@ -507,9 +527,8 @@ function get_cache_ref( @error "Cannot listen to flow_rate of $node_id, the node type must be one of $conservative_nodetypes." CacheRef() else - # Index in the state vector (inflow) - idx = get_state_index(state_ranges, node_id) - CacheRef(; idx, from_du = true) + type = cache_type_for_node(node_id.type) + CacheRef(; type, node_id.idx) end else type = if node_id.type == NodeType.Pump @@ -534,7 +553,7 @@ end Set references to all variables that are listened to by discrete/continuous control """ function set_listen_cache_refs!(p_independent::ParametersIndependent)::Nothing - (; discrete_control, continuous_control, state_ranges) = p_independent + (; discrete_control, continuous_control) = p_independent compound_variable_sets = [discrete_control.compound_variables..., continuous_control.compound_variable] errors = false @@ -546,7 +565,6 @@ function set_listen_cache_refs!(p_independent::ParametersIndependent)::Nothing ref, error = get_cache_ref( subvariable.listen_node_id, subvariable.variable, - state_ranges, ) if !error subvariables[j] = @set subvariable.cache_ref = ref @@ -627,14 +645,13 @@ function set_target_ref!( target_ref::Vector{CacheRef}, node_id::Vector{NodeID}, controlled_variable::Vector{String}, - state_ranges::StateTuple{UnitRange{Int}}, graph::MetaGraph, )::Nothing errors = false for (i, (id, variable)) in enumerate(zip(node_id, controlled_variable)) controlled_node_id = only(outneighbor_labels_type(graph, id, LinkType.control)) ref, error = - get_cache_ref(controlled_node_id, variable, state_ranges; listen = false) + get_cache_ref(controlled_node_id, variable; listen = false) target_ref[i] = ref errors |= error end @@ -690,237 +707,189 @@ low_storage_factor_resistance_node(::Parameters, q::GradientTracer, ::NodeID, :: relaxed_root(x::GradientTracer, threshold::Real) = x get_level_from_storage(basin::Basin, state_idx::Int, storage::GradientTracer) = storage -"Create a NamedTuple of the node IDs per state component in the state order" -function state_node_ids( - p::Union{ParametersIndependent, NamedTuple}, - )::StateTuple{Vector{NodeID}} - return (; - tabulated_rating_curve = p.tabulated_rating_curve.node_id, - pump = p.pump.node_id, - outlet = p.outlet.node_id, - user_demand_inflow = [ - p.user_demand.node_id[i] for i in eachindex(p.user_demand.node_id) for - _ in 1:length(p.user_demand.inflow_links[i]) - ], - user_demand_outflow = p.user_demand.node_id, - linear_resistance = p.linear_resistance.node_id, - manning_resistance = p.manning_resistance.node_id, - evaporation = p.basin.node_id, - infiltration = p.basin.node_id, - integral = p.pid_control.node_id, - ) +@kwdef struct BackTracking + c_1::Float64 = 1.0e-4 + ρ_hi::Float64 = 0.5 + ρ_lo::Float64 = 0.1 + iter_max::Int = 100 end -"Create the axis of the state vector" -function count_state_ranges(u_ids::StateTuple{Vector{NodeID}})::StateTuple{UnitRange{Int}} - return StateTuple{UnitRange{Int}}(ranges(map(length, collect(u_ids)))) -end +""" +Find the minimum of the loss function between 0 and 1 +by approximating the loss function by cubic polynomials +""" +function (ls::BackTracking)(loss) + (; c_1, iter_max, ρ_lo, ρ_hi) = ls -function build_state_vector(p_independent::ParametersIndependent) - # It is assumed that the horizontal flow states come first in - # p_independent.state_inflow_link and p_independent.state_outflow_link - (; state_ranges) = p_independent - u_ids = state_node_ids(p_independent) - data = zeros(length(p_independent.node_id)) - u = CVector(data, state_ranges) - # Ensure p_independent.node_id, state_ranges and u have the same length and order - ranges = (getproperty(state_ranges, x) for x in propertynames(state_ranges)) - @assert length(u) == length(p_independent.node_id) == mapreduce(length, +, ranges) - @assert keys(u_ids) == state_components - return u -end + loss_0 = loss(0.0) + ϵ = sqrt(eps()) + dloss_0 = (loss(ϵ) - loss_0) / ϵ -function build_reltol_vector(u0::CVector, reltol::Float64) - reltolv = fill(reltol, length(u0)) - mask = trues(length(u0)) - # Mask the non-cumulative states - for (node, range) in pairs(getaxes(u0)) - if node in (:integral,) - mask[range] .= false - end + success = true + + αₙ₋₁ = 1.0 + loss_αₙ₋₁ = loss(αₙ₋₁) + + iterfinite = 0 + iterfinitemax = -log2(eps(Float64)) + while !isfinite(loss_αₙ₋₁) && (iterfinite < iterfinitemax) + iterfinite += 1 + αₙ₋₁ /= 2 + loss_αₙ₋₁ = loss(αₙ₋₁) end - return reltolv, mask -end - -function reduce_state!(u_reduced, u, p_independent)::Nothing - (; basin, link_to_state_idx) = p_independent - (; inflow_ids, outflow_ids) = basin - (; combined_cumulative_flows) = u_reduced - state_ranges = getaxes(u) - u_reduced .= 0 - - for i in eachindex(basin.node_id) - basin_id = basin.node_id[i] - for inflow_id in inflow_ids[i] - # Flow on the link (inflow_id → basin). For UserDemand outflow this is - # the single user_demand_outflow state (1:1 per node). Link-based lookup - # correctly resolves per-link states if we ever get them upstream too. - state_idx = get_state_index( - state_ranges, - link_to_state_idx, - (inflow_id, basin_id), - ) - if isnothing(state_idx) - state_idx = get_state_index(state_ranges, inflow_id; inflow = false) - end - isnothing(state_idx) && continue - combined_cumulative_flows[i] += u[state_idx] + + αₙ = αₙ₋₁ + loss_αₙ = loss_αₙ₋₁ + + iterations = 0 + + # Check convergence criterion + while loss_αₙ > loss_0 + c_1 * αₙ * dloss_0 + iterations += 1 + if iterations > iter_max + # If the maximum number of iterations is exceeded, + # the line search has failed + success = false + break end - for outflow_id in outflow_ids[i] - # Flow on the link (basin → outflow_id). Must be link-based because a - # UserDemand can have multiple inflow-link states, one per source basin. - state_idx = get_state_index( - state_ranges, - link_to_state_idx, - (basin_id, outflow_id), - ) - if isnothing(state_idx) - state_idx = get_state_index(state_ranges, outflow_id; inflow = true) + if iterations == 1 + # For the first iteration, approximate the loss function with a quadratic + # polynomial based on loss(0), loss'(0) and loss(α_start) and find its minimum + αₙ₊₁ = -dloss_0 / (2 * (loss_αₙ₋₁ - loss_0 - dloss_0)) + else + # For the later iterations, approximate the loss function with a cubic + # polynomial based on loss(0), loss'(0), loss(αₙ₋₁) and loss(αₙ) and find its minimum + rhs_1 = loss_αₙ - loss_0 - dloss_0 * αₙ + rhs_2 = loss_αₙ₋₁ - loss_0 - dloss_0 * αₙ₋₁ + det = (αₙ * αₙ₋₁)^2 * (αₙ - αₙ₋₁) + a = (rhs_1 * αₙ₋₁^2 - rhs_2 * αₙ^2) / det + b = (rhs_2 * αₙ^3 - rhs_1 * αₙ₋₁^3) / det + + αₙ₊₁ = if isapprox(a, 0.0; atol = eps(Float64)) + -dloss_0 / (2 * b) + else + # Discriminant + d = max(0.0, b^2 - 3a * dloss_0) + (-b + sqrt(d)) / (3a) end - isnothing(state_idx) && continue - combined_cumulative_flows[i] -= u[state_idx] end - combined_cumulative_flows[i] -= u.evaporation[i] - combined_cumulative_flows[i] -= u.infiltration[i] + # Avoid reductions that are too big or too small + αₙ₊₁ = NaNMath.min(αₙ₊₁, αₙ * ρ_hi) + αₙ₊₁ = NaNMath.max(αₙ₊₁, αₙ * ρ_lo) + + αₙ₋₁ = αₙ + loss_αₙ₋₁ = loss_αₙ + + αₙ = αₙ₊₁ + loss_αₙ = loss(αₙ) end - u_reduced.integral .= u.integral - return nothing + # As a last check, say the line search failed when the loss value actually increased + if success && (loss_αₙ > loss(1.0)) + success = false + end + + # If the line search failed, this NaN will propagate through OrdinaryDiffEq + # and causes the time step to be rejected + return success ? αₙ : NaN end """ -Create vectors state_inflow_link and state_outflow_link which give for each state -in the state vector in order the metadata of the link that is associated with that state. -Only for horizontal flows, which are assumed to come first in the state vector. +Compute the residual of the non-linear solver, i.e. a measure of the +error in the solution to the implicit equation defined by the solver algorithm """ -function get_state_flow_links( - graph::MetaGraph, - nodes::NamedTuple, - )::Tuple{Vector{LinkMetadata}, Vector{LinkMetadata}} - (; user_demand) = nodes - state_inflow_link = LinkMetadata[] - state_outflow_link = LinkMetadata[] - - placeholder_link = - LinkMetadata(0, LinkType.flow, (NodeID(:Terminal, 0, 0), NodeID(:Terminal, 0, 0))) - - for node_name in state_components - if hasproperty(nodes, node_name) - node::AbstractParameterNode = getproperty(nodes, node_name) - for id in node.node_id - inflow_ids_ = collect(inflow_ids(graph, id)) - outflow_ids_ = collect(outflow_ids(graph, id)) - - inflow_link = if length(inflow_ids_) == 0 - placeholder_link - elseif length(inflow_ids_) == 1 - inflow_id = only(inflow_ids_) - graph[inflow_id, id] - else - error("Multiple inflows not supported") - end - push!(state_inflow_link, inflow_link) +function residual(z::CVector, integrator, nlsolver, f::TF) where {TF} + (; uprev, t, p, dt, opts, isdae) = integrator + (; tmp, ztmp, γ, α, cache, method) = nlsolver + (; ustep, atmp, tstep, k, invγdt) = cache + @assert !isdae + b, ustep2 = _compute_rhs!(tmp, ztmp, ustep, γ, α, tstep, k, invγdt, method, p, dt, f, z) + calculate_residuals!( + atmp, + b, + uprev, + ustep2, + opts.abstol, + opts.reltol, + opts.internalnorm, + t, + ) + ndz = opts.internalnorm(atmp, t) + return ndz +end - outflow_link = if length(outflow_ids_) == 0 - placeholder_link - elseif length(outflow_ids_) == 1 - outflow_id = only(outflow_ids_) - graph[id, outflow_id] - else - error("Multiple outflows not supported") - end - push!(state_outflow_link, outflow_link) - end - elseif startswith(String(node_name), "user_demand") - if node_name == :user_demand_inflow - # One state per (UserDemand, inflow link) pair. - for links in user_demand.inflow_links - for link_meta in links - push!(state_inflow_link, link_meta) - push!(state_outflow_link, placeholder_link) - end - end - elseif node_name == :user_demand_outflow - placeholder_links = fill(placeholder_link, length(user_demand.node_id)) - append!(state_inflow_link, placeholder_links) - append!(state_outflow_link, user_demand.outflow_link) - end +function OrdinaryDiffEqNonlinearSolve.relax!( + dz, + nlsolver::AbstractNLSolver, + integrator::DEIntegrator, + f, + linesearch::BackTracking, + ) + atmp = nlsolver.cache.atmp + + let dz = dz, + integrator = integrator, + nlsolver = nlsolver, + f = f, + linesearch = linesearch + + function loss(α) + @. atmp = nlsolver.z - dz * α + return residual(atmp, integrator, nlsolver, f) end - end - return state_inflow_link, state_outflow_link -end + α = linesearch(loss) + @. dz *= α -""" -Build a mapping from a (from_node, to_node) link tuple to the index of that link's state. -Only inflow-link states are covered (horizontal flow components, which come first in the -state vector). Placeholder links (both node values == 0) are skipped. -""" -function build_link_to_state_idx( - state_inflow_link::Vector{LinkMetadata}, - )::Dict{Tuple{NodeID, NodeID}, Int} - link_to_state_idx = Dict{Tuple{NodeID, NodeID}, Int}() - for (idx, link_meta) in enumerate(state_inflow_link) - from_node, to_node = link_meta.link - from_node.value == 0 && to_node.value == 0 && continue - link_to_state_idx[link_meta.link] = idx + return dz end - return link_to_state_idx end -""" -Get the index of the state vector corresponding to the given NodeID. -Use the inflow Boolean argument to disambiguite for node types that have multiple states. -Can return nothing for node types that do not have a state, like Terminal. -""" -function get_state_index( - state_ranges::StateTuple{UnitRange{Int}}, - id::NodeID; - inflow::Bool = true, - )::Union{Int, Nothing} - component_name = if id.type == NodeType.UserDemand - inflow ? :user_demand_inflow : :user_demand_outflow - else - snake_case(id) - end +"Create a NamedTuple of the node IDs per state component in the state order" +function state_node_ids( + p::Union{ParametersIndependent, NamedTuple}, + )::StateTuple{Vector{NodeID}} + return (; + basin = p.basin.node_id, + integral = p.pid_control.node_id, + ) +end - if hasproperty(state_ranges, component_name) - state_range = getproperty(state_ranges, component_name) - return state_range[id.idx] - else - return nothing - end +"Create the axis of the state vector" +function count_state_ranges(u_ids::StateTuple{Vector{NodeID}})::StateTuple{UnitRange{Int}} + return StateTuple{UnitRange{Int}}(ranges(map(length, collect(u_ids)))) end -""" -Get the state index for a flow link. +function build_state_vector(p_independent::ParametersIndependent) + (; state_ranges) = p_independent + u_ids = state_node_ids(p_independent) + n_states = sum(length, u_ids) + data = zeros(n_states) + u = CVector(data, state_ranges) + @assert keys(u_ids) == state_components + return u +end -When the destination node has multiple inflow-link states (UserDemand with multiple -source links), the per-link index must be resolved via `link_to_state_idx`. Otherwise -this falls back to the to-node's (or from-node's) state. -""" -function get_state_index( - state_ranges::StateTuple{UnitRange{Int}}, - link_to_state_idx::Dict{Tuple{NodeID, NodeID}, Int}, - link::Tuple{NodeID, NodeID}, - )::Union{Int, Nothing} - idx = get(link_to_state_idx, link, nothing) - isnothing(idx) || return idx - idx = get_state_index(state_ranges, link[2]) - return isnothing(idx) ? get_state_index(state_ranges, link[1]; inflow = false) : idx +"Initialize the state vector with basin storages" +function initialize_state_vector!(u::CVector, p_independent::ParametersIndependent) + (; basin) = p_independent + u.basin .= basin.storage0 + # Initialize prev-saveat storage so the first water balance check + # correctly computes storage_rate = (s_now - storage0) / Δt + basin.Δstorage_prev_saveat .= basin.storage0 + # Initialize balance correction storage tracking + p_independent.balance_correction.storage_prev .= basin.storage0 + return nothing end """ Check whether any storages are negative given the state u. +Storage states are directly in u.basin. """ function isoutofdomain(u, p, t) - (; current_storage) = p.state_and_time_dependent_cache - (; u_reduced) = p.p_independent - reduce_state!(u_reduced, u, p.p_independent) - formulate_storages!(u_reduced, p, t) - return any(<(0), current_storage) + return any(<(0), u.basin) end function get_demand(user_demand, id, demand_priority_idx, t)::Float64 @@ -1174,7 +1143,7 @@ Check whether the inputs u and t are different from the previous call of water_b update the boolean flags in p_mutable. In several parts of the calculations in water_balance!, caches are only updated if the data they depend on is different from the previous water_balance! call. """ -function check_new_input!(p::Parameters, u_reduced::CVector, t::Number)::Nothing +function check_new_input!(p::Parameters, u::CVector, t::Number)::Nothing (; state_and_time_dependent_cache, time_dependent_cache, p_mutable) = p # Whether the time dependent cache must be renewed @@ -1195,18 +1164,18 @@ function check_new_input!(p::Parameters, u_reduced::CVector, t::Number)::Nothing ) new_u_state_and_time_dependent_cache = any( - i -> !isassigned(state_and_time_dependent_cache.u_reduced_prev_call, i), - eachindex(u_reduced), + i -> !isassigned(state_and_time_dependent_cache.u_prev_call, i), + eachindex(u), ) || any( i -> !( - u_reduced[i] == state_and_time_dependent_cache.u_reduced_prev_call[i] && - ForwardDiff.partials(u_reduced[i]) == ForwardDiff.partials( - state_and_time_dependent_cache.u_reduced_prev_call[i], + u[i] == state_and_time_dependent_cache.u_prev_call[i] && + ForwardDiff.partials(u[i]) == ForwardDiff.partials( + state_and_time_dependent_cache.u_prev_call[i], ) ), - eachindex(u_reduced), + eachindex(u), ) - state_and_time_dependent_cache.u_reduced_prev_call .= u_reduced + state_and_time_dependent_cache.u_prev_call .= u state_and_time_dependent_cache.t_prev_call[1] = t p_mutable.new_state_and_time_dependent_cache = new_t_state_and_time_dependent_cache || new_u_state_and_time_dependent_cache @@ -1336,3 +1305,17 @@ function add_substance_mass!( end return nothing end + +function get_link_index( + link::Tuple{NodeID, NodeID}, + flow_links::Vector{LinkMetadata}, + )::Union{Int64, Nothing} + return findfirst(l -> l.link == link, flow_links) +end + +function get_link_index( + link::Tuple{NodeID, NodeID}, + flow_link_lookup::Dict{Tuple{NodeID, NodeID}, Int}, + )::Union{Int64, Nothing} + return get(flow_link_lookup, link, nothing) +end diff --git a/core/src/write.jl b/core/src/write.jl index 4d3861cf6..9d75da58f 100644 --- a/core/src/write.jl +++ b/core/src/write.jl @@ -344,9 +344,7 @@ function basin_state_data(model::Model; table::Bool = true) (; current_level) = p.state_and_time_dependent_cache # ensure the levels are up-to-date - (; u_reduced) = p.p_independent - reduce_state!(u_reduced, u, p.p_independent) - set_current_basin_properties!(u_reduced, p, t) + set_current_basin_properties!(u, p, t) return (; node_id = Int32.(p.p_independent.basin.node_id), level = current_level) end @@ -354,8 +352,6 @@ end "Create the basin result table from the saved data" function basin_data(model::Model; table::Bool = true) (; saved) = model - (; u) = model.integrator - state_ranges = getaxes(u) # The last timestep is not included; there is no period over which to compute flows. data = get_storages_and_levels(model) @@ -369,25 +365,14 @@ function basin_data(model::Model; table::Bool = true) inflow_rate = FlatVector(saved.flow.saveval, :inflow) outflow_rate = FlatVector(saved.flow.saveval, :outflow) drainage = FlatVector(saved.flow.saveval, :drainage) - infiltration = zeros(nrows) - evaporation = zeros(nrows) + evaporation = FlatVector(saved.flow.saveval, :evaporation) + infiltration = FlatVector(saved.flow.saveval, :infiltration) precipitation = FlatVector(saved.flow.saveval, :precipitation) surface_runoff = FlatVector(saved.flow.saveval, :surface_runoff) storage_rate = FlatVector(saved.flow.saveval, :storage_rate) balance_error = FlatVector(saved.flow.saveval, :balance_error) relative_error = FlatVector(saved.flow.saveval, :relative_error) - convergence = FlatVector(saved.flow.saveval, :basin_convergence) - - idx_row = 0 - for saved_flow in saved.flow.saveval - saved_evaporation = view(saved_flow.flow, state_ranges.evaporation) - saved_infiltration = view(saved_flow.flow, state_ranges.infiltration) - for (evaporation_, infiltration_) in zip(saved_evaporation, saved_infiltration) - idx_row += 1 - evaporation[idx_row] = evaporation_ - infiltration[idx_row] = infiltration_ - end - end + convergence = FlatVector(saved.flow.saveval, :convergence) time = data.time[begin:(end - 1)] node_id = Int32.(data.node_id) @@ -442,7 +427,7 @@ end function flow_data(model::Model; table::Bool = true) (; config, saved, integrator) = model (; t, saveval) = saved.flow - (; u, p) = integrator + (; p) = integrator (; p_independent) = p (; graph) = p_independent (; internal_flow_links, external_flow_links, flow_link_map) = graph[] @@ -460,19 +445,16 @@ function flow_data(model::Model; table::Bool = true) nflow = length(unique_link_ids_flow) ntsteps = length(t) flow_rate = zeros(nflow * ntsteps) - flow_rate_conv = zeros(Union{Missing, Float64}, nflow * ntsteps) + flow_rate_conv = zeros(nflow * ntsteps) internal_flow_rate = zeros(length(internal_flow_links)) - internal_flow_rate_conv = zeros(Union{Missing, Float64}, length(internal_flow_links)) + internal_flow_rate_conv = zeros(length(internal_flow_links)) - for (ti, cvec) in enumerate(saveval) - (; flow, flow_boundary, flow_convergence) = cvec - flow = CVector(flow, getaxes(u)) - convergence = CVector(flow_convergence, getaxes(u)) + for (ti, saved_flow) in enumerate(saveval) + (; flow, flow_boundary, flow_convergence) = saved_flow for (fi, link) in enumerate(internal_flow_links) internal_flow_rate[fi] = get_flow(flow, p_independent, 0.0, link.link; boundary_flow = flow_boundary) - - internal_flow_rate_conv[fi] = get_convergence(convergence, link.link) + internal_flow_rate_conv[fi] = flow_convergence[fi] end mul!( view(flow_rate, (1 + (ti - 1) * nflow):(ti * nflow)), diff --git a/core/test/allocation_physics_test.jl b/core/test/allocation_physics_test.jl index 4340fa618..1008d298b 100644 --- a/core/test/allocation_physics_test.jl +++ b/core/test/allocation_physics_test.jl @@ -98,7 +98,6 @@ end filter!(:link_id => ==(1), allocation_flow_table) filter!(:link_id => ==(1), flow_table) - @test allocation_flow_table.flow_rate ≈ flow_table.flow_rate atol = 8.0e-4 end diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 1bc152e6a..cbb0c9113 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -71,7 +71,7 @@ end # In this section the Basin leaves no supply for the UserDemand stage_1 = t .≤ 2Δt_allocation u_stage_1(τ) = storage[1] + (q + ϕ) * τ - @test storage[stage_1] ≈ u_stage_1.(t[stage_1]) rtol = 1.0e-10 + @test storage[stage_1] ≈ u_stage_1.(t[stage_1]) rtol = 1.0e-5 # In this section the Basin gets exactly what it needs to get to the target min # level of 1 m (equivalent to 1000 m^3) @@ -79,7 +79,7 @@ end u_stage_2(τ) = (3Δt_allocation - τ) / Δt_allocation * u_stage_1(2Δt_allocation) + min_storage * (τ - 2Δt_allocation) / Δt_allocation - @test storage[stage_2] ≈ u_stage_2.(t[stage_2]) rtol = 1.0e-10 + @test storage[stage_2] ≈ u_stage_2.(t[stage_2]) rtol = 1.0e-5 # In this section (and following sections) the basin has no longer a (positive) demand, # since precipitation provides enough water to get the basin to its target level @@ -87,14 +87,14 @@ end stage_3 = 3Δt_allocation .≤ t .≤ 15Δt_allocation stage_3_start_idx = findfirst(stage_3) u_stage_3(τ) = min_storage + (ϕ + q - d) * (τ - t[stage_3_start_idx]) - @test storage[stage_3] ≈ u_stage_3.(t[stage_3]) rtol = 1.0e-10 + @test storage[stage_3] ≈ u_stage_3.(t[stage_3]) rtol = 1.0e-5 # At the start of this section precipitation stops, and so the UserDemand # partly uses surplus water from the basin to fulfill its demand stage_4 = 15Δt_allocation .≤ t .≤ 27Δt_allocation stage_4_start_idx = findfirst(stage_4) u_stage_4(τ) = storage[stage_4_start_idx] + (q - d) * (τ - t[stage_4_start_idx]) - @test storage[stage_4] ≈ u_stage_4.(t[stage_4]) rtol = 1.0e-10 + @test storage[stage_4] ≈ u_stage_4.(t[stage_4]) rtol = 1.0e-5 # From this point the basin is in a dynamical equilibrium, # since the basin has no supply so the UserDemand abstracts precisely @@ -102,7 +102,7 @@ end stage_5 = 27Δt_allocation .<= t stage_5_start_idx = findfirst(stage_5) u_stage_5(τ) = min_storage - @test storage[stage_5] ≈ u_stage_5.(t[stage_5]) rtol = 1.0e-10 + @test storage[stage_5] ≈ u_stage_5.(t[stage_5]) rtol = 1.0e-5 # Isolated LevelDemand + Basin pair to test optional min_level (; problem) = allocation.allocation_models[2] @@ -137,7 +137,7 @@ end seconds_since.(df_user_3.time, model.config.starttime), ), ) ./ Δt_allocation - @test all(isapprox.(supplied_numeric[3:end], df_user_3.supplied[4:end], atol = 1.0e-3)) + @test all(isapprox.(supplied_numeric[3:end], df_user_3.supplied[4:end], atol = 1.0e-10)) end @testitem "Flow demand" setup = [Teamcity] begin @@ -406,7 +406,7 @@ end for (link_id, flow) in zip([2, 4, 6], [flow_1, flow_2, flow_3]) data = filter(:link_id => ==(link_id), flow_table) - @test all(isapprox.(data.flow_rate, flow[1:(end - 1)], atol = 1.0e-5)) + @test all(isapprox.(data.flow_rate, flow[1:(end - 1)], atol = 1.0e-2)) end end diff --git a/core/test/bmi_test.jl b/core/test/bmi_test.jl index 49aca12cb..e5a33806e 100644 --- a/core/test/bmi_test.jl +++ b/core/test/bmi_test.jl @@ -5,7 +5,7 @@ toml_path = normpath(@__DIR__, "../../generated_testmodels/basic/ribasim.toml") model = BMI.initialize(Ribasim.Model, toml_path) @test BMI.get_time_units(model) == "s" - dt0 = 2.8280652f-5 + dt0 = 1.07127 @test BMI.get_time_step(model) ≈ dt0 atol = 5.0e-3 @test BMI.get_start_time(model) === 0.0 @test BMI.get_current_time(model) === 0.0 @@ -54,7 +54,7 @@ end toml_path = normpath(@__DIR__, "../../generated_testmodels/basic/ribasim.toml") model = BMI.initialize(Ribasim.Model, toml_path) storage0 = BMI.get_value_ptr(model, "basin.storage") - @test storage0 ≈ ones(4) + @test storage0 ≈ ones(4) atol = 2.0e-2 @test_throws "Unknown variable foo" BMI.get_value_ptr(model, "foo") BMI.update_until(model, 86400.0) storage = BMI.get_value_ptr(model, "basin.storage") @@ -101,7 +101,7 @@ end inflow = BMI.get_value_ptr(model, "user_demand.cumulative_inflow") day = 86400.0 BMI.update_until(model, 2day) - @test inflow ≈ [2.0e-3 * day, 3.0e-3 * day] atol = 1.0e-2 + @test inflow ≈ [2.0e-3 * day, 3.0e-3 * day] demand[1] = 3.0e-3 BMI.update_until(model, 3day) @test inflow[1] ≈ 5.0e-3 * day atol = 2.0e-3 diff --git a/core/test/data/allocation_problems/linear_resistance_demand/allocation_problem_2.lp b/core/test/data/allocation_problems/linear_resistance_demand/allocation_problem_2.lp index 6028e5b9b..be5831422 100644 --- a/core/test/data/allocation_problems/linear_resistance_demand/allocation_problem_2.lp +++ b/core/test/data/allocation_problems/linear_resistance_demand/allocation_problem_2.lp @@ -10,14 +10,14 @@ linear_resistance_constraint_LinearResistance_#2_: -0.01 basin_storage_change_Ba flow_demand_allocated_sum_constraint_LinearResistance_#2_: 1 flow_(Basin_#1,_LinearResistance_#2)_ - 1 flow_demand_allocated_LinearResistance_#2,1_ - 1 flow_demand_extra_LinearResistance_#2_ = 0 average_flow_unit_error_constraint_1_: -1 flow_demand_error_LinearResistance_#2,1,first_ + 1 average_flow_unit_error_1_ = 0 Bounds --2 <= basin_storage_change_Basin_#1_ <= 8 --2 <= basin_storage_change_Basin_#3_ <= 8 +-1.8181818181818181 <= basin_storage_change_Basin_#1_ <= 7.2727272727272725 +-1.8181818181818181 <= basin_storage_change_Basin_#3_ <= 7.2727272727272725 0 <= low_storage_factor_Basin_#1_ <= 1 0 <= low_storage_factor_Basin_#3_ <= 1 --345.6 <= flow_(Basin_#1,_LinearResistance_#2)_ <= 345.6 --345.6 <= flow_(LinearResistance_#2,_Basin_#3)_ <= 345.6 --86400000 <= flow_demand_allocated_LinearResistance_#2,1_ <= 2 -0 <= flow_demand_extra_LinearResistance_#2_ <= 86400000 +-314.1818181818182 <= flow_(Basin_#1,_LinearResistance_#2)_ <= 314.1818181818182 +-314.1818181818182 <= flow_(LinearResistance_#2,_Basin_#3)_ <= 314.1818181818182 +-7.854545454545455e7 <= flow_demand_allocated_LinearResistance_#2,1_ <= 2 +0 <= flow_demand_extra_LinearResistance_#2_ <= 7.854545454545455e7 0 <= flow_demand_error_LinearResistance_#2,1,first_ <= 1 0 <= flow_demand_error_LinearResistance_#2,1,second_ <= 1 0 <= average_flow_unit_error_1_ <= 1 diff --git a/core/test/docs.toml b/core/test/docs.toml index 40be3c27b..e475c242a 100644 --- a/core/test/docs.toml +++ b/core/test/docs.toml @@ -45,7 +45,7 @@ dt = 60.0 # optional, remove for adaptive time stepping dtmin = 0.0 # optional, default 0.0 dtmax = 0.0 # optional, default length of simulation force_dtmin = false # optional, default false -abstol = 1e-5 # optional, default 1e-5 +abstol = 1e-4 # optional, default 1e-4 reltol = 1e-5 # optional, default 1e-5 water_balance_abstol = 1e-3 # optional, default 1e-3 water_balance_reltol = 1e-2 # optional, default 1e-2 diff --git a/core/test/io_test.jl b/core/test/io_test.jl index 47ba37a74..7bb73bf96 100644 --- a/core/test/io_test.jl +++ b/core/test/io_test.jl @@ -541,5 +541,5 @@ end (; p_independent, state_and_time_dependent_cache) = model.integrator.p (; current_storage) = state_and_time_dependent_cache storage2_begin = current_storage - @test storage1_end ≈ storage2_begin + @test storage1_end ≈ storage2_begin rtol = 1.0e-2 end diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index ec7860d02..34921a852 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -14,10 +14,8 @@ (; u, du) = model.integrator (; p_independent) = model.integrator.p - @test p_independent.node_id == [0, 6, 6] @test u isa CVector - @test filter(!isempty, getaxes(u)) == - (; tabulated_rating_curve = 1:1, evaporation = 2:2, infiltration = 3:3) + @test filter(!isempty, getaxes(u)) == (; basin = 1:1) # Open NetCDF result files flow_path = normpath(dirname(toml_path), "results/flow.nc") @@ -138,8 +136,6 @@ end @test basin.vertical_flux.precipitation == [0.0] @test basin.vertical_flux.drainage == [0.0] du = get_du(model.integrator) - @test du.evaporation == [0.0] - @test du.infiltration == [0.0] @test success(model) end @@ -164,9 +160,9 @@ end Ribasim.water_balance!(du, u, p, t) stor = state_and_time_dependent_cache.current_storage prec = basin.vertical_flux.precipitation - evap = du.evaporation + evap = basin.vertical_flux.potential_evaporation drng = basin.vertical_flux.drainage - infl = du.infiltration + infl = basin.vertical_flux.infiltration # The dynamic data has missings, but these are not set. @test prec == [0.0] @test evap == [0.0] @@ -178,13 +174,13 @@ end @test prec == [0.0] @test evap == [0.0] @test drng == [0.003] - @test infl == [0.0] + @test infl == [0.001] stor ≈ Float32[init_stor + 86400 * (0.003 * 1.5 - 0.001 * 0.5)] BMI.update_until(model, 2.5 * 86400) @test prec == [0.0] @test evap == [0.0] @test drng == [0.001] - @test infl == [0.0] + @test infl == [0.002] stor ≈ Float32[init_stor + 86400 * (0.003 * 2.0 + 0.001 * 0.5 - 0.001 - 0.002 * 0.5)] @test success(Ribasim.solve!(model)) end @@ -214,10 +210,10 @@ end @test p isa Ribasim.Parameters @test isconcretetype(typeof(p_independent)) @test all(isconcretetype, fieldtypes(typeof(p_independent))) - @test p_independent.node_id == [4, 5, 8, 7, 10, 12, 2, 1, 3, 6, 9, 1, 3, 6, 9] + # node_id field was removed; state vector now only has basin and integral axes @test success(model) - @test length(model.integrator.sol) == 2 # start and end + @test length(model.integrator.sol.t) == 2 # start and end @test state_and_time_dependent_cache.current_storage ≈ Float32[775.23576, 775.23365, 572.60102, 1130.005] skip = Sys.isapple() atol = 1.5 @@ -388,7 +384,6 @@ end @testitem "UserDemand" begin using Dates using DataFrames: DataFrame - using Ribasim: formulate_storages! import BasicModelInterface as BMI toml_path = normpath(@__DIR__, "../../generated_testmodels/user_demand/ribasim.toml") @@ -401,16 +396,12 @@ end day = 86400.0 - @test only(state_and_time_dependent_cache.current_storage) ≈ 1000.0 + @test only(u.basin) ≈ 1000.0 # constant UserDemand withdraws to 0.9m or 900m3 due to min level = 0.9 BMI.update_until(model, 150day) - (; u_reduced) = p.p_independent - Ribasim.reduce_state!(u_reduced, u, p_independent) - formulate_storages!(u_reduced, p, t) @test only(state_and_time_dependent_cache.current_storage) ≈ 900 atol = 5 # dynamic UserDemand withdraws to 0.5m or 500m3 due to min level = 0.5 BMI.update_until(model, 200day) - formulate_storages!(u_reduced, p, t) @test only(state_and_time_dependent_cache.current_storage) ≈ 500 atol = 2 # Transient return factor diff --git a/core/test/time_test.jl b/core/test/time_test.jl index 5e08d968b..e8fa3ae0b 100644 --- a/core/test/time_test.jl +++ b/core/test/time_test.jl @@ -111,17 +111,6 @@ end @test length(only(tstops)) == 404 end -@testitem "decrease tolerance" begin - toml_path = normpath(@__DIR__, "../../generated_testmodels/cyclic_time/ribasim.toml") - @test ispath(toml_path) - - model = Ribasim.run(toml_path) - @test model.integrator.opts.reltol isa Vector{Float64} - @test all(model.integrator.opts.reltol .<= model.integrator.p.p_independent.reltol) - @test model.integrator.u[1] >= 1.0e11 - @test model.integrator.opts.reltol[1] <= 1.0e-11 -end - @testitem "transient_pump_outlet" begin using DataFrames: DataFrame @@ -133,7 +122,7 @@ end storage = Ribasim.get_storages_and_levels(model).storage flowrate = Ribasim.get_flow @test storage[1, 1] ≈ 100.0f0 - @test storage[1, end] ≈ 110.0f0 + @test storage[1, end] ≈ 110.0f0 atol = 1.0e-2 # test the flowrate has been set correctly flow_data = DataFrame(Ribasim.flow_data(model)) diff --git a/core/test/utils_test.jl b/core/test/utils_test.jl index bb8918c59..2897ce709 100644 --- a/core/test/utils_test.jl +++ b/core/test/utils_test.jl @@ -294,7 +294,6 @@ end @testitem "Solver algorithm" begin using LinearSolve: KLUFactorization - using OrdinaryDiffEqNonlinearSolve: NLNewton using OrdinaryDiffEqBDF: QNDF model = @@ -303,7 +302,7 @@ end @test alg isa QNDF @test alg.step_limiter! == Ribasim.limit_flow! - @test alg.nlsolve == NLNewton() + @test alg.nlsolve.relax isa Ribasim.BackTracking @test alg.linsolve == Ribasim.config.RibasimLinearSolve(KLUFactorization(; check_pattern = false)) end diff --git a/core/test/validation_test.jl b/core/test/validation_test.jl index a5b950a13..afa0a2e84 100644 --- a/core/test/validation_test.jl +++ b/core/test/validation_test.jl @@ -413,9 +413,9 @@ end "Warning: Convergence bottlenecks in descending order of severity:", output, ) - @test occursin("Pump #12 = ", output) - @test occursin("Pump #32 = ", output) - @test occursin("Pump #52 = ", output) + @test occursin("Basin #11 = ", output) + @test occursin("Basin #31 = ", output) + @test occursin("Basin #51 = ", output) end @testitem "Missing demand priority when allocation is active" begin diff --git a/python/ribasim/ribasim/config.py b/python/ribasim/ribasim/config.py index 597476995..0475f5f77 100644 --- a/python/ribasim/ribasim/config.py +++ b/python/ribasim/ribasim/config.py @@ -135,9 +135,9 @@ class Solver(ChildModel): If a smaller dt than dtmin is needed to meet the set error tolerances, the simulation stops, unless force_dtmin = true (Optional, defaults to False) abstol : float - The absolute tolerance for adaptive timestepping (Optional, defaults to 1e-7) + The absolute tolerance for adaptive timestepping (Optional, defaults to 1e-4) reltol : float - The relative tolerance for adaptive timestepping (Optional, defaults to 1e-7) + The relative tolerance for adaptive timestepping (Optional, defaults to 1e-5) maxiters : int The total number of linear iterations over the whole simulation. (Defaults to 1e9, only needs to be increased for extremely long simulations) sparse : bool @@ -160,7 +160,7 @@ class Solver(ChildModel): dtmin: float | None = None dtmax: float | None = None force_dtmin: bool = False - abstol: float = 1e-06 + abstol: float = 1e-04 reltol: float = 1e-05 maxiters: int = 1000000000 sparse: bool = True diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index bdf8b82f6..9d32c91ce 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -760,11 +760,11 @@ def linear_resistance_demand_model(): model.basin.add( Node(1, Point(0, 0), subnetwork_id=2), - [basin.Profile(area=1e3, level=[0.0, 1.0]), basin.State(level=[1.0])], + [basin.Profile(area=1e3, level=[0.0, 1.1]), basin.State(level=[1.0])], ) model.basin.add( Node(3, Point(2, 0), subnetwork_id=2), - [basin.Profile(area=1e3, level=[0.0, 1.0]), basin.State(level=[1.0])], + [basin.Profile(area=1e3, level=[0.0, 1.1]), basin.State(level=[1.0])], ) model.linear_resistance.add(