diff --git a/Manifest.toml b/Manifest.toml index 8abd724..b47b133 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -12,15 +12,33 @@ version = "3.3.3" [[deps.ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +[[deps.ArnoldiMethod]] +deps = ["LinearAlgebra", "Random", "StaticArrays"] +git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae" +uuid = "ec485272-7323-5ecc-a04f-4719b315124d" +version = "0.2.0" + [[deps.ArrayInterface]] deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"] -git-tree-sha1 = "1d6835607e9f214cb4210310868f8cf07eb0facc" +git-tree-sha1 = "1ee88c4c76caa995a885dc2f22a5d548dfbbc0ba" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "3.1.34" +version = "3.2.2" + +[[deps.ArrayLayouts]] +deps = ["FillArrays", "LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "56c347caf09ad8acb3e261fe75f8e09652b7b05b" +uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" +version = "0.7.10" [[deps.Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +[[deps.BandedMatrices]] +deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "Random", "SparseArrays"] +git-tree-sha1 = "ce68f8c2162062733f9b4c9e3700d5efc4a8ec47" +uuid = "aae01518-5342-5314-be14-df237901396f" +version = "0.16.11" + [[deps.Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" @@ -30,6 +48,17 @@ git-tree-sha1 = "5e98d6a6aa92e5758c4d58501b7bf23732699fa3" uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b" version = "0.1.2" +[[deps.BoundaryValueDiffEq]] +deps = ["BandedMatrices", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "NLsolve", "Reexport", "SparseArrays"] +git-tree-sha1 = "fe34902ac0c3a35d016617ab7032742865756d7d" +uuid = "764a87c0-6b3e-53db-9096-fe964310641d" +version = "2.7.1" + +[[deps.CEnum]] +git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9" +uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" +version = "0.4.1" + [[deps.CPUSummary]] deps = ["Hwloc", "IfElse", "Static"] git-tree-sha1 = "849799453de85b55e78550fc7b0c8f442eb497ab" @@ -44,9 +73,9 @@ version = "0.5.1" [[deps.ChainRulesCore]] deps = ["Compat", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "7dd38532a1115a215de51775f9891f0f3e1bac6a" +git-tree-sha1 = "c9a6160317d1abe9c44b3beb367fd448117679ca" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.12.1" +version = "1.13.0" [[deps.ChangesOfVariables]] deps = ["ChainRulesCore", "LinearAlgebra", "Test"] @@ -60,6 +89,17 @@ git-tree-sha1 = "03dc838350fbd448fca0b99285ed4d60fc229b72" uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9" version = "0.1.5" +[[deps.CommonSolve]] +git-tree-sha1 = "68a0743f578349ada8bc911a5cbd5a2ef6ed6d1f" +uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" +version = "0.2.0" + +[[deps.CommonSubexpressions]] +deps = ["MacroTools", "Test"] +git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" +uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" +version = "0.3.0" + [[deps.Compat]] deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582" @@ -81,6 +121,12 @@ git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15" uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" version = "4.1.1" +[[deps.DEDataArrays]] +deps = ["ArrayInterface", "DocStringExtensions", "LinearAlgebra", "RecursiveArrayTools", "SciMLBase", "StaticArrays"] +git-tree-sha1 = "31186e61936fbbccb41d809ad4338c9f7addf7ae" +uuid = "754358af-613d-5f8d-9788-280bf1605d4c" +version = "0.2.0" + [[deps.DataAPI]] git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8" uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" @@ -107,6 +153,12 @@ version = "1.0.0" deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +[[deps.DelayDiffEq]] +deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "LinearAlgebra", "Logging", "NonlinearSolve", "OrdinaryDiffEq", "Printf", "RecursiveArrayTools", "Reexport", "UnPack"] +git-tree-sha1 = "ceb3463f2913eec2f0af5f0d8e1386fb546fdd32" +uuid = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb" +version = "5.34.0" + [[deps.DelimitedFiles]] deps = ["Mmap"] uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" @@ -117,12 +169,54 @@ git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b" uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d" version = "0.4.0" +[[deps.DiffEqBase]] +deps = ["ArrayInterface", "ChainRulesCore", "DEDataArrays", "DataStructures", "Distributions", "DocStringExtensions", "FastBroadcast", "ForwardDiff", "FunctionWrappers", "IterativeSolvers", "LabelledArrays", "LinearAlgebra", "Logging", "MuladdMacro", "NonlinearSolve", "Parameters", "PreallocationTools", "Printf", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "Requires", "SciMLBase", "Setfield", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "ZygoteRules"] +git-tree-sha1 = "433291c9e63dcfc1a0e42c6aeb6bb5d3e5ab1789" +uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" +version = "6.81.4" + +[[deps.DiffEqCallbacks]] +deps = ["DataStructures", "DiffEqBase", "ForwardDiff", "LinearAlgebra", "NLsolve", "OrdinaryDiffEq", "Parameters", "RecipesBase", "RecursiveArrayTools", "SciMLBase", "StaticArrays"] +git-tree-sha1 = "e57ecaf9f7875714c164ccca3c802711589127cf" +uuid = "459566f4-90b8-5000-8ac3-15dfb0a30def" +version = "2.20.1" + +[[deps.DiffEqJump]] +deps = ["ArrayInterface", "Compat", "DataStructures", "DiffEqBase", "FunctionWrappers", "Graphs", "LinearAlgebra", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "StaticArrays", "TreeViews", "UnPack"] +git-tree-sha1 = "e30f058eb600407e3fd4ea082e2527e3a3671238" +uuid = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12" +version = "8.2.1" + +[[deps.DiffEqNoiseProcess]] +deps = ["DiffEqBase", "Distributions", "LinearAlgebra", "Optim", "PoissonRandom", "QuadGK", "Random", "Random123", "RandomNumbers", "RecipesBase", "RecursiveArrayTools", "Requires", "ResettableStacks", "SciMLBase", "StaticArrays", "Statistics"] +git-tree-sha1 = "d6839a44a268c69ef0ed927b22a6f43c8a4c2e73" +uuid = "77a26b50-5914-5dd7-bc55-306e6241c503" +version = "5.9.0" + +[[deps.DiffResults]] +deps = ["StaticArrays"] +git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805" +uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +version = "1.0.3" + [[deps.DiffRules]] deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] git-tree-sha1 = "dd933c4ef7b4c270aacd4eb88fa64c147492acf0" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" version = "1.10.0" +[[deps.DifferentialEquations]] +deps = ["BoundaryValueDiffEq", "DelayDiffEq", "DiffEqBase", "DiffEqCallbacks", "DiffEqJump", "DiffEqNoiseProcess", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEq", "Random", "RecursiveArrayTools", "Reexport", "SteadyStateDiffEq", "StochasticDiffEq", "Sundials"] +git-tree-sha1 = "3f3db9365fedd5fdbecebc3cce86dfdfe5c43c50" +uuid = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +version = "7.1.0" + +[[deps.Distances]] +deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"] +git-tree-sha1 = "3258d0659f812acde79e8a74b11f17ac06d0ca04" +uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +version = "0.10.7" + [[deps.Distributed]] deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" @@ -143,16 +237,45 @@ version = "0.8.6" deps = ["ArgTools", "LibCURL", "NetworkOptions"] uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +[[deps.DualNumbers]] +deps = ["Calculus", "NaNMath", "SpecialFunctions"] +git-tree-sha1 = "84f04fe68a3176a583b864e492578b9466d87f1e" +uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" +version = "0.6.6" + +[[deps.ExponentialUtilities]] +deps = ["ArrayInterface", "LinearAlgebra", "Printf", "Requires", "SparseArrays"] +git-tree-sha1 = "3e1289d9a6a54791c1ee60da0850f4fd71188da6" +uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18" +version = "1.11.0" + [[deps.ExprTools]] git-tree-sha1 = "56559bbef6ca5ea0c0818fa5c90320398a6fbf8d" uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" version = "0.1.8" +[[deps.FastBroadcast]] +deps = ["LinearAlgebra", "Polyester", "Static"] +git-tree-sha1 = "0f8ef5dcb040dbb9edd98b1763ac10882ee1ff03" +uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +version = "0.1.12" + +[[deps.FastClosures]] +git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef" +uuid = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" +version = "0.3.2" + [[deps.FillArrays]] deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"] -git-tree-sha1 = "4c7d3757f3ecbcb9055870351078552b7d1dbd2d" +git-tree-sha1 = "deed294cde3de20ae0b2e0355a6c4e1c6a5ceffc" uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "0.13.0" +version = "0.12.8" + +[[deps.FiniteDiff]] +deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"] +git-tree-sha1 = "ec299fdc8f49ae450807b0cb1d161c6b76fd2b60" +uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" +version = "2.10.1" [[deps.Formatting]] deps = ["Printf"] @@ -160,10 +283,27 @@ git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8" uuid = "59287772-0a20-5a39-b81b-1366585eb4c0" version = "0.4.2" +[[deps.ForwardDiff]] +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"] +git-tree-sha1 = "1bd6fc0c344fc0cbee1f42f8d2e7ec8253dda2d2" +uuid = "f6369f11-7733-5829-9624-2563aa707210" +version = "0.10.25" + +[[deps.FunctionWrappers]] +git-tree-sha1 = "241552bc2209f0fa068b6415b1942cc0aa486bcc" +uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" +version = "1.1.2" + [[deps.Future]] deps = ["Random"] uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" +[[deps.Graphs]] +deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] +git-tree-sha1 = "57c021de207e234108a6f1454003120a1bf350c4" +uuid = "86223c79-3864-5bf0-83f7-82e725a168b6" +version = "1.6.0" + [[deps.HostCPUFeatures]] deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"] git-tree-sha1 = "3965a3216446a6b020f0d48f1ba94ef9ec01720d" @@ -182,11 +322,22 @@ git-tree-sha1 = "d8bccde6fc8300703673ef9e1383b11403ac1313" uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8" version = "2.7.0+0" +[[deps.HypergeometricFunctions]] +deps = ["DualNumbers", "LinearAlgebra", "SpecialFunctions", "Test"] +git-tree-sha1 = "65e4589030ef3c44d3b90bdc5aac462b4bb05567" +uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" +version = "0.3.8" + [[deps.IfElse]] git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" version = "0.1.1" +[[deps.Inflate]] +git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c" +uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" +version = "0.1.2" + [[deps.InteractiveUtils]] deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" @@ -207,6 +358,12 @@ git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151" uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" version = "0.1.1" +[[deps.IterativeSolvers]] +deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"] +git-tree-sha1 = "1169632f425f79429f245113b775a0e3d121457c" +uuid = "42fd0dbc-a981-5370-80f2-aaf504508153" +version = "0.9.2" + [[deps.IteratorInterfaceExtensions]] git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856" uuid = "82899510-4779-5014-852e-03e436cf321d" @@ -218,6 +375,36 @@ git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" version = "1.4.1" +[[deps.JSON]] +deps = ["Dates", "Mmap", "Parsers", "Unicode"] +git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "0.21.3" + +[[deps.KLU]] +deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse_jll"] +git-tree-sha1 = "cae5e3dfd89b209e01bcd65b3a25e74462c67ee0" +uuid = "ef3ab10e-7fda-4108-b977-705223b18434" +version = "0.3.0" + +[[deps.Krylov]] +deps = ["LinearAlgebra", "Printf", "SparseArrays"] +git-tree-sha1 = "6333cc5b848295895f3b23eb763d020fc8e05867" +uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" +version = "0.7.12" + +[[deps.KrylovKit]] +deps = ["LinearAlgebra", "Printf"] +git-tree-sha1 = "0328ad9966ae29ccefb4e1b9bfd8c8867e4360df" +uuid = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" +version = "0.5.3" + +[[deps.LabelledArrays]] +deps = ["ArrayInterface", "ChainRulesCore", "LinearAlgebra", "MacroTools", "StaticArrays"] +git-tree-sha1 = "3e6a4c07ea78db18f885e474c7de466ce257de85" +uuid = "2ee39098-c373-598a-b85f-a56591580800" +version = "1.7.3" + [[deps.LayoutPointers]] deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static"] git-tree-sha1 = "6dd77ee76188b0365f7d882d674b95796076fa2c" @@ -243,10 +430,22 @@ uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" [[deps.Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +[[deps.LineSearches]] +deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"] +git-tree-sha1 = "f27132e551e959b3667d8c93eae90973225032dd" +uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +version = "7.1.1" + [[deps.LinearAlgebra]] deps = ["Libdl", "libblastrampoline_jll"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +[[deps.LinearSolve]] +deps = ["ArrayInterface", "DocStringExtensions", "IterativeSolvers", "KLU", "Krylov", "KrylovKit", "LinearAlgebra", "RecursiveFactorization", "Reexport", "Requires", "SciMLBase", "Setfield", "SparseArrays", "SuiteSparse", "UnPack"] +git-tree-sha1 = "f27bb8e4eabdb93ed3703c55025b111e045ffe81" +uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +version = "1.12.0" + [[deps.LogExpFunctions]] deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"] git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1" @@ -257,10 +456,10 @@ version = "0.3.6" uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[deps.LoopVectorization]] -deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "Requires", "SLEEFPirates", "Static", "ThreadingUtilities", "UnPack", "VectorizationBase"] -git-tree-sha1 = "c2c1a765d943267ffc01fd6a127fcb482e80f63a" +deps = ["ArrayInterface", "CPUSummary", "ChainRulesCore", "CloseOpenIntervals", "DocStringExtensions", "ForwardDiff", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "SIMDDualNumbers", "SLEEFPirates", "SpecialFunctions", "Static", "ThreadingUtilities", "UnPack", "VectorizationBase"] +git-tree-sha1 = "534aa24fae56f5f0956134d8789ab30d6fe2f615" uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.12.82" +version = "0.12.102" [[deps.MacroTools]] deps = ["Markdown", "Random"] @@ -296,6 +495,18 @@ version = "1.0.2" [[deps.Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" +[[deps.ModiaBase]] +deps = ["DataFrames", "DiffRules", "Measurements", "MonteCarloMeasurements", "OrderedCollections", "StaticArrays", "Unitful"] +path = "D:\\otter\\home\\.julia\\dev\\ModiaBase" +uuid = "ec7bf1ca-419d-4510-bbab-199861c55244" +version = "0.9.2" + +[[deps.ModiaResult]] +deps = ["DataFrames", "Measurements", "MonteCarloMeasurements", "OrderedCollections", "Pkg", "Tables", "Unitful"] +git-tree-sha1 = "8fc8801bf89b9f2efd39952386bda67b3724b0f3" +uuid = "16a87621-1533-42f6-8e19-4a825980cec2" +version = "0.4.1" + [[deps.MonteCarloMeasurements]] deps = ["Distributed", "Distributions", "LinearAlgebra", "MacroTools", "Random", "RecipesBase", "Requires", "SLEEFPirates", "StaticArrays", "Statistics", "StatsBase", "Test"] git-tree-sha1 = "bce0a32d6c64165388eec2573b68000ed06f39c1" @@ -305,14 +516,37 @@ version = "1.0.7" [[deps.MozillaCACerts_jll]] uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +[[deps.MuladdMacro]] +git-tree-sha1 = "c6190f9a7fc5d9d5915ab29f2134421b12d24a68" +uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +version = "0.2.2" + +[[deps.NLSolversBase]] +deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] +git-tree-sha1 = "50310f934e55e5ca3912fb941dec199b49ca9b68" +uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" +version = "7.8.2" + +[[deps.NLsolve]] +deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"] +git-tree-sha1 = "019f12e9a1a7880459d0173c182e6a99365d7ac1" +uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +version = "4.5.1" + [[deps.NaNMath]] -git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5" +git-tree-sha1 = "b086b7ea07f8e38cf122f5016af580881ac914fe" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" -version = "1.0.0" +version = "0.3.7" [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +[[deps.NonlinearSolve]] +deps = ["ArrayInterface", "FiniteDiff", "ForwardDiff", "IterativeSolvers", "LinearAlgebra", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "SciMLBase", "Setfield", "StaticArrays", "UnPack"] +git-tree-sha1 = "b61c51cd5b9d8b197dfcbbf2077a0a4e1505278d" +uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +version = "0.3.14" + [[deps.OffsetArrays]] deps = ["Adapt"] git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed" @@ -333,26 +567,56 @@ git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" version = "0.5.5+0" +[[deps.Optim]] +deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"] +git-tree-sha1 = "bc0a748740e8bc5eeb9ea6031e6f050de1fc0ba2" +uuid = "429524aa-4258-5aef-a3af-852621145aeb" +version = "1.6.2" + [[deps.OrderedCollections]] git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" version = "1.4.1" +[[deps.OrdinaryDiffEq]] +deps = ["Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "ExponentialUtilities", "FastClosures", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "LinearSolve", "Logging", "LoopVectorization", "MacroTools", "MuladdMacro", "NLsolve", "NonlinearSolve", "Polyester", "PreallocationTools", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"] +git-tree-sha1 = "df82fa0f9f90f669cc3cf9e3f0400e431e0704ac" +uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +version = "6.6.6" + [[deps.PDMats]] deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "ee26b350276c51697c9c2d88a072b339f9f03d73" +git-tree-sha1 = "7e2166042d1698b6072352c74cfd1fca2a968253" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.5" +version = "0.11.6" + +[[deps.Parameters]] +deps = ["OrderedCollections", "UnPack"] +git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe" +uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a" +version = "0.12.3" + +[[deps.Parsers]] +deps = ["Dates"] +git-tree-sha1 = "13468f237353112a01b2d6b32f3d0f80219944aa" +uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +version = "2.2.2" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +[[deps.PoissonRandom]] +deps = ["Random", "Statistics", "Test"] +git-tree-sha1 = "44d018211a56626288b5d3f8c6497d28c26dc850" +uuid = "e409e4f3-bfea-5376-8464-e040bb5c01ab" +version = "0.4.0" + [[deps.Polyester]] deps = ["ArrayInterface", "BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "ManualMemory", "PolyesterWeave", "Requires", "Static", "StrideArraysCore", "ThreadingUtilities"] -git-tree-sha1 = "de33c49a06d7eb1eef40b83fe873c1c2cba25623" +git-tree-sha1 = "2232d3865bc9a098e664f69cbe340b960d48217f" uuid = "f517fe37-dbe3-4b94-8317-1923a5111588" -version = "0.6.4" +version = "0.6.6" [[deps.PolyesterWeave]] deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"] @@ -366,11 +630,23 @@ git-tree-sha1 = "db3a23166af8aebf4db5ef87ac5b00d36eb771e2" uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720" version = "1.4.0" +[[deps.PositiveFactorizations]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20" +uuid = "85a6dd25-e78a-55b7-8502-1745935b8125" +version = "0.2.4" + +[[deps.PreallocationTools]] +deps = ["Adapt", "ArrayInterface", "ForwardDiff", "LabelledArrays"] +git-tree-sha1 = "e4cb8d4a2edf9b3804c1fb2c2de57d634ff3f36e" +uuid = "d236fae5-4411-538c-8e31-a6e3d9e00b46" +version = "0.2.3" + [[deps.Preferences]] deps = ["TOML"] -git-tree-sha1 = "2cf929d64681236a2e074ffafb8d568733d2e6af" +git-tree-sha1 = "de893592a221142f3db370f48290e3a2ef39998f" uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.2.3" +version = "1.2.4" [[deps.PrettyTables]] deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"] @@ -396,11 +672,29 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" deps = ["SHA", "Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +[[deps.Random123]] +deps = ["Libdl", "Random", "RandomNumbers"] +git-tree-sha1 = "0e8b146557ad1c6deb1367655e052276690e71a3" +uuid = "74087812-796a-5b5d-8853-05524746bad3" +version = "1.4.2" + +[[deps.RandomNumbers]] +deps = ["Random", "Requires"] +git-tree-sha1 = "043da614cc7e95c703498a491e2c21f58a2b8111" +uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143" +version = "1.5.3" + [[deps.RecipesBase]] git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d" uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" version = "1.2.1" +[[deps.RecursiveArrayTools]] +deps = ["Adapt", "ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"] +git-tree-sha1 = "736699f42935a2b19b37a6c790e2355ca52a12ee" +uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" +version = "2.24.2" + [[deps.RecursiveFactorization]] deps = ["LinearAlgebra", "LoopVectorization", "Polyester", "StrideArraysCore", "TriangularSolve"] git-tree-sha1 = "7ad4c2ef15b7aecd767b3921c0d255d39b3603ea" @@ -418,6 +712,12 @@ git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" uuid = "ae029012-a4dd-5104-9daa-d747884805df" version = "1.3.0" +[[deps.ResettableStacks]] +deps = ["StaticArrays"] +git-tree-sha1 = "256eeeec186fa7f26f2801732774ccf277f05db9" +uuid = "ae5879a3-cd67-5da8-be7f-38c6eb64a37b" +version = "1.1.1" + [[deps.Rmath]] deps = ["Random", "Rmath_jll"] git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f" @@ -433,6 +733,12 @@ version = "0.3.0+0" [[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +[[deps.SIMDDualNumbers]] +deps = ["ForwardDiff", "IfElse", "SLEEFPirates", "VectorizationBase"] +git-tree-sha1 = "62c2da6eb66de8bb88081d20528647140d4daa0e" +uuid = "3cdde19b-5bb0-4aaf-8931-af3e248e098b" +version = "0.1.0" + [[deps.SIMDTypes]] git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c" uuid = "94e857df-77ce-4151-89e5-788b33177be4" @@ -444,13 +750,31 @@ git-tree-sha1 = "61a96d8b89083a53fb2b745f3b59a05359651bbe" uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa" version = "0.6.30" +[[deps.SciMLBase]] +deps = ["ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "RecipesBase", "RecursiveArrayTools", "StaticArrays", "Statistics", "Tables", "TreeViews"] +git-tree-sha1 = "8ff1bf96965b3878ca5d235752ff1daf519e7a26" +uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +version = "1.26.3" + [[deps.Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" +[[deps.Setfield]] +deps = ["ConstructionBase", "Future", "MacroTools", "Requires"] +git-tree-sha1 = "38d88503f695eb0301479bc9b0d4320b378bafe5" +uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" +version = "0.8.2" + [[deps.SharedArrays]] deps = ["Distributed", "Mmap", "Random", "Serialization"] uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" +[[deps.SimpleTraits]] +deps = ["InteractiveUtils", "MacroTools"] +git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231" +uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d" +version = "0.9.4" + [[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" @@ -464,23 +788,29 @@ version = "1.0.1" deps = ["LinearAlgebra", "Random"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +[[deps.SparseDiffTools]] +deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "Graphs", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays", "VertexSafeGraphs"] +git-tree-sha1 = "87efd1676d87706f4079e8e717a7a5f02b6ea1ad" +uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804" +version = "1.20.2" + [[deps.SpecialFunctions]] deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "85e5b185ed647b8ee89aa25a7788a2b43aa8a74f" +git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.1.3" +version = "2.1.4" [[deps.Static]] deps = ["IfElse"] -git-tree-sha1 = "a8f30abc7c64a39d389680b74e749cf33f872a70" +git-tree-sha1 = "7f5a513baec6f122401abfc8e9c074fdac54f6c1" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.3.3" +version = "0.4.1" [[deps.StaticArrays]] deps = ["LinearAlgebra", "Random", "Statistics"] -git-tree-sha1 = "6354dfaf95d398a1a70e0b28238321d5d17b2530" +git-tree-sha1 = "74fb527333e72ada2dd9ef77d98e4991fb185f04" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.4.0" +version = "1.4.1" [[deps.Statistics]] deps = ["LinearAlgebra", "SparseArrays"] @@ -499,21 +829,49 @@ uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" version = "0.33.16" [[deps.StatsFuns]] -deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] -git-tree-sha1 = "f35e1879a71cca95f4826a14cdbf0b9e253ed918" +deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] +git-tree-sha1 = "25405d7016a47cf2bd6cd91e66f4de437fd54a07" uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" -version = "0.9.15" +version = "0.9.16" + +[[deps.SteadyStateDiffEq]] +deps = ["DiffEqBase", "DiffEqCallbacks", "LinearAlgebra", "NLsolve", "Reexport", "SciMLBase"] +git-tree-sha1 = "3e057e1f9f12d18cac32011aed9e61eef6c1c0ce" +uuid = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" +version = "1.6.6" + +[[deps.StochasticDiffEq]] +deps = ["Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqJump", "DiffEqNoiseProcess", "DocStringExtensions", "FillArrays", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "Logging", "MuladdMacro", "NLsolve", "OrdinaryDiffEq", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"] +git-tree-sha1 = "5f88440e7470baad99f559eed674a46d2b6b96f7" +uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" +version = "6.44.0" [[deps.StrideArraysCore]] deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "ManualMemory", "Requires", "SIMDTypes", "Static", "ThreadingUtilities"] -git-tree-sha1 = "2025a5d6564a93fa5b499dd216c0bc44537fb0d4" +git-tree-sha1 = "e0a02838565c4600ecd1d8874db8cfe263aaa6c7" uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da" -version = "0.2.11" +version = "0.2.12" [[deps.SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" +[[deps.SuiteSparse_jll]] +deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"] +uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" + +[[deps.Sundials]] +deps = ["CEnum", "DataStructures", "DiffEqBase", "Libdl", "LinearAlgebra", "Logging", "Reexport", "SparseArrays", "Sundials_jll"] +git-tree-sha1 = "76d881c22a2f3f879ad74b5a9018c609969149ab" +uuid = "c3572dad-4567-51f8-b174-8c6c989267f4" +version = "4.9.2" + +[[deps.Sundials_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "OpenBLAS_jll", "Pkg", "SuiteSparse_jll"] +git-tree-sha1 = "04777432d74ec5bc91ca047c9e0e0fd7f81acdb6" +uuid = "fb77eaff-e24c-56d4-86b1-d163f2edb164" +version = "5.2.1+0" + [[deps.TOML]] deps = ["Dates"] uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" @@ -540,9 +898,9 @@ uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [[deps.ThreadingUtilities]] deps = ["ManualMemory"] -git-tree-sha1 = "884539ba8c4584a3a8173cb4ee7b61049955b79c" +git-tree-sha1 = "f8629df51cab659d70d2e5618a430b4d3f37f2c3" uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5" -version = "0.4.7" +version = "0.5.0" [[deps.TimerOutputs]] deps = ["ExprTools", "Printf"] @@ -550,11 +908,17 @@ git-tree-sha1 = "97e999be94a7147d0609d0b9fc9feca4bf24d76b" uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" version = "0.5.15" +[[deps.TreeViews]] +deps = ["Test"] +git-tree-sha1 = "8d0d7a3fe2f30d6a7f833a5f19f7c7a5b396eae6" +uuid = "a2a6695c-b41b-5b7d-aed9-dbfdeacea5d7" +version = "0.3.0" + [[deps.TriangularSolve]] deps = ["CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "LoopVectorization", "Polyester", "Static", "VectorizationBase"] -git-tree-sha1 = "c3ab8b77b82fd92e2b6eea8a275a794d5a6e4011" +git-tree-sha1 = "5cbc1a4551fcf8afe8f80bb4f1f13e3271ee2656" uuid = "d5829a12-d9aa-46ab-831f-fb7c9ab06edf" -version = "0.1.9" +version = "0.1.10" [[deps.UUIDs]] deps = ["Random", "SHA"] @@ -580,10 +944,22 @@ git-tree-sha1 = "e9a35d501b24c127af57ca5228bcfb806eda7507" uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f" version = "0.21.24" +[[deps.VertexSafeGraphs]] +deps = ["Graphs"] +git-tree-sha1 = "8351f8d73d7e880bfc042a8b6922684ebeafb35c" +uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" +version = "0.2.0" + [[deps.Zlib_jll]] deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +[[deps.ZygoteRules]] +deps = ["MacroTools"] +git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0" +uuid = "700de1a5-db45-46bc-99cf-38207098b444" +version = "0.2.2" + [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl", "OpenBLAS_jll"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" diff --git a/Project.toml b/Project.toml index 1639114..91d0559 100644 --- a/Project.toml +++ b/Project.toml @@ -1,18 +1,15 @@ name = "ModiaBase" uuid = "ec7bf1ca-419d-4510-bbab-199861c55244" authors = ["Hilding Elmqvist ", "Martin Otter "] -version = "0.9.2" +version = "0.10.0" [deps] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] @@ -21,8 +18,6 @@ DiffRules = "1" Measurements = "2" MonteCarloMeasurements = "1" OrderedCollections = "1" -RecursiveFactorization = "0.2" StaticArrays = "1" -TimerOutputs = "0.5" Unitful = "1" julia = "1.7" diff --git a/docs/src/index.md b/docs/src/index.md index 7e07341..7eef6b6 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -77,6 +77,16 @@ julia> ]add ModiaBase ## Release Notes +### Version 0.10.0 + +**Non-backwards** compatible changes + +- EquationAndStateInfo.jl and StateSelection.jl moved to Modia (ModiaLang is merged into Modia), because + the AST generation in these files depends on details of CodeGeneration.jl of Modia/ModiaLang. + +- TestLinearEquations.jl also moved to Modia/ModiaLang. + + ### Version 0.9.2 - Minor (efficiency) improvement of linear equation system if iteration variables are SVectors. diff --git a/src/EquationAndStateInfo.jl b/src/EquationAndStateInfo.jl deleted file mode 100644 index 8380f61..0000000 --- a/src/EquationAndStateInfo.jl +++ /dev/null @@ -1,1026 +0,0 @@ -# License for this file: MIT (expat) -# Copyright 2020-2021, DLR Institute of System Dynamics and Control -# Author: Martin Otter, DLR-SR -# -# Provide information about the structure of the ODE or DAE equations. - -import DataFrames -import OrderedCollections -using LinearAlgebra -using StaticArrays -import TimerOutputs -import RecursiveFactorization - - -export StateCategory, XD, XALG, XLAMBDA, XMUE -export ResidualCategory, FD, FC_ALG, FC_LOW_HIGH, FC_MUE - -export LinearEquations - -export EquationInfo, StateElementInfo, update_equationInfo!, get_stateNames, get_x_table -export initialStateVector, updateEquationInfo! - -export EquationInfoStatus, MANUAL, CODE_GENERATION, SOLVER_MODEL - - - -""" - isFixedLengthStartOrInit(startOrInit::Any, name::AbstractString) - -Returns true, if `startOrInit` (so start value, init value, `nothing`) characterizes a value -with fixed length, so `length(v_startOrInit)` is fixed after compilation (= either a number or a StaticArray). -Note, if a start-value is not defined (`startOrInit = nothing`) a scalar default is used. - -An error is triggered, if `startOrInit` is neither a `Number` nor an `AbstractVector` nor `nothing`. -In this case argument `name` is used in the error message (currently, EquationAndStateInfo is restricted -to numbers and vectors. In a subsequent version, two and more-dimensional arrays might be supported). - -# Examples -``` -isFixedLengthStartOrInit(1.0) # = true -isFixedLengthStartOrInit(SVector{3}(1,2,3)) # = true -isFixedLengthStartOrInit([1,2,3]) # = false (size can be changed after compilation) -isFixedLengthStartOrInit(nothing) # = true (is assumed to be a scalar with value zero) -isFixedLengthStartOrInit(missing) # error -``` -""" -isFixedLengthStartOrInit(startOrInit, name::AbstractString) = begin - if !(startOrInit isa Nothing || startOrInit isa Number || startOrInit isa AbstractVector) - error("Start or init value $startOrInit of variable $name\nis neither a subtype of Number, subtype of AbstractVector, nor is nothing") - end - !(startOrInit isa AbstractVector && !(startOrInit isa StaticArray)) -end - - -""" - @enum StateCategory - -Type of a state. - -| value | description | -|:----------|:----------------------------------------------------------------------------------| -| `XD` | Differential variable (derivative of variable appears in equations) | -| `XALG` | Algebraic variable that is included in x | -| `XLAMBDA` | Algebraic variable that is included in der_x | -| `XMUE` | Algebraic variable used for stabilization and included in der_x (exact value = 0) | -""" -@enum StateCategory XD=1 XALG XLAMBDA XMUE - - -""" - @enum ResidualCategory - -Type of a residual (see also [`StateCategory`](@ref) for the type of a state). - -| value | description | -|:--------------|:-------------------------------------------------------------------------| -| `FD` | Differential equation (depends on XD, XLAMBDA, XMUE variables) | -| `FC_ALG` | Constraint equation that is purely algebraic (depends on XALG variables) | -| `FC_LOW_HIGH` | Constraint equation on lowest and on highest derivative level | -| `FC_MUE` | Constraint equation that is associated with a stabilizer XMUE | -""" -@enum ResidualCategory FD=1 FC_ALG FC_LOW_HIGH FC_MUE - - - -const niter_max = 20 # Maximum number of fixed-point iterations to solve A*x = b - -""" - leq = LinearEquations{FloatType}(x_names::Vector{String}, x_vec_julia_names::AbstractVector, - x_lengths::Vector{Int}, nx_fixedLength::Int, A_is_constant::Bool; - useRecursiveFactorizationUptoSize = 0) - -Define linear equation system "A*x=b" with `x::Vector{FloatType}`. - -- `x` is constructed from a set of scalar or vector variables i with names `x_names[i]` - and length `x_length[i]` (that is `length(x) = sum(x_lengths)`). - x_names[1:nx_fixedLength] are elements with fixed lengths (dimensions do not change after compilation). - x_names[nx_fixedLength+1:end] are vector-valued elements where the dimensions may changer after compilation. - -- `x_vec_julia_names` are the Julia names of the vector-valued elements. - -- If `A_is_constant = true` then `A` is a matrix that is constant after initialization - -- If length(x) <= useRecursiveFactorizationUptoSize, then linear equation systems will be solved with - `RecursiveFactorization.jl` instead of the default `lu!(..)` and `ldiv!(..)`. - -For details how to use this constructor, see [`LinearEquationsIteration!`](@ref). -""" -mutable struct LinearEquations{FloatType <: Real} - odeMode::Bool # Set from the calling function after LinearEquations was instantiated (default: true) - # = true (standard mode): Compute "x" from equation "residuals = A*x - b" - # = false (DAE solver): During events (including initialization) - # compute "x" as for odeMode=true. Outside of events: - # (1) "x" is set from the outside (= der(x_dae) provided by DAE solver) - # (2) Compute "residuals" from "residuals := A*x - b" - # (3) From the outside copy "residuals" into "residuals_dae" of the DAE solver. - - A_is_constant::Bool # = true, if A-matrix is constant - x_names::Vector{String} # Names of the x-variables - x_vec_julia_names::AbstractVector # Julia names of the vector-valued x-variables - # (needed to generate code for - # if _m.storeResult - # x_vec_julia_names[1] = deepcopy(x_vec_julia_names[1]) - # ...) - x_lengths::Vector{Int} # Lengths of the x-variables (sum(x_lengths) = length(x)) - nx_fixedLength::Int # x_names[1:nx_fixedLength] are elements with fixed lengths (after compilation) - # x_names[nx_fixedLength+1:end] are vector elements where the size may change after compilation - x_vec::Vector{Vector{FloatType}} # x_vec[1] is the value vector of x_names[nx_fixedLength+1] - # x_vec[2] is the value vector of x_names[nx_fixedLength+2] - # <..> - - - A::Matrix{FloatType} - b::Vector{FloatType} - x::Vector{FloatType} # Values of iteration variables - pivots::Vector{Int} # Pivot vector if recursiveFactorization = true - residuals::Vector{FloatType} # Values of the residuals FloatType vector; length(residuals) = length(x) - - # Iteration status of for-loop - mode::Int # Operational mode (see function LinearEquationsIteration!) - niter::Int # Current number of iterations in the fix point iteration scheme - niter_max::Int # Maximum number of iterations - success::Bool # = true, if solution of A*x = b is successfully computed - # = false, if solution is not computed; continue with fix point iteration - - # For warning message if niter > niter_max - inconsistentPositive::Vector{String} - inconsistentNegative::Vector{String} - - # Constructed during initialization - useRecursiveFactorizationUptoSize::Int - useRecursiveFactorization::Bool # = true, if RecursiveFactorization.jl shall be used to solve the linear equation system - - luA::LU{FloatType,Array{FloatType,2}} # lu-Decomposition of A - - function LinearEquations{FloatType}(x_names::Vector{String}, x_vec_julia_names::AbstractVector, - x_lengths::Vector{Int}, nx_fixedLength::Int, A_is_constant::Bool; - useRecursiveFactorizationUptoSize::Int = 0) where {FloatType <: Real} - @assert(length(x_names) > 0) - @assert(length(x_names) == length(x_lengths)) - nx = sum(x_lengths) - @assert(nx > 0) - useRecursiveFactorization = nx <= useRecursiveFactorizationUptoSize - - # Allocate value storage for vector elements - nx_vec = length(x_names) - nx_fixedLength - x_vec = fill(FloatType[], nx_vec) - j = nx_fixedLength - for i = 1:nx_vec - j += 1 - x_vec[i] = zeros(FloatType, x_lengths[j]) - end - - new(true, A_is_constant, x_names, Any[], x_lengths, nx_fixedLength, x_vec, - zeros(FloatType,nx,nx), zeros(FloatType,nx), zeros(FloatType,nx), fill(0,nx), zeros(FloatType,nx), - -2, 0, niter_max, false, String[], String[], - useRecursiveFactorizationUptoSize, useRecursiveFactorization) - end -end -LinearEquations(args...) = LinearEquations{Float64}(args...) - - -""" - b = copy_x_into_x_vec!(leq) - -If vector valued elements of x-vector, copy vector valued elements of leq.x -into vector leq.x_vec. The function returns true -""" -function copy_x_into_x_vec!(leq::LinearEquations{FloatType})::Bool where {FloatType} - i_vec = 0 - i_x = leq.nx_fixedLength - for i = leq.nx_fixedLength+1:length(leq.x_names) - i_vec += 1 - x_vec = leq.x_vec[i_vec] - for j = 1:leq.x_lengths[i] - x_vec[j] = leq.x[i_x+j] - end - i_x += leq.x_lengths[i] - end - return true -end - - - -""" - iterating = LinearEquationsIteration!(leq, isInitial, [solve, isStoreResult,] time, timer) - -This function solves a linear equation system in residual form "residual = A*x - b" -by iterating with a while loop over this system (arguments `solve, isStoreResult` are optional -and have only an effect if leq.odeMode=false, that is if the linear equation system is -solved in DAE mode): - -```julia -function getDerivatives!(_der_x, _x, _m, _time)::Nothing - _leq::Union{Nothing,LinearEquations{FloatType}} = nothing - ... - _leq = _m.linearEquations[] # leq::LinearEquations{FloatType} - _leq.mode = -3 # initializes the iteration - while LinearEquationsIteration!(_leq, _m.isInitial, _m.solve_leq, _m.storeResult, _m.time, _m.timer) - x1 = _leq.x[1] - x2 = _leq.x[2] - x3 = _leq.x_vec[1] - ... - v_solved1 = f(x1, x2, ..., positive(x1,.., _leq)) - v_solved2 = f(x1, x2, ..., v_solved1) - ... - appendResidual!(_leq.residuals, < getResiduals(x1,x2,...) > ) - end - ... -end -``` - -Note, A and b can be functions of linear event operators, such as `positive(c_i'*x - d_i)`. -In this case this system is solved by a fixed point iteration scheme, -that is `A*x = b` is solved until potentially present -`positive(c_i'*x - d_i)` calls are consistent to x -(`c_i'*x - d_i` and positive(..) must have the same sign). -Iteration is terminated after 20 iterations if no convergence is found. -In this case a warning message is triggered and simulation is continued -(it might be that simulation is still successful, even if `x` and -`positive(..)` are temporarily not consistent to each other). - -The current values of A,x,b,residuals are stored in leq. -If A is fixed after initialization (leq.A_is_constant = true), then A is computed -only once at initialization, the LU decomposition of A is stored in leq -and used in subsequent calls to solve the equation system. - - -# Input arguments - -- `leq::LinearEquations{FloatType}`: Instance of `LinearEquations`. -- `isInitial::Bool`: = true: Called during initialization. -- `solve::Bool`: = true: leq.x is computed by LinearEquationsIteration!. - = false: leq.x has been set by calling environment - (typically when called from DAE integrator). - Note, at events and at terminate, solve must be true). -- `isStoreResult::Bool`: = true: Called at a communication point (store results). -- `time`: Actual time instant (only used for warning/error messages). -- `timer::TimerOutputs.TimerOutput`: Timer used to measure the solution of the linear equation system. - - -# Output arguments - -- iterating::Bool: = true : Continue iteration of while-loop. - = false: Immediately break while-loop after this function was terminated. - - -# Enhancing efficiency - -Variable `leq.mode::Int` can be accessed (read-only) in the body of the while-loop to enhance efficiency. -For example, if residuals are computed in a function and the function evaluation is expensive. - -``` -leq.mode = -2: Residuals might be computed, but they are not used. -leq.mode = -1: Compute "residuals .= A*x - b" -leq.mode = 0: Compute "residuals .= A*0 - b" -leq.mode > 0: Compute "residuals .= A*e_i - b" # e_i = unit vector with i = leq.mode -``` - - -# Hidden argument `leq.mode::Int` on input - -```julia -leq.mode = -3 # LinearEquationsIteration! is called the first time - if leq.odeMode || solve - # ODE mode or DAE mode at an event (solve "x" from equation "residuals = A*x - b") - # Initialize fixed point iteration or continue fixed point iteration (if in DAE mode) - leq.niter = 0 # Number of event iterations - leq.success = false # Event iteration was not yet successful - leq.x .= 0 - leq.mode = 0 # Compute "residuals .= A*0 - b" - else # solve = false - # DAE mode (but not at an event) - leq.niter = 0 - leq.success = true - leq.x_is_zero = false - if isStoreResult - leq.mode = -2 # Residuals might be computed, but they are not used. - else - leq.mode = -1 # Compute "residuals .= A*x - b" - end - end - return true # Continue while-loop - -leq.mode = -2 # Terminate while-loop or initialize next event iteration - if leq.success # Set from leq.mode=-2 or leq.mode>=0 or from event operators, such as positive!(...) - return false # Terminate while-loop - elseif leq.niter > 20 - # Event iteration failed (variables are not consistent) - return false # Terminate while-loop - else - # Initialize next event iteration - leq.x .= 0 - leq.x_is_zero = true - leq.mode = 0 # Compute "residuals .= A*0 - b" - return true # Continue while-loop - end - -leq.mode = -1: @assert(!leq.odeMode && !solve) - # DAE mode, but not at an event - return false # Terminate while-loop - -leq.mode = 0: @assert(leq.odeMode || solve) - # ODE mode or DAE mode at an event (solve "x" from equation "residuals = A*x - b") - leq.b = -leq.residuals - leq.x[1] = 1.0 - leq.mode = 1 # Compute "residuals := A*e_1 - b" - return true # Continue while-loop - -leq.mode > 0: @assert(leq.odeMode || solve) - j = leq.mode - leq.A[:,j] = leq.residuals + leq.b - leq.x[j] = 0.0 - if j < nx - leq.x[j+1] = 1.0 - leq.mode = j+1 # Compute A[:,j+1] after the next iteration - return true # Continue while-loop - elseif j == nx - leq.x = solve(leq.A,leq.b) # Solve linear equation system A*x = b for x. - leq.mode = -1 # Compute all variables IN the next iteration as function of x - leq.niter += 1 # Increment number of iterations to solve A*x = b - leq.success = true # Terminate for-loop at the beginning of the next iteration, - # provided no positive(..) call changes its value. - # (leq.success is set to false in positive(..), if the return value changes). - return true # Continue while-loop - end -``` -""" -LinearEquationsIteration!(leq::LinearEquations, isInitial, time, timer) = LinearEquationsIteration!(leq, isInitial, true, false, time, timer) -function LinearEquationsIteration!(leq::LinearEquations{FloatType}, isInitial::Bool, solve::Bool, - isStoreResult::Bool, time, timer)::Bool where {FloatType} - mode = leq.mode - nx = length(leq.x) - - if mode == -3 - # LinearEquationsIteration! is called the first time in the current model evaluation - leq.niter = 0 # Number of event iterations - empty!(leq.inconsistentPositive) - empty!(leq.inconsistentNegative) - - if leq.odeMode || solve - # ODE mode or DAE mode at an event (solve "x" from equation "residuals = A*x - b") - # Initialize fixed point iteration or continue fixed point iteration (if in DAE mode) - leq.success = false # Event iteration was not yet successful - leq.x .= 0 - leq.mode = 0 # Initialize leq and compute "residuals .= A*0 - b" - else - # DAE mode (but not at an event) - leq.success = true - if isStoreResult - leq.mode = -2 # Residuals might be computed, but they are not used. - else - leq.mode = -1 # Initialize leq and compute "residuals .= A*x - b" - end - end - empty!(leq.residuals) - return copy_x_into_x_vec!(leq) # Continue while-loop - - elseif mode == -2 - # Terminate while-loop or initialize next event iteration - if leq.success # Set from leq.mode=-2 or leq.mode>=0 or from event operators, such as positive!(...) - return false # Terminate while-loop - elseif leq.niter > niter_max - str = "" - if length(leq.inconsistentPositive) > 0 - str = str * "positive(expr) is inconsistent for expr = $(leq.inconsistentPositive)." - end - if length(leq.inconsistentNegative) > 0 - if length(str) > 0 - str = str * "\n" - end - str = str * "negative(expr) is inconsistent for expr = $(leq.inconsistentNegative)." - end - @warn "At time = $time, no consistent solution found for mixed linear equation system.\n" * - "Simulation is continued although some variables might not be correct at this time instant.\n$str" - return false # Terminate while-loop - end - leq.x .= 0 - leq.mode = 0 # Compute "residuals .= A*0 - b" - empty!(leq.residuals) - return copy_x_into_x_vec!(leq) # Continue while-loop - - elseif mode < -3 || mode > nx - @goto ERROR - end - - x = leq.x - A = leq.A - b = leq.b - residuals = leq.residuals - - if length(residuals) != length(x) - error("Function LinearEquationsIteration! wrongly used:\n", - "length(leq.residuals) = ", length(leq.residuals), ", length(leq.x) = ", length(leq.x)) - end - - if mode == -1 - @assert(!leq.odeMode && !solve) - return false # Terminate while-loop (leq.residuals must be copied into DAE residuals) - - elseif !leq.A_is_constant || isInitial # A is not constant or A is constant and isInitial = true - if mode == 0 - # residuals = A*x - b -> b = -residuals) - for i = 1:nx - b[i] = -residuals[i] - end - leq.mode = 1 - x[1] = max(FloatType(1e-3), abs(b[1])) #convert(FloatType, 1) - empty!(leq.residuals) - return copy_x_into_x_vec!(leq) - end - - # residuals = A*e_j*x_j - b -> A[:,j] = (residuals + b)/x[j] - j = mode - for i = 1:nx - A[i,j] = (residuals[i] + b[i])/x[j] - end - x[j] = 0 - - if j < nx - leq.mode += 1 - x[leq.mode] = max(FloatType(1e-3), abs(b[leq.mode])) # convert(FloatType, 1) - empty!(leq.residuals) - return copy_x_into_x_vec!(leq) - end - - # Solve linear equation system - if nx == 1 - x[1] = b[1]/A[1,1] - if !isfinite(x[1]) - error("Linear scalar equation system is singular resulting in: ", leq.x_names[1], " = ", x[1]) - end - else - x .= b - if leq.useRecursiveFactorization - ModiaBase.TimerOutputs.@timeit timer "ModiaBase LinearEquationsIteration! (solve A*x=b Rec.Fac.)" begin - leq.luA = RecursiveFactorization.lu!(A, leq.pivots) - ldiv!(leq.luA, x) - end - else - ModiaBase.TimerOutputs.@timeit timer "ModiaBase LinearEquationsIteration! (solve A*x=b)" begin - leq.luA = lu!(A) - ldiv!(leq.luA, x) - end - end - end - - elseif leq.A_is_constant && !isInitial # isInitial=false, LU decomposition of A is available in leq.luA - for i = 1:nx - x[i] = -residuals[i] - end - - # Solve linear equation system - if nx == 1 - x[1] = x[1]/A[1,1] - if !isfinite(x[1]) - error("Linear scalar equation system is singular resulting in: ", leq.x_names[1], " = ", x[1]) - end - else - ldiv!(leq.luA, x) - end - - else - @goto ERROR - end - - # Linear equation system solved - leq.mode = -2 # Compute all variables IN the next iteration as function of x - leq.niter += 1 # Increment number of iterations to solve A*x = b - leq.success = true # Terminate for-loop at the beginning of the next iteration, - # provided no positive(..) call changes its value. - # (leq.success is set to false in positive(..), if the return value changes). - - empty!(leq.residuals) - return copy_x_into_x_vec!(leq) - - @label ERROR - error("Should not occur (Bug in file ModiaBase/src/EquationAndStateInfo.jl):\n", - " time = $time\n", - " isInitial = $isInitial,\n", - " solve = $solve,\n", - " isStoreResult = $isStoreREsult,\n", - " leq.odeMode = $(leq.odeMode),\n", - " leq.mode = $(leq.mode),\n", - " leq.x_names = $(leq.x_names),\n", - " leq.x_lengths = $(leq.x_lengths),\n") -end - - - -""" - @enum EquationInfoStatus - -Status of an EquationInfo instance: - -- `MANUAL`: Is defined manually. The following variables in x_info::Vector{StateElementInfo} - are **not** defined: - x_names_julia, der_x_names_julia, length, unit, startIndex. - Also variables nx and x_infoByIndex are not defined. - With function [`update_equationInfo!`](@ref), the variables length, unit, startIndex - in x_info::Vector{StateElementInfo} are computed, as well as nx and x_infoByIndex. - -- `CODE_GENERATION`: Is constructed during code generation with getSortedAndSolvedAST(..). - The following variables in x_info::Vector{StateElementInfo} - are **not** defined: startIndex. - Also variables nx and x_infoByIndex are not defined. - With function [`update_equationInfo!`](@ref), the variables startIndex - in x_info::Vector{StateElementInfo} are computed, as well as nx and x_infoByIndex. - -- `SOLVER_MODEL`: Is used during simulation in a SolverModel. With function - [`update_equationInfo!`](@ref), missing variables are constructed depending - on the information given with `MANUAL` or `CODE_GENERATION` and the actual - modelValues instance (here unit and length information is available). -""" -@enum EquationInfoStatus MANUAL=1 CODE_GENERATION SOLVER_MODEL - - - -""" - xe_info = StateElementInfo(...) - -Return an instance of the mutable struct `StateElementInfo` that defines the information -for one element of the state vector. There are three constructors: - -- Default constructor (all variables under section Arguments are given; - used to write/read the complete information). - -- All variables under section Arguments are given with exception of startIndex - (used for `EquationInfoStatus = CODE_GENERATION`). - -- StateElementInfo(x_name, der_x_name, stateCategory; fixed=nothing, nominal=NaN, unbounded=false) - (used for `EquationInfoStatus = MANUAL`): - -# Arguments -- x_name: Name of x-element or "" if no name (if stateCatebory = XLAMBDA or XMUE) -- x_name_julia: Julia name of x-element in getDerivatives!/getResiduals! function - or "" if not needed (since no code generation). -- der_x_name: Name of der_x-element or "" if either `der(x_name)` or if no name - (if stateCategory = XALG). -- der_x_name_julia: Julia name of der_x-element in getDerivatives!/getResiduals! function - or "" if not needed (since no code generation). -- stateCategory::StateCategory: Category of the state -- length: length of x-element (or -1 if not yet known) -- unit: unit of XD, XALG (x_name) or XLAMBDA, XMUE (der_x_name) or "" if not yet known) -- fixed: false (= guess value) or true (= not changed by initialization). - Only relevant for ode=false, otherwise ignored. -- nominal: Nominal value (NaN if determined via start value) -- unbounded: false or true -- startIndex: start index of x-element with respect to x-vector - or -1 if not yet known. -""" -mutable struct StateElementInfo - x_name::String # Modia name of x-element or "" if no name (for ) - x_name_julia # Julia name of x-element in getDerivatives! function - # or "", if not needed (since no code generation). - der_x_name::String # Modia name of der_x-element or "" if either "der(x_name)" or if no name, - der_x_name_julia # Julia name of der_x-element in getDerivatives! function - # or not needed (since no code generation) - stateCategory::StateCategory # category of the state - unit::String # unit of x-element as string (or "" if not yet known) - startOrInit::Any # start or init value or nothing depending on fixed - # (init : if value and fixed=true - # start : if value and fixed=false - # nothing: neither start nor init value; fixed=false) - fixed::Bool # - nominal::Float64 # nominal value (NaN if determined via start value) - unbounded::Bool # false or true - fixedLength::Bool # = true, if length of value is fixed at compile time (= scalar or StaticArray) - # = false, if length of value is determined before initialization - scalar::Bool # = true, if scalar - # = false, if vector - length::Int # length of x-element - startIndex::Int # start index of element with respect to x-vector (if fixedLength=true) - # or with respect to x_vec-vector (if fixedLength=false) -end - - - -# Constructor for code-generation -StateElementInfo(x_name, x_name_julia, der_x_name, der_x_name_julia, - stateCategory, unit, startOrInit, fixed, nominal, unbounded) = StateElementInfo( - x_name, x_name_julia, der_x_name, der_x_name_julia, - stateCategory, unit, startOrInit, fixed, nominal, unbounded, - isFixedLengthStartOrInit(startOrInit, x_name), !(startOrInit isa AbstractArray), - startOrInit isa Nothing ? 1 : length(startOrInit), -1) - - -# Constructor for reading StateElementInfo after code-generation (x_name_julia and der_x_name_julia are not included -#StateElementInfo(x_name, der_x_name, -# stateCategory, unit, startOrInit, fixed, nominal, unbounded) = StateElementInfo( -# x_name, :(), der_x_name, :(), -# stateCategory, unit, startOrInit, fixed, nominal, unbounded, -1) - -# Constructor for manually generated equation info. -#StateElementInfo(x_name, der_x_name, stateCategory=XD; fixed=false, nominal=NaN, unbounded=false) = -# StateElementInfo(x_name, :(), der_x_name, :(), stateCategory, "", nothing, fixed, nominal, unbounded, -1) - - -function Base.show(io::IO, xe_info::StateElementInfo) - print(io, "ModiaBase.StateElementInfo(") - show( io, xe_info.x_name) - #print(io, ",") - #show( io, xe_info.x_name_julia) - print(io, ",") - show( io, xe_info.der_x_name) - #print(io, ",") - #show( io, xe_info.der_x_name_julia) - print(io, ",ModiaBase.", xe_info.stateCategory) - print(io, ",") - show( io, xe_info.unit) - print(io, ",", xe_info.startOrInit) - print(io, ",", xe_info.fixed) - print(io, ",", xe_info.nominal) - print(io, ",", xe_info.unbounded) - print(io, ",", xe_info.fixedLength) - print(io, ",", xe_info.scalar) - print(io, ",", xe_info.length) - print(io, ",", xe_info.startIndex) - print(io, ")") - return nothing -end - - -""" - x_table = get_x_table(x_info::Vector{StateElementInfo}) - -Return the state element info as DataFrames.DataFrame table, for example -to print it as: - -``` -show(x_table, allrows=true, allcols=true, summary=false, eltypes=false) -``` -""" -function get_x_table(x_info::Vector{StateElementInfo}) - x_table = DataFrames.DataFrame(name=String[], unit=String[], startOrInit=[], fixed=Bool[], nominal=Float64[], unbounded=Bool[]) - - for xe_info in x_info - push!(x_table, (xe_info.x_name, xe_info.unit, xe_info.startOrInit, xe_info.fixed, xe_info.nominal, xe_info.unbounded)) - end - - return x_table -end - - -""" - nx = stateVectorLength(x_info::Vector{StateElementInfo}) - -Return the length of the state vector -""" -function stateVectorLength(x_info::Vector{StateElementInfo})::Int - nx = 0 - for xe_info in x_info - nx += xe_info.length - end - return nx -end - - - -""" - eqInfo = EquationInfo(; - status = MANUAL, - ode = true, - nz = 0, - x_info = StateElementInfo[], - nx = -1, - nxFixedLength = nx, - residualCategories = ResidualCategory[], - linearEquations = Tuple{Vector{String},AbstractVector,Vector{Int},Int,Bool}[], - vSolvedWithFixedTrue = String[], - defaultParameterAndStartValues = nothing, - ResultType = nothing, - ResultTypeHasFloatType = false) - -Return instance `eqInfo` that defines the information for the equation system. - -# Arguments -- status::EquationInfoStatus: Defines the variables that have a value. -- ode: = true if ODE-interface (`getDerivatives!`), - = false if DAE-interface (`getResiduals!`). -- nz: Number of zero crossing functions. -- x_info: Vector of StateElementInfo elements provding info for every x-element -- residualCategories: If ode=true, length(residualCategories) = 0. - If ode=false: residualCategories[i] is the `ResidualCategory`](@ref) of x-element "i". -- `linearEquations::Vector{Tuple{Vector{String},AbstractVector,Vector{Int},Int,Bool}}`: - linearEquations[i] defines a - ModiaBase.LinearEquations system, where the first tuple value - is a vector of the names of the unknowns, - the second tuple value is a vector of the Julia names of the vector-valued elements, - the third tuple value is a vector with the lengths of the unknowns, - the fourth tuple value is the number of residuals and the fifth tuple value - defines whether the coefficient matrix A - has only constant entries (=true) or not (=false). -- `vSolvedWithFixedTrue::Vector{String}`: Vector of variables that are computed - from other variables and have `fixed=true`. During initialization - it is checked whether the calculated values and the start values of - these variables agree. If this is not the case, an error is triggered. -- `defaultParameterAndStartValues`: Dictionary of default paramter and default start values. -- `ResultType::AbstractModelValues`: ModelValues type that shall be used for the result generation - (= ModelValues struct without parameters). -- `ResultTypeHasFloatType::Bool`: = false, if `ResultType` is not parameterized. - = true, if `ResultType` has parameter `FloatType`. -""" -mutable struct EquationInfo - status::EquationInfoStatus - ode::Bool # = true if ODE model interface, otherwise DAE model interface - nz::Int # Number of crossing functions - x_info::Vector{StateElementInfo} - residualCategories::Vector{ResidualCategory} # If ode=true, length(residualCategories) = 0 - # If ode=false, residualCategories[j] is the ResidualCategory of residual[j] - linearEquations::Vector{Tuple{Vector{String},AbstractVector,Vector{Int},Int,Bool}} - vSolvedWithFixedTrue::Vector{String} - nx::Int # = length(x) or -1 if not yet known - nxFixedLength::Int # x_info[1:nxFixedLengt] are states with fixed length (does not change after compilation) - x_infoByIndex::Vector{Int} # i = x_infoByIndex[j] -> x_info[i] - # or empty vector, if not yet known. - x_dict::OrderedCollections.OrderedDict{String,Int} # x_dict[x_name] returns the index of x_name with respect to x_info - der_x_dict::OrderedCollections.OrderedDict{String,Int} # der_x_dict[der_x_name] returns the index of der_x_name with respect to x_info - defaultParameterAndStartValues::Union{AbstractDict,Nothing} - ResultType - ResultTypeHasFloatType::Bool -end - - -EquationInfo(; status = MANUAL, - ode = true, - nz = 0, - x_info = StateElementInfo[], - residualCategories = ResidualCategory[], - linearEquations = Tuple{Vector{String},AbstractVector,Vector{Int},Int,Bool}[], - vSolvedWithFixedTrue = String[], - nxFixedLength = -1, - x_infoByIndex = Int[], - defaultParameterAndStartValues::Union{AbstractDict,Nothing} = nothing, - ResultType = nothing, - ResultTypeHasFloatType = false) = EquationInfo(status, ode, nz, x_info, - residualCategories, linearEquations, - vSolvedWithFixedTrue, stateVectorLength(x_info), -1, x_infoByIndex, - OrderedCollections.OrderedDict{String,Int}(), - OrderedCollections.OrderedDict{String,Int}(), - defaultParameterAndStartValues, - ResultType, ResultTypeHasFloatType) - - -""" - initEquationInfo!(eqInfo::EquationInfo) - -Set eqInfo.x_dict, eqInfo.der_x_dict, eqInfo.nx and eqInfo.x_info[:].startIndex -""" -function initEquationInfo!(eqInfo::EquationInfo)::Nothing - x_dict = eqInfo.x_dict - der_x_dict = eqInfo.der_x_dict - startIndex = 1 - for (i, xi_info) in enumerate(eqInfo.x_info) - x_dict[xi_info.x_name] = i - der_x_dict[xi_info.der_x_name] = i - xi_info.startIndex = startIndex - startIndex += xi_info.length - end - eqInfo.nx = startIndex - 1 - return nothing -end - - -""" - x_init = initialStateVector(eqInfo::EquationInfo, FloatType) - -Return initial state vector `Vector{FloatType}` with stripped units. -""" -function initialStateVector(eqInfo::EquationInfo, FloatType::Type)::Vector{FloatType} - x_start = fill(FloatType(0), eqInfo.nx) - startIndex = 1 - for xe_info in eqInfo.x_info - if xe_info.scalar - if !(xe_info.startOrInit isa Nothing) - x_start[startIndex] = FloatType(ustrip(xe_info.startOrInit)) - startIndex += 1 - end - else - for i = 1:xe_info.length - x_start[startIndex] = FloatType(ustrip(xe_info.startOrInit[i])) - startIndex += 1 - end - end - end - return x_start -end - - -""" - x_start = updateEquationInfo!(eqInfo::EquationInfo, FloatType) - -Set eqInfo.x_dict, eqInfo.der_x_dict, eqInfo.nx and eqInfo.x_info[:].startIndex -""" -function updateEquationInfo!(eqInfo::EquationInfo, FloatType::Type)::Vector{FloatType} - nxFixedLength = eqInfo.nxFixedLength - if nxFixedLength == 0 - startIndex = 1 - else - xi_info = eqInfo.x_info[nxFixedLength] - startIndex = xi_info.startIndex + xi_info.length - end - - x_info = eqInfo.x_info - for i = nxFixedLength+1:length(x_info) - xi_info = x_info[i] - xi_info.length = length(xi_info.startOrInit) - xi_info.startIndex = startIndex - startIndex += xi_info.length - end - eqInfo.nx = startIndex - 1 - - return initialStateVector(eqInfo, FloatType) -end - - -get_x_startIndexAndLength(eqInfo::EquationInfo, name::String) = begin - xe_info = eqInfo.x_info[ eqInfo.x_dict[name] ] - (xe_info.startIndex, xe_info.length) -end - - -""" - nvec = get_equationSizes(equationInfo) - -Return the sizes of the linear equations as `nvec::Vector{Int}`. -(e.g. `nvec=[10,3]` means that there are two linear equation systems with -size 10 and size 3). -""" -get_equationSizes(eqInfo::EquationInfo) = Int[leq[3] for leq in eqInfo.linearEquations] - - -function Base.show(io::IO, eqInfo::EquationInfo; indent=4) - indentation = repeat(" ", indent) - indentation2 = repeat(" ", 2*indent) - indentation3 = repeat(" ", 3*indent) - println(io, "ModiaBase.EquationInfo(") - println(io, indentation2, "status = ModiaBase.", eqInfo.status, ",") - println(io, indentation2, "ode = ", eqInfo.ode, ",") - if eqInfo.nz > 0 - println(io, indentation2, "nz = ", eqInfo.nz, ",") - end - print(io, indentation2, "x_info = ModiaBase.StateElementInfo[") - for (i, xe_info) in enumerate(eqInfo.x_info) - if i == 1 - print(io, "\n", indentation3) - else - print(io, ",\n", indentation3) - end - show(io, xe_info) - end - print(io,"]") - - if length(eqInfo.residualCategories) > 0 - print(io, ",\n", indentation2, "residualCategories = [") - for (i, rcat) in enumerate(eqInfo.residualCategories) - if i > 1 - print(io, ",") - end - show(io, rcat) - end - print(io,"]") - end - - leqs = eqInfo.linearEquations - if length(leqs) > 0 - println(io, ",\n", indentation2, "linearEquations = [") - for (i, leq) in enumerate(leqs) - print(io, indentation3, "(", leq[1], ",\n", - indentation3, " ", leq[2], ",\n", - indentation3, " ", leq[3], ", ", leq[4], ")") - if i < length(leqs) - print(io, ",\n") - end - end - print(io, "]") - end - - - if length(eqInfo.vSolvedWithFixedTrue) > 0 - print(io, ",\n", indentation2, "vSolvedWithFixedTrue = ") - show(io, eqInfo.vSolvedWithFixedTrue) - end - - if eqInfo.nx > 0 - print(io, ",\n", indentation2, "nx = ", eqInfo.nx) - end - - if length(eqInfo.x_infoByIndex) > 0 - print(io, ",\n", indentation2, "x_infoByIndex = ") - show(io, eqInfo.x_infoByIndex) - end - - if !isnothing(eqInfo.defaultParameterAndStartValues) - print(io, ",\n", indentation2, "defaultParameterAndStartValues = ") - show(io, eqInfo.defaultParameterAndStartValues, indent=12, finalLineBreak=false) - end - - if !isnothing(eqInfo.ResultType) - print(io, ",\n", indentation2, "ResultType = ") - show(io, eqInfo.ResultType) - - if eqInfo.ResultTypeHasFloatType - print(io, ",\n", indentation2, "ResultTypeHasFloatType = ") - show(io, eqInfo.ResultTypeHasFloatType) - end - end - - - println(io, "\n", indentation, ")") - return nothing -end - - - -""" - names = get_stateNames(equationInfo::EquationInfo) - -Return the names of the states defined in `equationInfo` as a Vector of strings. -""" -get_stateNames(eqInfo::EquationInfo) = String[xi_info.x_name for xi_info in eqInfo.x_info] - - -""" - names = get_xNames(eqInfo::EquationInfo) - -Return the names of all elements of the x-vector as a vector of strings. -""" -function get_xNames(eqInfo::EquationInfo)::Vector{String} - xNames = Vector{String}(undef, eqInfo.nx) - for xe_info in eqInfo.x_info - if xe_info.length == 1 - xNames[xe_info.startIndex] = xe_info.x_name - else - for i = 1:xe_info.length - xNames[xe_info.startIndex+i-1] = xe_info.x_name*"["*string(i)*"]" - end - end - end - return xNames -end - - - -#= -""" - update_equationInfo!(eqInfo::EquationInfo) - -Update equationInfo as needed for [`EquationInfoStatus`](@ref)` = SOLVER_MODEL`. -""" -function update_equationInfo!(eqInfo::EquationInfo, modelValues::ModiaBase.AbstractModelValues)::Nothing - x_info = eqInfo.x_info - - if eqInfo.status == MANUAL - # Determine length and unit - count = 0 - for xi_info in x_info - if xi_info.x_name != "" # XD or XALG - (component,name) = ModiaBase.get_modelValuesAndName(modelValues, xi_info.x_name) - elseif xi_info.der_x_name != "" # XLAMBDA or XMUE - (component,name) = ModiaBase.get_modelValuesAndName(modelValues, xi_info.der_x_name) - else - # Pure algebraic with dummy equation der(x) = -x, x(0) = 0 - xi_info.length = 1 - count = count+1 - continue - end - - xi_type = fieldtype(typeof(component), name) - xi_numberType = ModiaBase.numberType(xi_type) - xi_info.unit = replace(string(unit(xi_numberType)), " " => "*") - xi_info.length = Base.length(xi_type) # Base.length(::Type{<:Number})=1 is defined in ModiaBase, since not defined in Base - end - if count == 1 && length(x_info) > 1 || count > 1 - error("Error in update_equationInfo!: x_info is wrong.\n", - "x_info = $x_info") - end - end - - # Set startIndex in x_info and compute nx - startIndex = 1 - for xi_info in x_info - xi_info.startIndex = startIndex - startIndex += xi_info.length - end - nx = startIndex - 1 - eqInfo.nx = nx - - # Set x_infoByIndex::Vector{Int} - eqInfo.x_infoByIndex = zeros(Int,nx) - i1 = 0 - i2 = 0 - for (i,xi_info) in enumerate(x_info) - i1 = xi_info.startIndex - i2 = i1 + xi_info.length - 1 - for j = i1:i2 - eqInfo.x_infoByIndex[j] = i - end - end - - eqInfo.status = SOLVER_MODEL - - return nothing -end -=# - diff --git a/src/ModiaBase.jl b/src/ModiaBase.jl index d5d704c..c8b4dfa 100644 --- a/src/ModiaBase.jl +++ b/src/ModiaBase.jl @@ -9,18 +9,13 @@ Main module of ModiaBase. module ModiaBase const path = dirname(dirname(@__FILE__)) # Absolute path of package directory -const Version = "0.9.2" -const Date = "2022-02-23" +const Version = "0.10.0" +const Date = "2022-03-01" #println("\nImporting ModiaBase Version $Version ($Date)") -using Unitful -using StaticArrays - - -# append! as needed in EquationAndStateInfo.jl and in ModiaLang/src/CodeGeneration.jl -appendVariable!(v1::Vector{FloatType}, s::FloatType) where {FloatType} = push!(v1,s) -appendVariable!(v1::Vector{FloatType}, v2) where {FloatType} = append!(v1,v2) +using Unitful +using StaticArrays include("LinearIntegerEquations.jl") @@ -41,7 +36,4 @@ using .Simplify include("Symbolic.jl") using .Symbolic -include("EquationAndStateInfo.jl") -include("StateSelection.jl") - end diff --git a/src/StateSelection.jl b/src/StateSelection.jl deleted file mode 100644 index 9ad7c8d..0000000 --- a/src/StateSelection.jl +++ /dev/null @@ -1,1579 +0,0 @@ -# License for this file: MIT (expat) -# Copyright 2017-2021, DLR Institute of System Dynamics and Control -# Author: Martin Otter, DLR-SR - - -#= -This file contains functions to transform the ODAE generated by the Pantelides -algorithm to an ODE, if this is reasonably possible -(e.g. without solving nonlinear equation systems inside the model or introducing stabilizing terms. -The approach is based on the dummy derivative method from -[Mattsson and Söderlind (1992](https://ieeexplore.ieee.org/document/274429), -the transformation algorithm from -[Otter and Elmqvist (2017)](https://modelica.org/events/modelica2017/proceedings/html/submissions/ecp17132565_OtterElmqvist.pdf), -and on some new ideas that are not yet published. - -# Functions exported by this module: - -- [`getSortedAndSolvedAST`](@ref): Given the BLT - (after Pantelides algorithm), - this function returns the Abstract Syntax Tree (AST) of - the sorted and solved equations and selects the ODE states. -=# - -export getSortedAndSolvedAST -export StateSelectionFunctions -export showCodeWithoutComments - - -#using OrderedCollections -using ModiaBase -using ModiaBase.BLTandPantelides -using StaticArrays - - - -const INFORMATION = 0 -const WARNING = 1 -const ERROR = 2 - - - -""" - stateSelectionFunctions = StateSelection.StateSelectionFunctions(; - var_name = v -> nothing, - var_julia_name = v -> nothing, - var_unit = v -> "", - var_startInitFixed = v -> (nothing,false), # (startOrInit, fixed) - var_is_state = v_original -> false, - equation = e -> "" - isSolvableEquation = (e_original,v_original) -> false, - isLinearEquation = (e_original,v_original) -> (false, false), - getSolvedEquationAST = (e, v) -> nothing, - getResidualEquationAST = (e, residualName) -> nothing, - showMessage = (message; severity=0,from="???",details="",variables=Int[],equations=Int[]) -> println("...") - ) - -Generate an instance of an immutable struct that holds the callback functions -as needed by [`getSortedAndSolvedAST`](@ref). - -The functions need to have the following arguments: - -- `var_name(v::Int)::String`:\\ - Return full (Modia) path name of `v` as String (e.g. `"a.b.c"` or `"a.b.der(c)"`). - -- `var_julia_name(v::Int)::String`:\\ - Return full Julia path name of variable `v` as `Symbol` - (e.g. Symbol("m.a.b.c") or Symbol("a_b_der_c")). - -- `var_unit(v::Int)::String`:\\ - Return unit of variable `v` as `String` (for example "N*m") or `""`, if no unit is defined. - -- `var_startInitFixed(v::Int)::Tuple(Any,Bool)`:\\ - Return `(startOrInit, fixed)`, where `startOrInit` is the `start` or `init` value or `nothing` - (if neither `start` nor `init` defined) and `fixed` is `true`, if an `init` value is defined. - -- `var_is_state(v_original::Int)::Bool`:\\ - Return true, if variable `v_original` is defined to be a state. - -- `equation(e::Int)`:\\ - Return equation as string or "", if equation is not provided. - -- `isSolvableEquation(e_original::Int, v_original::Int)::Bool`:\\ - Return `true` if original, undifferentiated equation `e_original` - can be uniquely solved for variable `v_original`. - -- `(isLinear, hasConstantCoefficient) = isLinearEquation(e_original::Int, v_original::Int)`:\\ - Return the information whether the original, undifferentiated equation `e` - is **linear** with respect to variable `v` (`isLinear=true`) and if this is the case whether - the coefficient of this linear term is **constant** after initialization - (`hasConstantCoefficient=true`). If `v` is an array, then the array equation `e` must - be linear with respect to all elements of array `v` in order to return `isLinear=true`. - -- `AST = getSolvedEquationAST(e:Int, v:Int)`:\\ - Return Abstract Syntax Tree `AST::Expr` of equation `e` solved for variable `v` - as instance of `Expr`. - -- `AST = getResidualEquationAST(e::Int, residualName::Symbol)`:\\ - Return Abstract Syntax Tree `AST::Expr` of eqation `e` as residual equation with ustrip(..) \\ - `residualName = ustrip( < residual form of equation e> )` - -- `showMessage(message; severity=0,from="???",details="",variables=Int[],equations=Int[])`:\\ - Print information, warning or error message. -""" -struct StateSelectionFunctions - var_name::Function - var_julia_name::Function - var_unit::Function - var_startInitFixed::Function - var_nominal::Function - var_unbounded::Function - var_is_state::Function - equation::Function - isSolvableEquation::Function - isLinearEquation::Function - getSolvedEquationAST::Function - getResidualEquationAST::Function - showMessage::Function - - StateSelectionFunctions(; - var_name = v -> nothing, - var_julia_name = v -> nothing, - var_unit = v -> "", - var_startInitFixed = v -> (nothing,false), # (startOrInit, fixed) - var_nominal = v_original -> NaN, - var_unbounded = v_original -> false, - var_is_state = v_original -> false, - equation = e -> "", - isSolvableEquation = (e_original,v_original) -> false, - isLinearEquation = (e_original,v_original) -> (false, false), - getSolvedEquationAST = (e, v) -> nothing, - getResidualEquationAST = (e, residualName) -> nothing, - showMessage = (message; severity=0,from="???",details="",variables=Int[],equations=Int[]) -> - println("\n!!! Message from $from with severity=$severity:\n$message,\ndetails\nvariables=$variables,\nequations=$equations") - ) = new(var_name, - var_julia_name, - var_unit, - var_startInitFixed, - var_nominal, - var_unbounded, - var_is_state, - equation, - isSolvableEquation, - isLinearEquation, - getSolvedEquationAST, - getResidualEquationAST, - showMessage) -end - - -function printVector(vec) - for v in vec - println(" ", v) - end -end - - - -""" - showMessage(message, severity, from, details, variables, equations) - -Print message. - -# Arguments: - -- `severity::Int`: Severity of message: 0=Information, 1=warning, 2=Error -- `from::String`: Name of function that throws error message -- `message::String`: One line description of issue -- `details::String`: Detailed description of issue -- `variables::Vector{String}`: List of involved variables -- `equations::Vector{String}`: List of involved equations -""" -function showMessage(message::String, severity::Int, from::String, details::String, - variables::AbstractVector, equations::AbstractVector)::Nothing - if severity == 0 - println("\nInformation message from $(from):") - elseif severity == 1 - println("\nWarning message from $(from):") - else - println("\nError message from $(from):") - end - println(message) - if length(details) > 0 - println(details) - end - if length(variables) > 0 - println("Involved variables:") - printVector(variables) - end - if length(equations) > 0 - println("Involved equations:"), - printVector(equations) - end - println("") - return nothing -end - - - -showCodeWithoutComments(code) = println("code = ", replace(sprint(show,code), r"( *#= .*=#\n)|(#= .*=#)" => "") ) - - - -""" - (eConstraints,vConstraints) = getConstraintSets(BLT_i, assignRev, Arev, Brev, - var_is_state, showMessage) - -Determines the set of equations/constraints and their unknowns that are associated with -the highest derivative level BLT block i (BLT_i), according to the dummy derivative -method of [Mattsson and Söderlind (1992](https://ieeexplore.ieee.org/document/274429), -(note, eConstraints[j] and vConstraints[j] are of type Vector{Int}): - -- eConstraints[end] = BLT_i, so the equations of BLT block i (must be highest derivative equations) -- eConstraints[j]: Constraint equations that form the integral of eConstraints[j+1]. -- eConstraints[1]: Lowest order constraint equations. All equations in eConstraints[1] - are original, non-differentiated equations. - -vConstraints[j] are the unknowns of eConstraints[j]. - -# Arguments: - -- `BLT_i::Vector{Int}`: Vector of equation indices of BLT block i. - -- `assignRev::Vector{Int}`: Reverted assign vector: vi = assignRev[ei] is variable vi assigned to equation ei - -- `Arev::Vector{Int}`: Reverted variable association: `Arev[i] = der(V[k]) == V[i] ? k : 0`. - -- `Brev::Vector{Int}`: Reverted equation association: `Erev[i] = der(E[k]) == E[i] ? k : 0`. - -- `var_is_state::Function`: var_is_state(`vi`) returns true, if `vi` is defined to be a definite state. - -- showMessage: showMessage2 below to print error/warning message - -""" -function getConstraintSets(BLT_i::Vector{Int}, - assignRev::Vector{Int}, - Arev::Vector{Int}, - Brev::Vector{Int}, - var_is_state::Function, - showMessage::Function) - # Determine unknowns of BLT block i - BLT_i_unknowns = fill(0, length(BLT_i)) - for j in eachindex(BLT_i) - BLT_i_unknowns[j] = assignRev[BLT_i[j]] - end - - # Initialize constraint sets - eConstraints = [BLT_i] - vConstraints = [BLT_i_unknowns] - - # Check that number of equations and number of unknowns agree - if length(BLT_i) != length(BLT_i_unknowns) - showMessage("Number of equations and number of unknowns do not agree for BLT block"; - severity = ERROR, - variables = BLT_i_unknowns, - equations = BLT_i) - return (eConstraints, vConstraints, false) - end - - while true - # Determine constraints at one differentiation order less - ceq = fill(0, 0) - for eq in eConstraints[1] - if Brev[eq] > 0 - push!(ceq, Brev[eq]) - end - end - if length(ceq) == 0; break; end - - pushfirst!(eConstraints, ceq) # move ceq to the beginning of the constraints vector - - # Determine unknowns of ceq - veq = fill(0,0) - veq_states = fill(0,0) - for vc in vConstraints[1] - if Arev[vc] > 0 - if var_is_state( Arev[vc] ) - push!(veq_states, Arev[vc]) - else - push!(veq, Arev[vc]) - end - end - end - if length(veq) < length(ceq) - nRemoveStates = length(veq_states) - (length(ceq) - length(veq)) - showMessage("$nRemoveStates of the followig defined states cannot be states!"; - severity = ERROR, - variables = veq_states, - equations = ceq) - return (eConstraints, vConstraints, false) - end - pushfirst!(vConstraints, veq) # move veq to the beginning of the constraints vector - end - - @assert(length(eConstraints) == length(vConstraints)) - return (eConstraints, vConstraints, true) -end - - - -""" - names = getNames(var_name::Function, v::Vector{Int}) - names = getNames(eq::EquationGraph , v::Vector{Int}) - -Get names of variables `v` as a string. -""" -function getNames(var_name::Function, v)::String - if length(v) > 0 - str = var_name(v[1]) - for i in 2:length(v) - str = str * ", " * var_name(v[i]) - end - return str - else - return "" - end -end -getNames(eq, v) = getNames(eq.fc.var_name, v) - - - -function printVariables(eq, v; indent=4) - for vi in v - println(lpad(vi,indent), ": ", eq.fc.var_name(vi)) - end - return nothing -end - - -""" - printEquations(equation::Function, e::Vector{Int}; indent=4) - printEquations(eq::EquationGraph, e::Vector{Int}; indent=4) - printEquations(eq::EquationGraph, expr; indent=4) - -Print equations `e` as a string or as vector of equations. -""" -function printEquations(equation::Function, e::Vector{Int}; indent=4)::Nothing - for ei in e - println(lpad(ei,indent), ": ", equation(ei)) - end - return nothing -end -printEquations(eq, e; indent=4) = printEquations(eq.fc.equation, e; indent=indent) - -function printEquations(expr::Expr; indent=4) - indentation = repeat(" ",indent) - str = string(expr) - str = replace(str, r".*#= .*\n" => "") - println(indentation, str) -end - -#= - -""" - Arev = revertAssociation(A::Vector{Int}) - -Revert the association Vector `A[i] = j`, such that -`Arev[j] = i` (`A[i]=0` is allowed and is ignored). -""" -function revertAssociation(A::Vector{Int})::Vector{Int} - nArev = max(maximum(A), length(A)) - Arev = fill(0, nArev) - for i in eachindex(A) - if A[i] != 0 - Arev[ A[i] ] = i - end - end - return Arev -end -=# - - -""" - tearedEquations = TearedEquations(..) - -Return an immutable struct that contains the information about a teared system of equations. -""" -struct TearedEquations - isLinear::Bool # = true if equation system is linear - hasConstantCoefficients::Bool # = true if linear equation system has a constant coefficient matrix A - vSolved::Vector{Int} # Variables that are explicitly solved by tearing - vTear::Vector{Int} # Iteration variables of tearing - eSolved::Vector{Int} # Equations that are explicitly solved by tearing - eResiduals::Vector{Int} # Residual equations of the teared equations - AST::Expr # AST of equation system -end - - -""" - eq = EquationGraph(G, BLT, assign, A, B, log, showMessage::Function) - -Construct the representation of the equation graph as needed by state selection. - - -# Arguments - -- `G`: G[i] is the Integer vector of variable indices of equation i. - A variable i can be a scalar or an array. - G is of type Vector{ Vector{Int} } or Vector{ Any }. - -- `BLT::Vector{Vector{Int}}`: BLT[i] is the vector of equations belonging to BLT-block i. - In principal, **BLT[i] must include only highest derivative equations**. - In order to ease the generation of BLT[i], to not be forced to map indices forth and back - when calling the function that computes the BLT, the following assumption is - made: If BLT[j] contains only **one equation with B[j][1] > 0** then this BLT block is - a lower-derivative (dummy) block and is ignored. - -- `assign::Vector{Int}`: ei = assign[vi] is equation ei assigned to variable vi - -- `A`: A-Vector of Pantelides algorithm: - `A[i] = if der(v[i]) == v[k] then k else 0` - where `v[i]` is variable `i`. - -- `B`: B-Vector of Pantelides algorithm: - `B[i] = if der(e[i]) == e[k] then k else 0` - where `e[i]` is equation `i`. - -- `log`: = true, if logging is on - -- `showMessage`: Function used to print warning/error messages. -""" -mutable struct EquationGraph - success::Bool # = true, if operation was successful, otherwise false - log::Bool # = true, if logging is on - logDetails::Bool # = true, if detailled logging is on - showMessage::Function - fc::StateSelectionFunctions # Callback functions - G # Structure of equations - A::Vector{Int} - B::Vector{Int} - Gunknowns::Vector{Vector{Int}} - vOriginal::Vector{Bool} # vOriginal[v] = true, if variable v appears in an original (non-differentiated) equation - vActive::Vector{Bool} # vActive[v] = false, if v is a state - vSolvedWithFixedTrue::Vector{Int} # Explicitly solved variables with fixed=true - - # Reverted info vectors - Arev::Vector{Int} # Reverted vector A - Brev::Vector{Int} # Reverted vector B - assignRev::Vector{Int} # Reverted vector assign; vi = assignRev[ei] if equation ei is assigned to variable vi - - # Equation/Constraint sets - eConstraintsVec::Vector{Vector{Vector{Int}}} # eConstraintsVec[i] is the constraint set of highest derivative BLT block i - vConstraintsVec::Vector{Vector{Vector{Int}}} # vConstraintsVec[i] are the assigned variables of eConstraintsVec[i] - - # Tearing - td::TearingSetup # Data structure for equation tearing to traverse a set of equations - # that are represented as a DAG (Directed Acyclic Graph). - vSolved::Vector{Int} # Variables that are explicitly solved by tearing - vTear::Vector{Int} # Iteration variables of tearing - eSolved::Vector{Int} # Equations that are explicitly solved by tearing - eResiduals::Vector{Int} # Residual equations of the teared equations - vTearWithoutStart::Vector{Int} # Summary of all vTear variables that have no start value (to print warning message at the end) - - der_vSolvedLower::Vector{Int} # Initialize tearing with variables that are already explicitly solved - # (this is known from the constraints one differentiation level less) - der_vTearLower::Vector{Int} # Highest priority variables to solve for - # (= differentiated tearing variables from one differentiation level less) - # (this is known from the constraints one differentiation level less) - der_eSolvedLower::Vector{Int} # Initialize tearing with equations that are already explicitly solved - # (this is known from the constraints one differentiation level less) - vc::Vector{Int} # Variables of the current component that are unknowns of the tearing - ec::Vector{Int} # Equations of the current component that shall be teared - vcCandidates::Vector{Vector{Int}} # Try to solve equations ec in the following order for vc: - # vcCandidates[1], ..., [6] - vLinActive::Vector{Bool} # Utility vector (vLinActive[i]=true, if variable is inspected) - equationInfo::ModiaBase.EquationInfo # Info about the structure of the equations - # (equationInfo.x_names, equationInfo.linearEquations) - - # Sorted equations - fullAssignRev::Vector{Int} # fullAssignRev[e] = v: If e is a scalar equation system, it is solved for v. - # If v = -1, the scalar equation system is not solved directly, but - # as linear equation system (with tearing). - tearedEquations::Vector{TearedEquations} # Info about teared equation system - tearedEquations_indices::Vector{Int} # teq = tearedEquations_indices[e] is tearedEquations[teq] - eAST::Vector{Expr} # eAST[e] is the AST of equation e - AST_aux::Vector{Expr} # Abstract Syntax Tree of some code (auxiliary vector) - AST::Vector{Expr} # Abstract Syntax Tree of the sorted and solved equations as vector of Expr - # Statements are evaluted in the order: AST[1], AST[2], .... - - function EquationGraph(G, BLT, assign, A, B, stateSelectionFunctions, log, logDetails, showMessage::Function) - - # Check input arguments - @assert(length(B) == length(G)) - - # Revert association vectors - Arev = revertAssociation(A) - Brev = revertAssociation(B) - assignRev = revertAssociation(assign) - - # Initialize vector of variables of the non-differentiated, original equations - vOriginal = fill(false, length(A)) - - # Initialize the active variables - vActive = fill(true, length(A)) - for v = 1:length(vActive) - if Arev[v] > 0 - vActive[ Arev[v] ] = false # Arev[v] is a potential state - end - end - - # Define empty constraint sets - eConstraintsVec = Vector{Vector{Int}}[] - vConstraintsVec = Vector{Vector{Int}}[] - - # Gunknowns has the same dimension as G. - # Gunknowns[i] is not empty, if it is an original, non-differentiated - # equation. In this case it is a subset of G[i] containing only the - # unknown variables of the corresponding constraint set as needed by tearing. - # Note, tearing is only applied on non-differentiated equations - # Tearing information of differentiated equations is - # deduced from the tearing of the non-differentiated equations. - Gunknowns = fill(Int[], length(G)) - - # Construct constraint sets of the dummy derivative method from the BLT - success = true - for (i,BLT_i) in enumerate(BLT) - # Check that BLT block contains only higher derivative equations - # (with one exception due to an agreed convention with Modia: - # If BLT[j] has only one equation and B( BLT[j][1] ) > 0, so it is - # a lower derivative equation, then this BLT block is ignored). - highest = true - for eqj in BLT_i - if B[eqj] != 0 # One equation of block BLT_i appears also differentiated - highest = false - if length(BLT_i) > 1 - success = false - showMessage("Equation $eqj in BLT[$i]=$BLT_i is not a highest derivative equation."; - details = "This means that there is an error in argument BLT.", - severity = ERROR, - equations = [eqj]) - end - break - end - end - if !success - break - end - - if highest - # Get all equation sets eConstraints and their corresponding unknowns vConstraints - # from lowest to highest differentiation order (eConstraints[end] is c) - (eConstraints, vConstraints, success2) = getConstraintSets(BLT_i, assignRev, Arev, Brev, - stateSelectionFunctions.var_is_state, showMessage) - if !success2 - success=false - break - end - - # Initialize vActive - #for vc in vConstraints - # for v in vc - # if Arev[v] > 0 - # vActive[ Arev[v] ] = false # Arev[v] is a potential state - # end - # end - #end - - push!(eConstraintsVec, eConstraints) - push!(vConstraintsVec, vConstraints) - - # Determine unknown variables of G as used by tearing - for i in eachindex(eConstraints) - for eq in eConstraints[i] - # Gunknowns[eq] contains the variables of G[eq] that are unknowns of vConstraints[i] - Gunknowns[eq] = intersect(G[eq], vConstraints[i]) - - if Brev[eq] == 0 - # Original, non-differentiated equation - for v in Gunknowns[eq] - # v is a variable from the original, non-differentiated equations - vOriginal[v] = true - end - end - end - end - end - end - - new(success, log, logDetails, showMessage, stateSelectionFunctions, G, A, B, Gunknowns, vOriginal, vActive, - Int[], Arev, Brev, assignRev, - eConstraintsVec, vConstraintsVec, - TearingSetup(Gunknowns,length(Arev)), - Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[], Int[], - Vector{Int}[Int[], Int[], Int[], Int[], Int[], Int[]], fill(false, length(Arev)), - ModiaBase.EquationInfo(status=ModiaBase.CODE_GENERATION), - fill(0, length(B)), TearedEquations[], fill(0, length(B)), - Vector{Expr}(undef, length(B)), Expr[], Expr[]) - end -end - - -function undifferentiated(eq::EquationGraph, v::Int)::Int - while eq.Arev[v] > 0 - v = eq.Arev[v] - end - return v -end - - -""" - isLinearOriginalOrDifferentiatedEquation(eq, e, v) -""" -function isLinearOriginalOrDifferentiatedEquation(eq::EquationGraph, e::Int, v::Int) - # Determine original, undifferentiated equation - while eq.Brev[e] > 0 - e = eq.Brev[e] - v = eq.Arev[v] - end - - # Check whether original, undifferentiated equation is linear - return eq.fc.isLinearEquation(e, v) -end - - -""" - (isLinear, hasConstantCoefficients) = isLinearEquationSystem!(eq, es, vs) - -Returns `isLinear = true`, if equations `es::Vector{Int}` are linear with respect to -variables `vs::Vector{Int}`. If **all** coefficients of the linear factors are additionally constant, -`hasConstantCoefficients=true` is returned, otherwise false. -`eq::EquationGraph` is the data structure holding all information about the graph. -""" -function isLinearEquationSystem!(eq::EquationGraph, es::Vector{Int}, vs::Vector{Int})::Tuple{Bool,Bool} - # It is assumed, that on input, eq.vLinActive[i] = false - # On return, eq.vLinActive[i] = false again - - # Mark that vs are active variables - eq.vLinActive[vs] .= true - - # Check whether es is linear with respect to vs - isLinear::Bool = true - hasConstantCoefficients::Bool = true - for e in es - for v in eq.G[e] - if eq.vLinActive[v] # Only check variables vs in G[e] - (e_isLinear, e_hasConstantCoefficient) = isLinearOriginalOrDifferentiatedEquation(eq, e, v) - if !e_isLinear - isLinear = false - hasConstantCoefficients = false - break - else - hasConstantCoefficients = hasConstantCoefficients && e_hasConstantCoefficient - end - end - end - if !isLinear - break - end - end - - # Unmark that vs are active variables - eq.vLinActive[vs] .= false - - return (isLinear, hasConstantCoefficients) -end - - - -""" - tearEquationsWithCandidates!(eq::EquationGraph) - -Tear equations (eq.der_eSolvedLower, eq.ec) with respect to variables -(eq.der_vSolvedLower, eq.vc) and initialize tearing with the -Directed-Acylic-Graph (eq.der_eSolvedLower, eq.der_vSolvedLower). - -It is tried to solve the equations eq.ec with respect to variables `eq.vc` -in the following order: - -1. eq.der_vTearLower\\ - (These are differentiated variables from the previous differentiation level that - have been used as tearing variables at the previous differentation level. - In general, it is adviceable to use non-differentiated - variables as states and by trying to eliminate der_vTearLower first, it is tried to solve - for all differentiated variables first; note these variables have no start values). - -2. Variable has start=nothing and the variable does not appear differentiated - in the set of original equations. - (Variables without a start value should be explicitly computed during initialization and - therefore these variables should be solved for; - however, if a variable appears also differentiated in the original equations, - then it should be tried to use it as a state and therefore, it should be first tried - to not solve for this variable). - -3. Variable has start=nothing and the variable appears differentiated - in the set of original equations. - (Variables without a start value should be explicitly computed during initialization and - therefore these variables should be solved for). - -4. Variable has fixed=false\\ - (The start value is a guess value, so it is safe to solve for the variable and - hereyby ignore the start value) - -5. Variable has fixed=nothing\\ - (The intention is to respect the start value if possible, that is, it should be avoided - to solve for such a variable) - -6. Variable has fixed=true\\ - (The start value must be respected, that is, if any possible no equation should be solved - for this variable; if it is solved for, an error is triggered during initialization if the - start value and the computed value do not agree) - - -# Example - -Assume the constraint `phi1 = phi2` and differentiated constraints -`der(phi1) = der(phi2), w1 = der(phi1), w2 = der(phi2)` and start values are defined -for `phi1` and `w1`. - -- Since no start value is defined for `phi2` (but a start value is defined for `phi1`), - equation `phi1 = phi2` is solved for `phi2`, so `phi2 := phi1` and `phi1` is the tearing variable. - This means that `der_vSolvedLower = der(phi1); der_vTearLower = der(phi2)`. - -- On the next differentiation level, the differentiated equation `der(phi1) = der(phi2)` - is solved in the same way as on the previous level due to `der_eSolvedLower, der_vSolvedLower`, that is - `der(phi2) := der(phi1)`. - -- Since `der_vTearLower = der(phi1)`, it is first tried to solve for `der(phi1)`. There is only - one equation left with `der(phi1)` and therefore this equation is solved for it: - `der(phi1) := w1`. - -- Since a start value is defined for `w1` and no start value for `w2`, it is tried to solve - for `w2`, that is: `w2 := der(phi2)` and the result is that `w1` is used as tearing variable. - -Resulting equation set with states = `phi1, w1`: - -```julia -phi2 := phi1 -der(phi1) := w1 -der(phi2) := der(phi1) -w2 := der(phi2) -``` -""" -function tearEquationsWithCandidates!(eq::EquationGraph)::Nothing - # Empty all candidate vectors - for c in eq.vcCandidates - empty!(c) - end - - # Determine priorities to solve for variables eq.vc - # Assume equations: - # (1) s = ... - # (1') w = der(s) - # (1'') a = der(w) - # Then it is most natural to use "w" and "a" as tearing variables if "s" is a tearing variable - # and solve for "der(phi)" and "der(w)". - # This means, that eq.der_vTearLower should be tried first to be solved - eq.vcCandidates[1] = eq.der_vTearLower - - for v in setdiff(eq.vc, eq.der_vTearLower) - # Variable has no start value because: - # (a) is not from the original, non-differentiated equations - # (b) is a differentiated variable - # (c) no start or init value is defined - (startOrInit, fixed) = eq.fc.var_startInitFixed(v) - hasNoStart = !eq.vOriginal[v] || eq.Arev[v] > 0 || startOrInit isa Nothing - if hasNoStart - # Has no start value - if eq.vOriginal[v] && (eq.A[v]>0 && eq.vOriginal[eq.A[v]]) - # Is from the original, non-differentiated equations and its derivative - # is also in the original equations (so try to use it as state). - push!(eq.vcCandidates[3], v) - else - # All other variables without a start value are tried - push!(eq.vcCandidates[2], v) - end - - else - if startOrInit isa Nothing - push!(eq.vcCandidates[5], v) - elseif fixed - push!(eq.vcCandidates[6], v) - else - push!(eq.vcCandidates[4], v) - end - end - end - - if eq.log && eq.logDetails - println(" eSolvedFixed = ", eq.der_eSolvedLower) - println(" vSolvedFixed = ", eq.der_vSolvedLower) - println(" vcCandidates = ", eq.vcCandidates) - end - - # Tear equations (ec are non-differentiated, original equations) - (eq.eSolved, eq.vSolved, eq.eResiduals, eq.vTear) = - tearEquations!(eq.td, eq.fc.isSolvableEquation, eq.ec, eq.vcCandidates; - eSolvedFixed = eq.der_eSolvedLower, - vSolvedFixed = eq.der_vSolvedLower, - log = eq.log && eq.logDetails) - return nothing -end - - - -""" - v = derivativeOfVlower(vLower::Vector{Int}, A::Vector{Int}) - -Given a set of variables vLower and the A-Vector of Pantelides, -return the derivative of vLower (so the differentiated variables of vLower). -""" -function derivativeOfVlower(vLower::Vector{Int}, A::Vector{Int})::Vector{Int} - v = fill(0, length(vLower)) - for i in eachindex(vLower) - derv = A[ vLower[i] ] - @assert(derv > 0) - v[i] = derv - end - return v -end - - - - -""" - addFixedVariable!(eq::EquationGraph, v::Int) - -If variable v is not differentiated AND has fixed=true store it in eq.equationInfo -(this is used to check during initialization, that the start value is respected) -""" -function addFixedVariable!(eq::EquationGraph, v::Int)::Nothing - (v_startOrInit,v_fixed) = eq.fc.var_startInitFixed(v) - if eq.Arev[v] == 0 && v_fixed - push!(eq.vSolvedWithFixedTrue, v) - push!(eq.equationInfo.vSolvedWithFixedTrue, eq.fc.var_name(v)) - end - return nothing -end - - - -""" - addSolvedEquations!(eq::EquationGraph, eSolved::Vector{Int}, vSolved::Vector{Int}) - -Solve equation eSolved[i] for variable vSolved[i] and push this information to -eq.SortedEquations. Furthermore, generate the AST. -""" -function addSolvedEquations!(eq::EquationGraph, eSolved::Vector{Int}, vSolved::Vector{Int})::Nothing - if eq.log - println("Julia code:") - end - - for i in eachindex(eSolved) - # Push solved equation to eAST - v = vSolved[i] - e = eSolved[i] - eq.eAST[e] = eq.fc.getSolvedEquationAST(e, v) # filter_lineno( ) - if eq.log - printEquations(eq.eAST[e]) - end - eq.fullAssignRev[e] = v - addFixedVariable!(eq,v) - end - - return nothing -end - - - -""" - addLinearEquations!(eq::EquationGraph, hasConstantCoefficients::Bool, unitless::Bool) - -A linear equation system is solved. Push this information to -eq.SortedEquations. Furthermore, generate the AST for a for-loop that builds -and solves a linear equation system A*x = b -from the solved equations and the residual equations, -and push this information to eq.SortedEquations.AST. -""" -function addLinearEquations!(eq::EquationGraph, hasConstantCoefficients::Bool, unitless::Bool)::Nothing - # Construct body of for-loop - empty!(eq.AST_aux) - while_body = eq.AST_aux - vTear_names = String[] - vTear_lengths = Int[] - - # Assign iteration variables - # unitless = false: - # v_i = leq.x[i] - # - # unitless = true: - # v_i = leq.x[i]*@u_str($v_unit) - # - # assign vector-valued variables vec_j at the end (storage is allocated inside LinearEquations) - # vec_j = leq.x_vec[j] - # - vAssigned_names = Any[] - v_vec_julia_names = Any[] - v_vec = Int[] - v_vec_startOrInit = Any[] - i1 = 0 - i2 = 0 - for (i,v) in enumerate(eq.vTear) - (v_startOrInit, v_fixed) = eq.fc.var_startInitFixed(undifferentiated(eq,v)) - v_name = eq.fc.var_name(v) - - if isFixedLengthStartOrInit(v_startOrInit, v_name) - # length(v_startOrInit) is fixed after compilation - v_julia_name = eq.fc.var_julia_name(v) - v_unit = unitless ? "" : eq.fc.var_unit(v) - v_length = isnothing(v_startOrInit) ? 1 : length(v_startOrInit) - i1 = i2 + 1 - i2 = i1 + v_length - 1 - if v_startOrInit isa AbstractVector - # v is a vector - if v_unit == "" - push!(while_body, :( $v_julia_name::ModiaBase.SVector{$v_length,_FloatType} = ModiaBase.SVector{$v_length,_FloatType}(_leq_mode.x[$i1:$i2])) ) - else - push!(while_body, :( $v_julia_name = ModiaBase.SVector{$v_length}(_leq_mode.x[$i1:$i2])::ModiaBase.SVector{$v_length,_FloatType}*@u_str($v_unit)) ) - end - elseif v_startOrInit isa Number || v_startOrInit isa Nothing - # v is a scalar or nothing (= assumed to be a scalar) - if v_unit == "" - push!(while_body, :( $v_julia_name = _leq_mode.x[$i1] ) ) - else - push!(while_body, :( $v_julia_name = _leq_mode.x[$i1]*@u_str($v_unit) ) ) - end - else - error("Should not occur: v_name = $v_name, v_startOrInit = $v_startOrInit") - end - push!(vAssigned_names, v_julia_name) - push!(vTear_lengths , v_length) - - # Store iteration variables in linearEquations to be used in error messages - push!(vTear_names, v_name) - else - # length(v_startOrInit) may change after compilation - push!(v_vec, v) - push!(v_vec_startOrInit, v_startOrInit) - end - end - - # Assignment code for vector-valued variables - nx_fixedLength = length(vAssigned_names) - for (i,v) in enumerate(v_vec) - v_startOrInit = v_vec_startOrInit[i] - v_julia_name = eq.fc.var_julia_name(v) - v_unit = unitless ? "" : eq.fc.var_unit(v) - if v_unit == "" - push!(while_body, :( $v_julia_name = _leq_mode.x_vec[$i] ) ) - else - push!(while_body, :( $v_julia_name = _leq_mode.x_vec[$i]*@u_str($v_unit) )) - end - push!(vAssigned_names , v_julia_name) - push!(vTear_lengths , length(v_startOrInit)) - push!(v_vec_julia_names, v_julia_name) - - # Store iteration variables in linearEquations to be used in error messages - push!(vTear_names, eq.fc.var_name(v)) - end - - # Add solved equations - for (i,e) in enumerate(eq.eSolved) - push!(while_body, eq.fc.getSolvedEquationAST(e, eq.vSolved[i])) - push!(vAssigned_names, eq.fc.var_julia_name(eq.vSolved[i])) - end - - # Add residual equations - # appendResidual!(leq.residuals, < equation in residual form >) - for (i,e) in enumerate(eq.eResiduals) - e_AST = eq.fc.getResidualEquationAST(e, :(_leq_mode.residuals[$i]) ) - if !isnothing(e_AST) - push!(while_body, e_AST) - end - end - - # Generate LinearEquations data structure - push!(eq.equationInfo.linearEquations, (vTear_names, v_vec_julia_names, vTear_lengths, nx_fixedLength, hasConstantCoefficients)) - - # Construct for-loop - leq_index = length(eq.equationInfo.linearEquations) - while_loop = quote - local $(vAssigned_names...) - _leq_mode = initLinearEquationsIteration!(_m, $leq_index) - ModiaBase.TimerOutputs.@timeit _m.timer "ModiaBase LinearEquationsIteration!" while ModiaBase.LinearEquationsIteration!(_leq_mode, _m.isInitial, _m.solve_leq, _m.storeResult, _m.time, _m.timer) - $(while_body...) - end - _leq_mode = nothing - end - - if eq.log - showCodeWithoutComments(while_loop) - end - - # Store information about linear equation system - teq = TearedEquations(true, hasConstantCoefficients, - copy(eq.vSolved), copy(eq.vTear), - copy(eq.eSolved), copy(eq.eResiduals), - while_loop) - push!(eq.tearedEquations, teq) - teq_index = length(eq.tearedEquations) - - # Reference teq_index for all equations - for eSolved in eq.eSolved - eq.tearedEquations_indices[eSolved] = teq_index - end - for eResiduals in eq.eResiduals - eq.tearedEquations_indices[eResiduals] = teq_index - end - - # Store fixed=true variables - for vTear in eq.vTear - addFixedVariable!(eq,vTear) - end - for vSolved in eq.vSolved - addFixedVariable!(eq,vSolved) - end - - # Store vTear without a start value - for vTear in eq.vTear - (vTear_startOrInit, vTear_fixed) = eq.fc.var_startInitFixed(vTear) - if vTear_startOrInit isa Nothing - push!(eq.vTearWithoutStart, vTear) - end - end - - return nothing -end - - - -""" - sortEquations!(eq:EquationGraph) - -Generate the sorted equations -""" -function sortEquations!(eq::EquationGraph, Goriginal)::Nothing - if eq.log - println("Sort equations (BLT on all equations under the assumption that the ODE states are known).") - # println("\nSorted equations (BLT on all equations under\nthe assumption that the ODE states are known):") - # - # code = quote - # $(eq.AST...) - # end - # showCodeWithoutComments(code) - end - - # Matching of all equations - assign = matching(Goriginal, length(eq.A), eq.vActive) - assignRev = revertAssociation(assign) - - # Sort equations and find strong components - blt = BLT(Goriginal,assign) - - # Generate AST of sorted equations - for blt_i in blt - @assert(length(blt_i) >= 1) - - solveLinearEquation = true - if length(blt_i) == 1 # One equation in one unknown - e = blt_i[1] - - # Check - v = eq.fullAssignRev[e] - if v != -1 - @assert(e == assign[v]) - - # Push AST - push!(eq.AST, eq.eAST[e]) - solveLinearEquation = false - end - end - - if solveLinearEquation # Equation system - # Check that all equations reference the same equation system - e1 = blt_i[1] - teq_index = eq.tearedEquations_indices[e1] - @assert(teq_index != 0) - for j=2:length(blt_i) - e = blt_i[j] - if eq.tearedEquations_indices[e] != teq_index - error("... should not occur: teq_index=$teq_index, e1=$e1, e_teq_index=", eq.tearedEquations_indices[e]) - end - end - - # Check that the assigned variables of the BLT block are the variables computed from the equation system - teq = eq.tearedEquations[teq_index] - vTearComputed = union(teq.vTear, teq.vSolved) - vBltComputed = Int[] - for e in blt_i - push!(vBltComputed, assignRev[e]) - end - vdiff = setdiff(vBltComputed, vTearComputed) - @assert(length(vdiff) == 0) - - # Push AST - push!(eq.AST, teq.AST) - end - end - - return nothing -end - - - -""" - (success, AST, equationInfo) = getSortedAndSolvedAST( - G, BLT, assign, A, B, stateSelectionFunctions; - log = false, logDetails=false, modelName = "???", - defaultParameterAndStartValues = nothing, - unitless=false) - -From the BLT on highest derivative level `(G, BLT, assign, A, B)`, -a set of callback functions `stateSelectionFunctions`, this function returns the -AST (Abstract Syntax Tree) of the sorted and solved equations as -a vector of `Expr` and information about the ODE equations -(especially names of the ODE states). - - -# Input arguments - -- `G`: G[i] is the Integer vector of variable indices of equation i. - A variable i can be a scalar or an array. - is of type Vector{ Vector{Int} } or Vector{ Any }. - -- `BLT::Vector{Vector{Int}}`: BLT[i] is the vector of equations belonging to BLT-block i. - In principal, **BLT[i] must include only highest derivative equations**. - In order to ease the generation of BLT[i], to not be forced to map indices forth and back - when calling the function that computes the BLT, the following assumption is - made: If BLT[j] contains only **one equation with B[j][1] > 0** then this BLT block is - a lower-derivative (dummy) block and is ignored. - -- `assign::Vector{Int}`: ei = assign[vi] is equation ei assigned to variable vi - -- `A`: A-Vector of Pantelides algorithm: - `A[i] = if der(v[i]) == v[k] then k else 0` - where `v[i]` is variable `i`. - -- `B`: B-Vector of Pantelides algorithm: - `B[i] = if der(e[i]) == e[k] then k else 0` - where `e[i]` is equation `i`. - -- `stateSelectionFunctions`: - An instance of immutable struct [`StateSelectionFunctions`](@ref) - that holds the information about the callback functions needed by - `getSortedAndSolvedAST`. - -- `log`: = true: Print debug information (to find bugs in the code). - -- `logDetails`: = true: Print detailed debug information, if `log = true`. - -- `modelName`: Name of model used in messages (otherwise not used). - - -# Output arguments - -A tuple with the following values: - -- `success::Bool`: = true, if generation was successful. = false, if - an error occured. Arguments `AST, equationInfo` hold the information that - have been collected until the error was detected. - -- `AST::Vector{Expr}`: The AST (Abstract Syntax Tree) of the sorted and solved equations - as Vector of `Expr`. - -- `equationInfo::`[`ModiaBase.EquationInfo`](@ref): Object that defines - the states of the ODE. - -When variable `v` has neither a `start` nor `init` attribute defined, -it is assumed to have `start=0.0`. - - -| defined | Description | -|:----------|:-----------------------------------------------------------| -| `start` | `start` is allowed to be changed during initialization | -| `init` | `init` is not allowed to be changed during initialization | - -Note, `start, init` influence the state selection. - -If an ODE state has neither a `start` nor an `init` value defined, a warning message is printed -(start/init value missing). - -If a solved variable has an `init` value, an information message is printed in case `log=true` -(init value has no effect). - - -# Sketch of the Algorithm - -Note, a statement *error xxx* means that an error occurs, -because it is not possible to transform the equations to an ODE. - - -## Generate equation sets - -From BLT blocks that are not further differentiated (= highest order BLT blocks), -generate the equation sets as described by the dummy derivative method of -[Mattsson and Söderlind (1992)](https://ieeexplore.ieee.org/document/274429). -See also sections 4.5 and 4.6 of the paper -[Otter and Elmqvist (2017)](https://modelica.org/events/modelica2017/proceedings/html/submissions/ecp17132565_OtterElmqvist.pdf). - - -## Define initial ODE states - -Initially, all differentiated variables are defined to be ODE states. -Whenever an initial ODE state is explicitly solved for in the course of the -algorithm, the variable is no longer an ODE state. - - -## Analyse equation sets - -For every equation set on every differentiation level perform the following actions: - -- If the equation set consists of *one* equation in *one* unknown, the equation - is solved for this unknown (*error, if this is not possible*). - -- If the equation set consists of *N linear* equations in *N* unknowns solve this - equation system (*error, if no linear system*):\\ - The equations are first teared to reduce the number of iteration variables and - afterwards the teared equation system is solved with a special iterator loop that - solves a linear equation system with an LU decomposition with column pivoting - (for details see [`ModiaBase.LinearEquationsIteration!`](@ref)).\\ - -- If the equation set consists of *N linear* equations in *M* unknowns (*M > N*) perform - the following actions (*error, if no linear system*).\\ - Note, due to the structure of BLT, this equation set cannot be - on highest derivative level; due to the Pantelides algorithm, - M < N is not possible). - - 1. Solve the equation system explicitly via tearing for N unknowns - (*error, if this is not possible*). Use start/init attributes of the M unknowns - to influence tearing (e.g. it is tried to eliminate first variables that have neither start not init values). - - 2. The N solved variables have been initially ODE states and are now defined to be - no ODE states (so called dummy states, that is algebraic variables). - -A BLT transformation of **all equations** is made under the assumption that -the selected ODE states are known. This phase is needed, because the dummy-derivative -ordering at the beginning does not necessarily provide already the right ordering. -In this phase the information from the previous phase is used (that is the already determined -tearing variables for systems of equations are utilized). -""" -function getSortedAndSolvedAST(Goriginal, # Typically ::Vector{Vector{Int}} - G, - BLT, # Typically ::Vector{Vector{Int}} - assign::Vector{Int}, - A::Vector{Int}, - B::Vector{Int}, - stateSelectionFunctions::StateSelectionFunctions; - log::Bool = false, - logDetails::Bool = false, - logStates::Bool = false, - modelName = "???", - unitless::Bool = false, - defaultParameterAndStartValues = nothing) - - if log - println("\n=== getSortedAndSolvedAST(...) started for $modelName.") - end - - # Define error/warning function - showMessage2(message; severity=0, details="",variables=Int[],equations=Int[]) = - stateSelectionFunctions.showMessage(message; from="getSortedAndSolvedAST for model $modelName", - severity=severity, details=details, variables=variables, equations=equations) - - # Construct data structure eq that holds all needed information for the state selection - # (this includes constructing the equation/constraint sets and initialization of tearing) - eq = EquationGraph(G, BLT, assign, A, B, stateSelectionFunctions, log, logDetails, showMessage2) - if !eq.success - @goto ERROR_RETURN - end - eq.equationInfo.defaultParameterAndStartValues = defaultParameterAndStartValues - - # Inspect every equation/constraint set in sequence from lowest-order to highest-order derivative - for j in eachindex(eq.eConstraintsVec) - # Equation/Constraint set and its unknowns that are inspected - eConstraints = eq.eConstraintsVec[j] - vConstraints = eq.vConstraintsVec[j] - - # Analyze all equation sets eConstraints[i] from lowest-order to highest-order derivatives - for i in eachindex(eConstraints) - if log - println("\n... Equation set $j.$i ..............................") - println("Equations: ") - printEquations(eq, eConstraints[i]) - println("Unknown variables: ") - printVariables(eq, vConstraints[i]) - end - - solveLinearEquation = true - if length(eConstraints[i]) == 1 && length(vConstraints[i]) == 1 - # One equation in one unknown - if i < length(eConstraints) - # Not on highest derivative level. Unknown is a dummy state. - v = vConstraints[i][1] - @assert(!eq.vActive[v]) - eq.vActive[v] = true - if log - println("One equation in one unknown variable: Unknown is a dummy state.") - end - else - if log - println("One equation in one unknown variable. Solve the equation:") - end - end - - # Solve equation for unknown v and add solved equation to AST - try - addSolvedEquations!(eq, eConstraints[i], vConstraints[i]) - solveLinearEquation = false - catch - if log - println("Not possible to solve the equation directly. Try so solve it as linear equation:") - end - solveLinearEquation = true - eq.fullAssignRev[eConstraints[i][1]] = -1 - end - end - - if solveLinearEquation - # N equations in M unknowns (M >= N) - @assert(length(vConstraints[i]) >= length(eConstraints[i])) - if log - N = length(eConstraints[i]) - M = length(vConstraints[i]) - println("$N equation(s) in $M unknown variable(s). Tear the system of equations:") - end - - # Tear systems of equations - if i == 1 - # On the lowest derivative level, initialize tearing with empty vectors - empty!(eq.der_eSolvedLower) - empty!(eq.der_vSolvedLower) - empty!(eq.der_vTearLower) - - # Variables/equations to be inspected by tearing - eq.vc = vConstraints[i] - eq.ec = eConstraints[i] - - else - # On derivative level i > 1, initialize tearing with the differentiated eSolved, vSolved, vTear - # from the previous level i-1. - # (eConstraints[i] is the derivative of eConstraints[i-1] + potentially additional equations) - eq.der_eSolvedLower = derivativeOfVlower(eq.eSolved, B) - eq.der_vSolvedLower = derivativeOfVlower(eq.vSolved, A) - eq.der_vTearLower = derivativeOfVlower(eq.vTear , A) - - # Variables/equations to be inspected by tearing (remove the parts that are known from level i-1) - eq.vc = setdiff(vConstraints[i], eq.der_vSolvedLower) - eq.ec = setdiff(eConstraints[i], eq.der_eSolvedLower) - end - - if log && i < length(eConstraints) - vWithStart = Int[] - for v in vConstraints[i] - (v_startOrInit,v_fixed) = eq.fc.var_startInitFixed(v) - if !(v_startOrInit isa Nothing) - push!(vWithStart, v) - end - end - println(" Unknowns with start or init values: ", getNames(eq, vWithStart)) - end - - tearEquationsWithCandidates!(eq) - if log - println(" Tearing variables: ", getNames(eq, eq.vTear)) - if length(eq.eResiduals) > 0 - println(" Residual equations: "); - printEquations(eq, eq.eResiduals; indent=8) - end - end - - # Select dummy states - if i < length(eConstraints) - # Not on highest derivative level. Solved variables are dummy states - for v in eq.vSolved - @assert(!eq.vActive[v]) - eq.vActive[v] = true - end - - if length(eConstraints[i]) == length(vConstraints[i]) - # N equations in N unknowns - tearing variables are dummy states - for v in eq.vTear - @assert(!eq.vActive[v]) - eq.vActive[v] = true - end - if log - println(" All unknowns are dummy states.") - end - - else - if log - println(" All solved unknowns are dummy states.") - end - end - elseif log - println(" All unknowns are solved.") - end - - # Further actions depending on type of equations - if length(eq.eResiduals) == 0 - # There are no residual equations, add solved equations to AST - addSolvedEquations!(eq, eq.eSolved, eq.vSolved) - - elseif length(eConstraints[i]) == length(vConstraints[i]) - # N equations in N unknowns with residual equations - - # Check that equation system is linear in the unknowns - (isLinear, hasConstantCoefficients) = isLinearEquationSystem!(eq, eConstraints[i], vConstraints[i]) - if !isLinear - if i == length(eConstraints) - # On highest derivative level: - # Assume that the equation system is linear, if at least one of the unknowns is a derivative - linearAssumption = true # temporarily - - #linearAssumption = false - #for v in vConstraints[i] - # if eq.Arev[v] > 0 - # linearAssumption = true - # end - #end - - if linearAssumption - isLinear = true - hasConstantCoefficients = false - #showMessage2("It is heuristically assumed that equation system is linear (although isLinearEquation returned isLinear=false)."; - # severity = WARNING, - # variables = vConstraints[i], - # equations = eConstraints[i]) - end - end - - if !isLinear - showMessage2("Cannot transform to ODE, because equation system is not linear."; - severity = ERROR, - variables = vConstraints[i], - equations = eConstraints[i]) - @goto ERROR_RETURN - end - end - - if log - println("Teared equation system is linear. Solve system with hasConstantCoefficients = $hasConstantCoefficients.") - end - - # Generate AST to build-up and solve linear equation system of the teared equations - addLinearEquations!(eq, hasConstantCoefficients, unitless) - - - else - # N equations in M unknowns with M > N and residual equations - @assert(length(vConstraints[i]) > length(eConstraints[i])) - - # Assertion, that this case is not possible on highest derivative level - @assert(i < length(eConstraints)) - - showMessage2("Cannot transform to ODE because constraint equations cannot be explicitly solved by tearing", - severity = ERROR, - variables = vConstraints[i], - equations = eConstraints[i]) - @goto ERROR_RETURN - end - end # if - end # for - end # for - - - # Determine ODE states - ODE_states = Int[] - for (ve,veActive) in enumerate(eq.vActive) - if !veActive - push!(ODE_states, ve) - end - end - - - # Store state information in equationInfo - x_info = eq.equationInfo.x_info - x_without_start = Int[] - x_vec = Int[] - x_vec_startOrInit = Any[] - x_vec_fixed = Bool[] - for v in ODE_states - (v_startOrInit, v_fixed) = eq.fc.var_startInitFixed(v) - if isnothing(v_startOrInit) - push!(x_without_start, v) - end - v_name = eq.fc.var_name(v) - - if isFixedLengthStartOrInit(v_startOrInit, v_name) - v_der_x = eq.A[v] - v_unit = unitless ? "" : eq.fc.var_unit(v) - v_julia_name = eq.fc.var_julia_name(v) - v_der_x_name = eq.fc.var_name(v_der_x) - v_der_x_julia_name = eq.fc.var_julia_name(v_der_x) - v_nominal = eq.fc.var_nominal(v) - v_unbounded = eq.fc.var_unbounded(v) - v_stateCategory = ModiaBase.XD - - push!(x_info, ModiaBase.StateElementInfo( - v_name, v_julia_name, - v_der_x_name, v_der_x_julia_name, - v_stateCategory, v_unit, v_startOrInit, v_fixed, - v_nominal, v_unbounded)) - else - # length(v_startOrInit) may change after compilation - push!(x_vec , v) - push!(x_vec_startOrInit, v_startOrInit) - push!(x_vec_fixed , v_fixed) - end - end - eq.equationInfo.nxFixedLength = length(x_info) - for (i,v) in enumerate(x_vec) - v_startOrInit = x_vec_startOrInit[i] - v_fixed = x_vec_fixed[i] - v_name = eq.fc.var_name(v) - v_der_x = eq.A[v] - v_unit = unitless ? "" : eq.fc.var_unit(v) - v_julia_name = eq.fc.var_julia_name(v) - v_der_x_name = eq.fc.var_name(v_der_x) - v_der_x_julia_name = eq.fc.var_julia_name(v_der_x) - v_nominal = eq.fc.var_nominal(v) - v_unbounded = eq.fc.var_unbounded(v) - v_stateCategory = ModiaBase.XD - - push!(x_info, ModiaBase.StateElementInfo( - v_name, v_julia_name, - v_der_x_name, v_der_x_julia_name, - v_stateCategory, v_unit, v_startOrInit, v_fixed, - v_nominal, v_unbounded)) - end - - # Handle systems with only algebraic variables, by introducing a dummy - # differential equation der_x[1] = -x[1]. - if length(ODE_states) == 0 - if log - println("Model has only algebraic variables.\n", - "Added a dummy differential equation der(_dummy_x) = -_dummy_x, _dummy_x(t0) = 0") - end - push!(eq.equationInfo.x_info, ModiaBase.StateElementInfo( - "_dummy_x", :(), "der(_dummy_x)", :(), XD, "", 0.0, true, NaN, false)) - end - - # Finalize equationInfo - initEquationInfo!(eq.equationInfo) - - # Print ODE states - if logStates - println("\nSelected ODE states: ") - x_table = ModiaBase.get_x_table(x_info) - show(stdout, x_table; allrows=true, allcols=true, summary=false, eltypes=false) - print("\n\n") - end - - # Sort equations and generate sorted AST (eq.AST) - sortEquations!(eq, Goriginal) - - # Check that all ODE states have a start value - if length(x_without_start) > 0 - showMessage2("Init/start values missing in the model for some ODE states."; - severity = WARNING, - variables = x_without_start) - end - - # Print info, if there are iteration variables without a start value - if length(eq.vTearWithoutStart) > 0 && !unitless - showMessage2("The following variables are iteration variables but have no start/init values defined.", - details = "If units are used in the model, start/init values with correct units should be defined\n" * - "to avoid unit errors during compilation.", - severity = INFORMATION, - variables = eq.vTearWithoutStart) - end - - # Print warning, if there are variables with fixed=true that are - # explicitly solved for and log=true - if log && length(eq.equationInfo.vSolvedWithFixedTrue) > 0 - showMessage2("The following variables have an 'init' initialization and are explicitly solved for.", - details = "Therefore, the 'init' values have no effect, but must exactly match the values,\n"* - "computed during initialization. Otherwise this gives a run-time error.\n"* - "It is adviced to use 'start' initialization or remove initialization for these variables." , - severity = WARNING, - variables = eq.vSolvedWithFixedTrue) - end - - - if log - println("\n=== getSortedAndSolvedAST(...) terminated for $modelName.\n") - end - return (true, eq.AST, eq.equationInfo) - - @label ERROR_RETURN - return (false, eq.AST, eq.equationInfo) -end diff --git a/test/Project.toml b/test/Project.toml deleted file mode 100644 index 0c36332..0000000 --- a/test/Project.toml +++ /dev/null @@ -1,2 +0,0 @@ -[deps] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/TestLinearEquations.jl b/test/TestLinearEquations.jl deleted file mode 100644 index 2920f22..0000000 --- a/test/TestLinearEquations.jl +++ /dev/null @@ -1,95 +0,0 @@ -""" - module TestLinearEquations - Test LinearEquations(..) in ModiaBase/src/EquationAndStateInfo.jl -""" -module TestLinearEquations - -using Test -using ModiaBase -using ModiaBase.TimerOutputs - -println("... Test TestLinearEquations") - -# Solve A*x = b as residual equation -const A = [1.0 2.0 3.0 4.0 5.0; - 0.0 2.0 0.0 0.0 0.0; - 0.0 0.0 3.0 0.0 0.0; - 0.0 0.0 0.0 4.0 0.0; - 0.0 0.0 0.0 0.0 5.0] -const b = [-1.0, -2.0, -3.0, -4.0, -5.0] -const isInitial = true # the first call must be performed with isInitial=true -const time = 0.0 # only for warning/error messages -const timer = TimerOutput() # timer to measure the time to solve the linear equation system -const leq = ModiaBase.LinearEquations(["x2", "x1"], [:x1], [1,5], 1, false) - -computeResidual1(x1) = A*x1-b -computeResidual2(x1,x2) = x1[1] + 2*x1[2] + 3*x1[3] + 4*x2 - -# Solution -const x1_sol = A\b -const x2_sol = -(x1_sol[1] + 2*x1_sol[2] + 3*x1_sol[3])/4.0 - -# Test function for ODE mode -function solveLinearEquation(leq,isInitial,time,timer) - local x1::Vector{Float64}, x2 - leq.mode = -3 - while ModiaBase.LinearEquationsIteration!(leq, isInitial, time, timer) - x1 = leq.x_vec[1] - x2 = leq.x[1] - append!(leq.residuals, computeResidual1(x1)) - append!(leq.residuals, computeResidual2(x1,x2)) - end - return (x1,x2) -end - -# Test function for DAE mode -function solveLinearEquation(leq,isInitial,solve_leq,isStoreResult,time,timer) - local x1, x2 - leq.mode = -3 - while ModiaBase.LinearEquationsIteration!(leq, isInitial, solve_leq, isStoreResult, time, timer) - x1 = leq.x_vec[1] - x2 = leq.x[1] - append!(leq.residuals, computeResidual1(x1)) - append!(leq.residuals, computeResidual2(x1,x2)) - end - return (x1,x2) -end - -@testset "Test TestLinearEquations" begin - # Test ODE mode (x is always explicitly computed by solving the linear equation system) - leq.odeMode = true - (x1,x2) = solveLinearEquation(leq,isInitial,time,timer) - - d1 = x1_sol-x1 - err1 = sqrt(d1'*d1) + abs(x2_sol - x2) - @test isapprox(err1, 0.0, atol=1e-13) - - - # Test DAE mode - leq.odeMode = false - - # x is computed explicitly at initialization or at an event - (x1,x2) = solveLinearEquation(leq,false,true,false,time,timer) - d2 = x1_sol-x1 - err2 = sqrt(d2'*d2) + abs(x2_sol - x2) - @test isapprox(err2, 0.0, atol=1e-13) - - # x is computed by DAE solver - leq.x = [6.0, 1.0, 2.0, 3.0, 4.0, 5.0] # Provide x from DAE solver - (x1,x2) = solveLinearEquation(leq,false,false,false,time,timer) - residuals = leq.residuals # Store residuals in DAE solver - d3 = leq.x - vcat(x2,x1) - d4 = vcat(A*leq.x[2:6] - b, leq.x[2] + 2*leq.x[3] + 3*leq.x[4] + 4*leq.x[1]) - residuals - err3 = sqrt(d3'*d3) - err4 = sqrt(d4'*d4) - @test isapprox(err3, 0.0, atol=1e-15) - @test isapprox(err4, 0.0, atol=1e-15) - - # x is computed by DAE solver at a communication point (leq.residuals is not provided) - leq.x = [6.0, 1.0, 2.0, 3.0, 4.0, 5.0] # Provide x from DAE solver - (x1,x2) = solveLinearEquation(leq,false,false,true,time,timer) - d5 = leq.x - vcat(x2,x1) - err5 = sqrt(d5'*d5) - @test isapprox(err5, 0.0, atol=1e-15) -end - -end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 34d4070..f6136cd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,7 +9,6 @@ using Test include("TestDifferentiate.jl") include("TestTearing.jl") include("TestLinearIntegerEquations.jl") - include("TestLinearEquations.jl") end end \ No newline at end of file