diff --git a/Manifest.toml b/Manifest.toml
index 94ceb1aefd9d24264a211efed335b80b954768db..29b2bbfdf36df5445385c51d504ae12052560b6a 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -1,14 +1,18 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.9.3"
+julia_version = "1.10.0-rc2"
 manifest_format = "2.0"
-project_hash = "76d77479c493e8b2f7f5c60b9dc879fe8436281c"
+project_hash = "4772e5e6caf338d2c755b5d4461c2d0c5c7482c5"
 
 [[deps.Adapt]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "195c5505521008abea5aee4f96930717958eac6f"
+deps = ["LinearAlgebra", "Requires"]
+git-tree-sha1 = "cde29ddf7e5726c9fb511f340244ea3481267608"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.4.0"
+version = "3.7.2"
+weakdeps = ["StaticArrays"]
+
+    [deps.Adapt.extensions]
+    AdaptStaticArraysExt = "StaticArrays"
 
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
@@ -16,9 +20,9 @@ version = "1.1.1"
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "f83ec24f76d4c8f525099b2ac475fc098138ec31"
+git-tree-sha1 = "247efbccf92448be332d154d6ca56b9fcdd93c31"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.4.11"
+version = "7.6.1"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
@@ -36,12 +40,6 @@ version = "7.4.11"
     StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 
-[[deps.ArrayInterfaceCore]]
-deps = ["LinearAlgebra", "SnoopPrecompile", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "e5f08b5689b1aad068e01751889f2f615c7db36d"
-uuid = "30b0a656-2188-435a-8636-2ec0e6a096e2"
-version = "0.1.29"
-
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
@@ -56,21 +54,19 @@ version = "0.1.5"
 
 [[deps.CPUSummary]]
 deps = ["CpuId", "IfElse", "PrecompileTools", "Static"]
-git-tree-sha1 = "89e0654ed8c7aebad6d5ad235d6242c2d737a928"
+git-tree-sha1 = "601f7e7b3d36f18790e2caf83a882d88e9b71ff1"
 uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
-version = "0.2.3"
+version = "0.2.4"
 
 [[deps.ChainRulesCore]]
-deps = ["Compat", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "e7ff6cadf743c098e08fca25c91103ee4303c9bb"
+deps = ["Compat", "LinearAlgebra"]
+git-tree-sha1 = "e0af648f0692ec1691b5d094b8724ba1346281cf"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.15.6"
+version = "1.18.0"
+weakdeps = ["SparseArrays"]
 
-[[deps.ChangesOfVariables]]
-deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
-git-tree-sha1 = "38f7a08f19d8810338d4f5085211c7dfa5d5bdd8"
-uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
-version = "0.1.4"
+    [deps.ChainRulesCore.extensions]
+    ChainRulesCoreSparseArraysExt = "SparseArrays"
 
 [[deps.CloseOpenIntervals]]
 deps = ["Static", "StaticArrayInterface"]
@@ -80,31 +76,43 @@ version = "0.1.12"
 
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
-git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
+git-tree-sha1 = "cd67fc487743b2f0fd4380d4cbd3a24660d0eec8"
 uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
-version = "0.7.0"
+version = "0.7.3"
 
 [[deps.CommonSolve]]
-git-tree-sha1 = "9441451ee712d1aec22edad62db1a9af3dc8d852"
+git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c"
 uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
-version = "0.2.3"
+version = "0.2.4"
 
 [[deps.Compat]]
-deps = ["Dates", "LinearAlgebra", "UUIDs"]
-git-tree-sha1 = "00a2cccc7f098ff3b66806862d275ca3db9e6e5a"
+deps = ["UUIDs"]
+git-tree-sha1 = "886826d76ea9e72b35fcd000e535588f7b60f21d"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.5.0"
+version = "4.10.1"
+weakdeps = ["Dates", "LinearAlgebra"]
+
+    [deps.Compat.extensions]
+    CompatLinearAlgebraExt = "LinearAlgebra"
 
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+0"
+version = "1.0.5+1"
 
 [[deps.ConstructionBase]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "fb21ddd70a051d882a1686a5a550990bbe371a95"
+git-tree-sha1 = "c53fc348ca4d40d7b371e71fd52251839080cbc9"
 uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
-version = "1.4.1"
+version = "1.5.4"
+
+    [deps.ConstructionBase.extensions]
+    ConstructionBaseIntervalSetsExt = "IntervalSets"
+    ConstructionBaseStaticArraysExt = "StaticArrays"
+
+    [deps.ConstructionBase.weakdeps]
+    IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [[deps.CpuId]]
 deps = ["Markdown"]
@@ -118,9 +126,9 @@ uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
 version = "4.1.1"
 
 [[deps.DataAPI]]
-git-tree-sha1 = "e8119c1a33d267e16108be441a287a6981ba1630"
+git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c"
 uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
-version = "1.14.0"
+version = "1.15.0"
 
 [[deps.DataValueInterfaces]]
 git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
@@ -147,24 +155,34 @@ uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 version = "1.6.0"
 
 [[deps.ExprTools]]
-git-tree-sha1 = "56559bbef6ca5ea0c0818fa5c90320398a6fbf8d"
+git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec"
 uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
-version = "0.1.8"
+version = "0.1.10"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
-git-tree-sha1 = "7be5f99f7d15578798f338f5433b6c432ea8037b"
+git-tree-sha1 = "299dc33549f68299137e51e6d49a13b5b1da9673"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.16.0"
+version = "1.16.1"
 
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
 [[deps.FillArrays]]
-deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
-git-tree-sha1 = "e17cc4dc2d0b0b568e80d937de8ed8341822de67"
+deps = ["LinearAlgebra", "Random"]
+git-tree-sha1 = "25a10f2b86118664293062705fd9c7e2eda881a2"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.2.0"
+version = "1.9.2"
+
+    [deps.FillArrays.extensions]
+    FillArraysPDMatsExt = "PDMats"
+    FillArraysSparseArraysExt = "SparseArrays"
+    FillArraysStatisticsExt = "Statistics"
+
+    [deps.FillArrays.weakdeps]
+    PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.Formatting]]
 deps = ["Printf"]
@@ -178,21 +196,27 @@ uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
 
 [[deps.GPUArraysCore]]
 deps = ["Adapt"]
-git-tree-sha1 = "6872f5ec8fd1a38880f027a26739d42dcda6691f"
+git-tree-sha1 = "2d6ca471a6c7b536127afccfa7564b5b39227fe0"
 uuid = "46192b85-c4d5-4398-a991-12ede77f4527"
-version = "0.1.2"
+version = "0.1.5"
 
 [[deps.HDF5]]
-deps = ["Compat", "HDF5_jll", "Libdl", "Mmap", "Random", "Requires", "UUIDs"]
-git-tree-sha1 = "b5df7c3cab3a00c33c2e09c6bd23982a75e2fbb2"
+deps = ["Compat", "HDF5_jll", "Libdl", "MPIPreferences", "Mmap", "Preferences", "Printf", "Random", "Requires", "UUIDs"]
+git-tree-sha1 = "26407bd1c60129062cec9da63dc7d08251544d53"
 uuid = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
-version = "0.16.13"
+version = "0.17.1"
+
+    [deps.HDF5.extensions]
+    MPIExt = "MPI"
+
+    [deps.HDF5.weakdeps]
+    MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 
 [[deps.HDF5_jll]]
-deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenSSL_jll", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "4cc2bb72df6ff40b055295fdef6d92955f9dede8"
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"]
+git-tree-sha1 = "38c8874692d48d5440d5752d6c74b0c6b0b60739"
 uuid = "0234f1f7-429e-5d53-9886-15a909be8d59"
-version = "1.12.2+2"
+version = "1.14.2+1"
 
 [[deps.HostCPUFeatures]]
 deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"]
@@ -200,6 +224,12 @@ git-tree-sha1 = "eb8fed28f4994600e29beef49744639d985a04b2"
 uuid = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0"
 version = "0.1.16"
 
+[[deps.Hwloc_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "ca0f6bf568b4bfc807e7537f081c81e35ceca114"
+uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
+version = "2.10.0+0"
+
 [[deps.IfElse]]
 git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
 uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
@@ -207,9 +237,9 @@ version = "0.1.1"
 
 [[deps.InlineStrings]]
 deps = ["Parsers"]
-git-tree-sha1 = "0cf92ec945125946352f3d46c96976ab972bde6f"
+git-tree-sha1 = "9cc2baf75c6d09f9da536ddf58eb2f29dedaf461"
 uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
-version = "1.3.2"
+version = "1.4.0"
 
 [[deps.InteractiveUtils]]
 deps = ["Markdown"]
@@ -217,20 +247,14 @@ uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [[deps.Intervals]]
 deps = ["Dates", "Printf", "RecipesBase", "Serialization", "TimeZones"]
-git-tree-sha1 = "f3c7f871d642d244e7a27e3fb81e8441e13230d8"
+git-tree-sha1 = "ac0aaa807ed5eaf13f67afe188ebc07e828ff640"
 uuid = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
-version = "1.8.0"
-
-[[deps.InverseFunctions]]
-deps = ["Test"]
-git-tree-sha1 = "49510dfcb407e572524ba94aeae2fced1f3feb0f"
-uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
-version = "0.1.8"
+version = "1.10.0"
 
 [[deps.IrrationalConstants]]
-git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
+git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
 uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
-version = "0.1.1"
+version = "0.2.2"
 
 [[deps.IteratorInterfaceExtensions]]
 git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
@@ -238,27 +262,33 @@ uuid = "82899510-4779-5014-852e-03e436cf321d"
 version = "1.0.0"
 
 [[deps.JLLWrappers]]
-deps = ["Preferences"]
-git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
+deps = ["Artifacts", "Preferences"]
+git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.4.1"
+version = "1.5.0"
 
 [[deps.Jacobi]]
 deps = ["Polynomials", "SpecialFunctions"]
-git-tree-sha1 = "c793da50a68cb78f33f77c19c293a4e8398d86d5"
+git-tree-sha1 = "9601b0f23e66ca61f0d48d59054f161c6c8f6f33"
 uuid = "83f21c0b-4282-5fbc-9e3f-f6da3d2e584c"
-version = "0.5.1"
+version = "0.6.0"
+
+[[deps.LLVMOpenMP_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "f689897ccbe049adb19a065c495e75f372ecd42b"
+uuid = "1d63c593-3942-5779-bab2-d838dc0a180e"
+version = "15.0.4+0"
 
 [[deps.LaTeXStrings]]
-git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
+git-tree-sha1 = "50901ebc375ed41dbf8058da26f9de442febbbec"
 uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
-version = "1.3.0"
+version = "1.3.1"
 
 [[deps.LayoutPointers]]
 deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface"]
-git-tree-sha1 = "88b8f66b604da079a627b6fb2860d3704a6729a1"
+git-tree-sha1 = "62edfee3211981241b57ff1cedf4d74d79519277"
 uuid = "10f19ff3-798f-405d-979b-55457f8fc047"
-version = "0.1.14"
+version = "0.1.15"
 
 [[deps.LazyArtifacts]]
 deps = ["Artifacts", "Pkg"]
@@ -267,55 +297,70 @@ uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
-version = "0.6.3"
+version = "0.6.4"
 
 [[deps.LibCURL_jll]]
 deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
 uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
-version = "7.84.0+0"
+version = "8.4.0+0"
 
 [[deps.LibGit2]]
-deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
+[[deps.LibGit2_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"]
+uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5"
+version = "1.6.4+0"
+
 [[deps.LibSSH2_jll]]
 deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
 uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
-version = "1.10.2+0"
+version = "1.11.0+1"
 
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
 [[deps.Libiconv_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "c7cb1f5d892775ba13767a87c7ada0b980ea0a71"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175"
 uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
-version = "1.16.1+2"
+version = "1.17.0+0"
 
 [[deps.LightXML]]
 deps = ["Libdl", "XML2_jll"]
-git-tree-sha1 = "e129d9391168c677cd4800f5c0abb1ed8cb3794f"
+git-tree-sha1 = "3a994404d3f6709610701c7dabfc03fed87a81f8"
 uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
-version = "0.9.0"
+version = "0.9.1"
 
 [[deps.LinearAlgebra]]
 deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[deps.LogExpFunctions]]
-deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "946607f84feb96220f480e0422d3484c49c00239"
+deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.19"
+version = "0.3.26"
+
+    [deps.LogExpFunctions.extensions]
+    LogExpFunctionsChainRulesCoreExt = "ChainRulesCore"
+    LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables"
+    LogExpFunctionsInverseFunctionsExt = "InverseFunctions"
+
+    [deps.LogExpFunctions.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+    InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 
 [[deps.Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
 [[deps.LoopVectorization]]
-deps = ["ArrayInterface", "ArrayInterfaceCore", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"]
-git-tree-sha1 = "c88a4afe1703d731b1c4fdf4e3c7e77e3b176ea2"
+deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"]
+git-tree-sha1 = "0f5648fbae0d015e3abe5867bca2b362f67a5894"
 uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
-version = "0.12.165"
+version = "0.12.166"
 
     [deps.LoopVectorization.extensions]
     ForwardDiffExt = ["ChainRulesCore", "ForwardDiff"]
@@ -326,11 +371,29 @@ version = "0.12.165"
     ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
     SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
+[[deps.MPICH_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
+git-tree-sha1 = "8a5b4d2220377d1ece13f49438d71ad20cf1ba83"
+uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4"
+version = "4.1.2+0"
+
+[[deps.MPIPreferences]]
+deps = ["Libdl", "Preferences"]
+git-tree-sha1 = "8f6af051b9e8ec597fa09d8885ed79fd582f33c9"
+uuid = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"
+version = "0.1.10"
+
+[[deps.MPItrampoline_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
+git-tree-sha1 = "6979eccb6a9edbbb62681e158443e79ecc0d056a"
+uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748"
+version = "5.3.1+0"
+
 [[deps.MacroTools]]
 deps = ["Markdown", "Random"]
-git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2"
+git-tree-sha1 = "9ee1618cbf5240e6d4e0371d6f24065083f60c48"
 uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-version = "0.5.10"
+version = "0.5.11"
 
 [[deps.ManualMemory]]
 git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd"
@@ -344,26 +407,32 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
-version = "2.28.2+0"
+version = "2.28.2+1"
+
+[[deps.MicrosoftMPI_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "b01beb91d20b0d1312a9471a36017b5b339d26de"
+uuid = "9237b28f-5490-5468-be7b-bb81f5f5e6cf"
+version = "10.1.4+1"
 
 [[deps.Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
 [[deps.Mocking]]
 deps = ["Compat", "ExprTools"]
-git-tree-sha1 = "c272302b22479a24d1cf48c114ad702933414f80"
+git-tree-sha1 = "4cc0c5a83933648b615c36c2b956d94fda70641e"
 uuid = "78c3b35d-d492-501b-9361-3d52fe80e533"
-version = "0.7.5"
+version = "0.7.7"
 
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
-version = "2022.10.11"
+version = "2023.1.10"
 
 [[deps.MutableArithmetics]]
 deps = ["LinearAlgebra", "SparseArrays", "Test"]
-git-tree-sha1 = "aa532179d4a643d4bd9f328589ca01fa20a0d197"
+git-tree-sha1 = "806eea990fb41f9b36f1253e5697aa645bf6a9f8"
 uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
-version = "1.1.0"
+version = "1.4.0"
 
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
@@ -371,25 +440,31 @@ version = "1.2.0"
 
 [[deps.OffsetArrays]]
 deps = ["Adapt"]
-git-tree-sha1 = "f71d8950b724e9ff6110fc948dff5a329f901d64"
+git-tree-sha1 = "2ac17d29c523ce1cd38e27785a7d23024853a4bb"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.12.8"
+version = "1.12.10"
 
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.21+4"
+version = "0.3.23+2"
 
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.1+0"
+version = "0.8.1+2"
+
+[[deps.OpenMPI_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "PMIx_jll", "TOML", "Zlib_jll", "libevent_jll", "prrte_jll"]
+git-tree-sha1 = "694458ae803b684f09c07f90459cb79655fb377d"
+uuid = "fe0851c0-eecd-5654-98d4-656369965a5c"
+version = "5.0.0+0"
 
 [[deps.OpenSSL_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "f6e9dba33f9f2c44e08a020b0caf6903be540004"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "cc6e1927ac521b659af340e0ca45828a3ffc748f"
 uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
-version = "1.1.19+0"
+version = "3.0.12+0"
 
 [[deps.OpenSpecFun_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
@@ -398,9 +473,15 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 version = "0.5.5+0"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.4.1"
+version = "1.6.3"
+
+[[deps.PMIx_jll]]
+deps = ["Artifacts", "Hwloc_jll", "JLLWrappers", "Libdl", "Zlib_jll", "libevent_jll"]
+git-tree-sha1 = "8b3b19351fa24791f94d7ae85faf845ca1362541"
+uuid = "32165bc3-0280-59bc-8c0b-c33b6203efab"
+version = "4.2.7+0"
 
 [[deps.Parameters]]
 deps = ["OrderedCollections", "UnPack"]
@@ -409,15 +490,15 @@ uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 version = "0.12.3"
 
 [[deps.Parsers]]
-deps = ["Dates", "SnoopPrecompile"]
-git-tree-sha1 = "6466e524967496866901a78fca3f2e9ea445a559"
+deps = ["Dates", "PrecompileTools", "UUIDs"]
+git-tree-sha1 = "a935806434c9d4c506ba941871b327b96d41f2bf"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.5.2"
+version = "2.8.0"
 
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-version = "1.9.2"
+version = "1.10.0"
 
 [[deps.PolyesterWeave]]
 deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"]
@@ -439,9 +520,9 @@ version = "1.2.0"
 
 [[deps.Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
+git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.3.0"
+version = "1.4.1"
 
 [[deps.PrettyTables]]
 deps = ["Crayons", "Formatting", "LaTeXStrings", "Markdown", "Reexport", "StringManipulation", "Tables"]
@@ -455,23 +536,23 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [[deps.ProgressMeter]]
 deps = ["Distributed", "Printf"]
-git-tree-sha1 = "d7a7aef8f8f2d537104f170139553b14dfe39fe9"
+git-tree-sha1 = "00099623ffee15972c16111bcf84c58a0051257c"
 uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
-version = "1.7.2"
+version = "1.9.0"
 
 [[deps.REPL]]
 deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 
 [[deps.Random]]
-deps = ["SHA", "Serialization"]
+deps = ["SHA"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[deps.RecipesBase]]
-deps = ["SnoopPrecompile"]
-git-tree-sha1 = "18c35ed630d7229c5584b945641a73ca83fb5213"
+deps = ["PrecompileTools"]
+git-tree-sha1 = "5c3d09cc4f31f5fc6af001c250bf1278733100ff"
 uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-version = "1.3.2"
+version = "1.3.4"
 
 [[deps.Reexport]]
 git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
@@ -486,9 +567,21 @@ version = "1.3.0"
 
 [[deps.Roots]]
 deps = ["ChainRulesCore", "CommonSolve", "Printf", "Setfield"]
-git-tree-sha1 = "a3db467ce768343235032a1ca0830fc64158dadf"
+git-tree-sha1 = "0f1d92463a020321983d04c110f476c274bafe2e"
 uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-version = "2.0.8"
+version = "2.0.22"
+
+    [deps.Roots.extensions]
+    RootsForwardDiffExt = "ForwardDiff"
+    RootsIntervalRootFindingExt = "IntervalRootFinding"
+    RootsSymPyExt = "SymPy"
+    RootsSymPyPythonCallExt = "SymPyPythonCall"
+
+    [deps.Roots.weakdeps]
+    ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+    IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
+    SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
+    SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c"
 
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
@@ -501,15 +594,15 @@ version = "0.1.0"
 
 [[deps.SLEEFPirates]]
 deps = ["IfElse", "Static", "VectorizationBase"]
-git-tree-sha1 = "4b8586aece42bee682399c4c4aee95446aa5cd19"
+git-tree-sha1 = "3aac6d68c5e57449f5b9b865c9ba50ac2970c4cf"
 uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
-version = "0.6.39"
+version = "0.6.42"
 
 [[deps.Scratch]]
 deps = ["Dates"]
-git-tree-sha1 = "f94f779c94e58bf9ea243e77a37e16d9de9126bd"
+git-tree-sha1 = "3bac05bc7e74a75fd9cba4295cde4045d9fe2386"
 uuid = "6c6a2e73-6563-6170-7368-637461726353"
-version = "1.1.1"
+version = "1.2.1"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -520,23 +613,23 @@ git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac"
 uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 version = "1.1.1"
 
-[[deps.SnoopPrecompile]]
-git-tree-sha1 = "f604441450a3c0569830946e5b33b78c928e1a85"
-uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c"
-version = "1.0.1"
-
 [[deps.Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
 [[deps.SparseArrays]]
 deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+version = "1.10.0"
 
 [[deps.SpecialFunctions]]
-deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "5d65101b2ed17a8862c4c05639c3ddc7f3d791e1"
+deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
+git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "1.8.7"
+version = "2.3.1"
+weakdeps = ["ChainRulesCore"]
+
+    [deps.SpecialFunctions.extensions]
+    SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
 
 [[deps.Static]]
 deps = ["IfElse"]
@@ -545,10 +638,10 @@ uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 version = "0.8.8"
 
 [[deps.StaticArrayInterface]]
-deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "Requires", "SnoopPrecompile", "SparseArrays", "Static", "SuiteSparse"]
-git-tree-sha1 = "33040351d2403b84afce74dae2e22d3f5b18edcb"
+deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"]
+git-tree-sha1 = "03fec6800a986d191f64f5c0996b59ed526eda25"
 uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718"
-version = "1.4.0"
+version = "1.4.1"
 weakdeps = ["OffsetArrays", "StaticArrays"]
 
     [deps.StaticArrayInterface.extensions]
@@ -556,31 +649,33 @@ weakdeps = ["OffsetArrays", "StaticArrays"]
     StaticArrayInterfaceStaticArraysExt = "StaticArrays"
 
 [[deps.StaticArrays]]
-deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"]
-git-tree-sha1 = "6954a456979f23d05085727adb17c4551c19ecd1"
+deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
+git-tree-sha1 = "5ef59aea6f18c25168842bded46b16662141ab87"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.5.12"
+version = "1.7.0"
+
+    [deps.StaticArrays.extensions]
+    StaticArraysStatisticsExt = "Statistics"
+
+    [deps.StaticArrays.weakdeps]
+    Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.StaticArraysCore]]
-git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a"
+git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
 uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-version = "1.4.0"
-
-[[deps.Statistics]]
-deps = ["LinearAlgebra", "SparseArrays"]
-uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-version = "1.9.0"
+version = "1.4.2"
 
 [[deps.StringManipulation]]
-git-tree-sha1 = "46da2434b41f41ac3594ee9816ce5541c6096123"
+deps = ["PrecompileTools"]
+git-tree-sha1 = "a04cabe79c5f01f4d723cc6704070ada0b9d46d5"
 uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e"
-version = "0.3.0"
+version = "0.3.4"
 
 [[deps.StructArrays]]
-deps = ["Adapt", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"]
-git-tree-sha1 = "b03a3b745aa49b566f128977a7dd1be8711c5e71"
+deps = ["Adapt", "ConstructionBase", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"]
+git-tree-sha1 = "0a3db38e4cce3c54fe7a71f831cd7b6194a54213"
 uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
-version = "0.6.14"
+version = "0.6.16"
 
 [[deps.SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
@@ -589,13 +684,19 @@ uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 [[deps.SuiteSparse_jll]]
 deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"]
 uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
-version = "5.10.1+6"
+version = "7.2.1+1"
 
 [[deps.TOML]]
 deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 version = "1.0.3"
 
+[[deps.TZJData]]
+deps = ["Artifacts"]
+git-tree-sha1 = "d39314cdbaf5b90a047db33858626f8d1cc973e1"
+uuid = "dc5dba14-91b3-4cab-a142-028a31da12f7"
+version = "1.0.0+2023c"
+
 [[deps.TableTraits]]
 deps = ["IteratorInterfaceExtensions"]
 git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
@@ -603,10 +704,10 @@ uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
 version = "1.0.1"
 
 [[deps.Tables]]
-deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
-git-tree-sha1 = "c79322d36826aa2f4fd8ecfa96ddb47b174ac78d"
+deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits"]
+git-tree-sha1 = "cb76cf677714c095e535e3501ac7954732aeea2d"
 uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "1.10.0"
+version = "1.11.1"
 
 [[deps.Tar]]
 deps = ["ArgTools", "SHA"]
@@ -624,16 +725,20 @@ uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5"
 version = "0.5.2"
 
 [[deps.TimeZones]]
-deps = ["Dates", "Downloads", "InlineStrings", "LazyArtifacts", "Mocking", "Printf", "RecipesBase", "Scratch", "Unicode"]
-git-tree-sha1 = "a92ec4466fc6e3dd704e2668b5e7f24add36d242"
+deps = ["Artifacts", "Dates", "Downloads", "InlineStrings", "LazyArtifacts", "Mocking", "Printf", "Scratch", "TZJData", "Unicode", "p7zip_jll"]
+git-tree-sha1 = "89e64d61ef3cd9e80f7fc12b7d13db2d75a23c03"
 uuid = "f269a46b-ccf7-5d73-abea-4c690281aa53"
-version = "1.9.1"
+version = "1.13.0"
+weakdeps = ["RecipesBase"]
+
+    [deps.TimeZones.extensions]
+    TimeZonesRecipesBaseExt = "RecipesBase"
 
 [[deps.TimerOutputs]]
 deps = ["ExprTools", "Printf"]
-git-tree-sha1 = "f2fd3f288dfc6f507b0c3a2eb3bac009251e548b"
+git-tree-sha1 = "f548a9e9c490030e545f72074a41edfd0e5bcdd7"
 uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
-version = "0.5.22"
+version = "0.5.23"
 
 [[deps.ToggleableAsserts]]
 git-tree-sha1 = "7e70bd2062df35258220ec30f454be1056efb431"
@@ -642,9 +747,9 @@ version = "0.1.0"
 
 [[deps.TranscodingStreams]]
 deps = ["Random", "Test"]
-git-tree-sha1 = "e4bdc63f5c6d62e80eb1c0043fcc0360d5950ff7"
+git-tree-sha1 = "9a6ae7ed916312b41236fcef7e0af564ef934769"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
-version = "0.9.10"
+version = "0.9.13"
 
 [[deps.UUIDs]]
 deps = ["Random", "SHA"]
@@ -665,38 +770,56 @@ version = "1.0.1"
 
 [[deps.VectorizationBase]]
 deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static", "StaticArrayInterface"]
-git-tree-sha1 = "b182207d4af54ac64cbc71797765068fdeff475d"
+git-tree-sha1 = "7209df901e6ed7489fe9b7aa3e46fb788e15db85"
 uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
-version = "0.21.64"
+version = "0.21.65"
 
 [[deps.WriteVTK]]
 deps = ["Base64", "CodecZlib", "FillArrays", "LightXML", "TranscodingStreams", "VTKBase"]
-git-tree-sha1 = "7b46936613e41cfe1c6a5897d243ddcab8feabec"
+git-tree-sha1 = "41f0dc2a8f6fd860c266b91fd5cdf4fead65ae69"
 uuid = "64499a7a-5c06-52f2-abe2-ccb03c286192"
-version = "1.18.0"
+version = "1.18.1"
 
 [[deps.XML2_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "93c41695bc1c08c46c5899f4fe06d6ead504bb73"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
+git-tree-sha1 = "801cbe47eae69adc50f36c3caec4758d2650741b"
 uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
-version = "2.10.3+0"
+version = "2.12.2+0"
 
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
-version = "1.2.13+0"
+version = "1.2.13+1"
+
+[[deps.libaec_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "eddd19a8dea6b139ea97bdc8a0e2667d4b661720"
+uuid = "477f73a3-ac25-53e9-8cc3-50b2fa2566f0"
+version = "1.0.6+1"
 
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.8.0+0"
+version = "5.8.0+1"
+
+[[deps.libevent_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "OpenSSL_jll"]
+git-tree-sha1 = "f04ec6d9a186115fb38f858f05c0c4e1b7fc9dcb"
+uuid = "1080aeaf-3a6a-583e-a51c-c537b09f60ec"
+version = "2.1.13+1"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
-version = "1.48.0+0"
+version = "1.52.0+1"
 
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
-version = "17.4.0+0"
+version = "17.4.0+2"
+
+[[deps.prrte_jll]]
+deps = ["Artifacts", "Hwloc_jll", "JLLWrappers", "Libdl", "PMIx_jll", "libevent_jll"]
+git-tree-sha1 = "5adb2d7a18a30280feb66cad6f1a1dfdca2dc7b0"
+uuid = "eb928a42-fffd-568d-ab9c-3f5d54fc65b9"
+version = "3.0.2+0"
diff --git a/Project.toml b/Project.toml
index f56d3a01e04d9fcada0f16b7ce7ac21a50464402..c1a04c28fb8806825c7243ee2f313e5f10eb63b9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -28,6 +28,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 julia = "1.9"
+PrettyTables = "<2.3.0"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/examples/euler_isentropic_flow/mono.toml b/examples/euler_isentropic_flow/mono.toml
new file mode 100644
index 0000000000000000000000000000000000000000..cf957f20cb2541368fa6ef988519f3065886266e
--- /dev/null
+++ b/examples/euler_isentropic_flow/mono.toml
@@ -0,0 +1,26 @@
+[EquationOfState]
+idealgas_gamma = 3.0
+eos = "idealgas"
+
+[EulerEq]
+id = "smooth_isentropic_flow"
+av_regularization = "mono"
+
+[Evolution]
+cfl = 0.75
+tend = 0.50
+
+[Mesh]
+range = [-1.0, 1.0]
+n = 5
+k = 680
+basis = "lgl"
+
+[Output]
+variables = ["rho", "q", "E"]
+aligned_ts = "$(collect(range(0.01,0.50,step=0.01)))"
+
+[HRSC]
+method = "av"
+av_method = "entropy"
+entropy_ce = 20.0
diff --git a/examples/euler_isentropic_flow/navier.toml b/examples/euler_isentropic_flow/navier.toml
new file mode 100644
index 0000000000000000000000000000000000000000..94e21dcc44726b79037a7786510f3a3c0539589e
--- /dev/null
+++ b/examples/euler_isentropic_flow/navier.toml
@@ -0,0 +1,26 @@
+[EquationOfState]
+idealgas_gamma = 3.0
+eos = "idealgas"
+
+[EulerEq]
+id = "smooth_isentropic_flow"
+av_regularization = "navierstokes"
+
+[Evolution]
+cfl = 0.05
+tend = 0.50
+
+[Mesh]
+range = [-1.0, 1.0]
+n = 5
+k = 680
+basis = "lgl"
+
+[Output]
+variables = ["rho", "q", "E"]
+aligned_ts = "$(collect(range(0.01,0.50,step=0.01)))"
+
+[HRSC]
+method = "av"
+av_method = "entropy"
+entropy_ce = 20.0
diff --git a/examples/sod_shock_tube/fv.toml b/examples/sod_shock_tube/fv.toml
new file mode 100644
index 0000000000000000000000000000000000000000..e9880d8a41397700e244a22ad09e683c8a0956b0
--- /dev/null
+++ b/examples/sod_shock_tube/fv.toml
@@ -0,0 +1,23 @@
+[EquationOfState]
+idealgas_gamma = 1.4
+eos = "idealgas"
+
+[EulerEq]
+id = "sod_shock_tube"
+bc = "from_id"
+
+[Evolution]
+cfl = 1.0
+tend = 0.5
+
+[Mesh]
+periodic = [false]
+range = [-1.0, 1.0]
+k = 1536
+scheme = "FV"
+
+[Output]
+every_iteration = 0
+every_dt        = 0.0
+aligned_ts      = "$(collect(range(0.01,0.5,step=0.01)))"
+variables       = [ "rho", "q", "E" ]
diff --git a/examples/sod_shock_tube/fv3.toml b/examples/sod_shock_tube/fv3.toml
new file mode 100644
index 0000000000000000000000000000000000000000..5fa047809259994b6763466b0f9a7cf3848ed95f
--- /dev/null
+++ b/examples/sod_shock_tube/fv3.toml
@@ -0,0 +1,23 @@
+[EquationOfState]
+idealgas_gamma = 1.4
+eos = "idealgas"
+
+[EulerEq]
+id = "sod_shock_tube"
+bc = "from_id"
+
+[Evolution]
+cfl = 1.0
+tend = 0.5
+
+[Mesh]
+periodic = [false]
+range = [-1.0, 1.0]
+k = 768
+scheme = "FV"
+
+[Output]
+every_iteration = 0
+every_dt        = 0.0
+aligned_ts      = "$(collect(range(0.01,0.5,step=0.01)))"
+variables       = [ "rho", "q", "E" ]
diff --git a/examples/sod_shock_tube/fv4.toml b/examples/sod_shock_tube/fv4.toml
new file mode 100644
index 0000000000000000000000000000000000000000..bc12f87c94f53e35ca1a0bce139866727e3daffa
--- /dev/null
+++ b/examples/sod_shock_tube/fv4.toml
@@ -0,0 +1,23 @@
+[EquationOfState]
+idealgas_gamma = 1.4
+eos = "idealgas"
+
+[EulerEq]
+id = "sod_shock_tube"
+bc = "from_id"
+
+[Evolution]
+cfl = 1.0
+tend = 0.5
+
+[Mesh]
+periodic = [false]
+range = [-1.0, 1.0]
+k = 1024
+scheme = "FV"
+
+[Output]
+every_iteration = 0
+every_dt        = 0.0
+aligned_ts      = "$(collect(range(0.01,0.5,step=0.01)))"
+variables       = [ "rho", "q", "E" ]
diff --git a/examples/sod_shock_tube/fv5.toml b/examples/sod_shock_tube/fv5.toml
new file mode 100644
index 0000000000000000000000000000000000000000..f75e9fff714a167396934c1d78c460732aa4072e
--- /dev/null
+++ b/examples/sod_shock_tube/fv5.toml
@@ -0,0 +1,23 @@
+[EquationOfState]
+idealgas_gamma = 1.4
+eos = "idealgas"
+
+[EulerEq]
+id = "sod_shock_tube"
+bc = "from_id"
+
+[Evolution]
+cfl = 1.0
+tend = 0.5
+
+[Mesh]
+periodic = [false]
+range = [-1.0, 1.0]
+k = 1280
+scheme = "FV"
+
+[Output]
+every_iteration = 0
+every_dt        = 0.0
+aligned_ts      = "$(collect(range(0.01,0.5,step=0.01)))"
+variables       = [ "rho", "q", "E" ]
diff --git a/examples/sod_shock_tube/fv6.toml b/examples/sod_shock_tube/fv6.toml
new file mode 100644
index 0000000000000000000000000000000000000000..e9880d8a41397700e244a22ad09e683c8a0956b0
--- /dev/null
+++ b/examples/sod_shock_tube/fv6.toml
@@ -0,0 +1,23 @@
+[EquationOfState]
+idealgas_gamma = 1.4
+eos = "idealgas"
+
+[EulerEq]
+id = "sod_shock_tube"
+bc = "from_id"
+
+[Evolution]
+cfl = 1.0
+tend = 0.5
+
+[Mesh]
+periodic = [false]
+range = [-1.0, 1.0]
+k = 1536
+scheme = "FV"
+
+[Output]
+every_iteration = 0
+every_dt        = 0.0
+aligned_ts      = "$(collect(range(0.01,0.5,step=0.01)))"
+variables       = [ "rho", "q", "E" ]
diff --git a/examples/sod_shock_tube/fv7.toml b/examples/sod_shock_tube/fv7.toml
new file mode 100644
index 0000000000000000000000000000000000000000..0e968a8dc7ad89576e4225d8d13500bc56244bc0
--- /dev/null
+++ b/examples/sod_shock_tube/fv7.toml
@@ -0,0 +1,23 @@
+[EquationOfState]
+idealgas_gamma = 1.4
+eos = "idealgas"
+
+[EulerEq]
+id = "sod_shock_tube"
+bc = "from_id"
+
+[Evolution]
+cfl = 1.0
+tend = 0.5
+
+[Mesh]
+periodic = [false]
+range = [-1.0, 1.0]
+k = 1792
+scheme = "FV"
+
+[Output]
+every_iteration = 0
+every_dt        = 0.0
+aligned_ts      = "$(collect(range(0.01,0.5,step=0.01)))"
+variables       = [ "rho", "q", "E" ]
diff --git a/examples/sod_shock_tube/mono_n3.toml b/examples/sod_shock_tube/mono_n3.toml
new file mode 100644
index 0000000000000000000000000000000000000000..cac7d55ec4643646c6d3d448f270aa8b1cac3d94
--- /dev/null
+++ b/examples/sod_shock_tube/mono_n3.toml
@@ -0,0 +1,33 @@
+[EquationOfState]
+idealgas_gamma = 1.4
+eos = "idealgas"
+
+[EulerEq]
+id = "sod_shock_tube"
+bc = "from_id"
+av_regularization = "mono"
+
+[Mesh]
+periodic = [ false ]
+range = [ -1.0, 1.0 ]
+cfl = 0.75
+n = 3
+k = 256
+basis = "lgl"
+
+[Output]
+every_iteration = 0
+every_dt        = 0.0
+aligned_ts      = "$(collect(range(0.01,0.5,step=0.01)))"
+variables       = [ "rho", "q", "E", "ldg_rho", "ldg_q", "ldg_E", "smoothed_mu" ]
+
+[Evolution]
+tend = 0.5
+
+[Log]
+progress_stdout = true
+
+[HRSC]
+method = "av"
+av_method = "entropy"
+entropy_ce = 5.0
diff --git a/examples/sod_shock_tube/mono_n4.toml b/examples/sod_shock_tube/mono_n4.toml
new file mode 100644
index 0000000000000000000000000000000000000000..756994abfe13df0bd5e030225c3f778311e4c3a3
--- /dev/null
+++ b/examples/sod_shock_tube/mono_n4.toml
@@ -0,0 +1,33 @@
+[EquationOfState]
+idealgas_gamma = 1.4
+eos = "idealgas"
+
+[EulerEq]
+id = "sod_shock_tube"
+bc = "from_id"
+av_regularization = "mono"
+
+[Mesh]
+periodic = [ false ]
+range = [ -1.0, 1.0 ]
+cfl = 0.75
+n = 4
+k = 256
+basis = "lgl"
+
+[Output]
+every_iteration = 0
+every_dt        = 0.0
+aligned_ts      = "$(collect(range(0.01,0.5,step=0.01)))"
+variables       = [ "rho", "q", "E", "ldg_rho", "ldg_q", "ldg_E", "smoothed_mu" ]
+
+[Evolution]
+tend = 0.5
+
+[Log]
+progress_stdout = true
+
+[HRSC]
+method = "av"
+av_method = "entropy"
+entropy_ce = 5.0
diff --git a/examples/sod_shock_tube/mono_n5.toml b/examples/sod_shock_tube/mono_n5.toml
new file mode 100644
index 0000000000000000000000000000000000000000..0a77c7af38e9a1c23d1de2cd5a5a3077de1515a4
--- /dev/null
+++ b/examples/sod_shock_tube/mono_n5.toml
@@ -0,0 +1,33 @@
+[EquationOfState]
+idealgas_gamma = 1.4
+eos = "idealgas"
+
+[EulerEq]
+id = "sod_shock_tube"
+bc = "from_id"
+av_regularization = "mono"
+
+[Mesh]
+periodic = [ false ]
+range = [ -1.0, 1.0 ]
+cfl = 0.75
+n = 5
+k = 256
+basis = "lgl"
+
+[Output]
+every_iteration = 0
+every_dt        = 0.0
+aligned_ts      = "$(collect(range(0.01,0.5,step=0.01)))"
+variables       = [ "rho", "q", "E", "ldg_rho", "ldg_q", "ldg_E", "smoothed_mu" ]
+
+[Evolution]
+tend = 0.5
+
+[Log]
+progress_stdout = true
+
+[HRSC]
+method = "av"
+av_method = "entropy"
+entropy_ce = 5.0
diff --git a/examples/sod_shock_tube/mono_n6.toml b/examples/sod_shock_tube/mono_n6.toml
new file mode 100644
index 0000000000000000000000000000000000000000..e006699f04ef244ccb9b9ac1ce329015a6ac38f2
--- /dev/null
+++ b/examples/sod_shock_tube/mono_n6.toml
@@ -0,0 +1,33 @@
+[EquationOfState]
+idealgas_gamma = 1.4
+eos = "idealgas"
+
+[EulerEq]
+id = "sod_shock_tube"
+bc = "from_id"
+av_regularization = "mono"
+
+[Mesh]
+periodic = [ false ]
+range = [ -1.0, 1.0 ]
+cfl = 0.75
+n = 6
+k = 256
+basis = "lgl"
+
+[Output]
+every_iteration = 0
+every_dt        = 0.0
+aligned_ts      = "$(collect(range(0.01,0.5,step=0.01)))"
+variables       = [ "rho", "q", "E", "ldg_rho", "ldg_q", "ldg_E", "smoothed_mu" ]
+
+[Evolution]
+tend = 0.5
+
+[Log]
+progress_stdout = true
+
+[HRSC]
+method = "av"
+av_method = "entropy"
+entropy_ce = 5.0
diff --git a/examples/sod_shock_tube/mono_n7.toml b/examples/sod_shock_tube/mono_n7.toml
new file mode 100644
index 0000000000000000000000000000000000000000..a3670c3e3696c26c845b3071f18274a30503f593
--- /dev/null
+++ b/examples/sod_shock_tube/mono_n7.toml
@@ -0,0 +1,33 @@
+[EquationOfState]
+idealgas_gamma = 1.4
+eos = "idealgas"
+
+[EulerEq]
+id = "sod_shock_tube"
+bc = "from_id"
+av_regularization = "mono"
+
+[Mesh]
+periodic = [ false ]
+range = [ -1.0, 1.0 ]
+cfl = 0.75
+n = 7
+k = 256
+basis = "lgl"
+
+[Output]
+every_iteration = 0
+every_dt        = 0.0
+aligned_ts      = "$(collect(range(0.01,0.5,step=0.01)))"
+variables       = [ "rho", "q", "E", "ldg_rho", "ldg_q", "ldg_E", "smoothed_mu" ]
+
+[Evolution]
+tend = 0.5
+
+[Log]
+progress_stdout = true
+
+[HRSC]
+method = "av"
+av_method = "entropy"
+entropy_ce = 5.0
diff --git a/src/EulerEq/EulerEq.jl b/src/EulerEq/EulerEq.jl
index 63a03ca6512def33222b9668f9b283b39ebd4e58..0459bac972fb2b1f502999ed6631aab9e389a235 100644
--- a/src/EulerEq/EulerEq.jl
+++ b/src/EulerEq/EulerEq.jl
@@ -49,6 +49,22 @@ dg1d.@parameters EulerEq begin
   bc = "periodic"
   @check bc in [ "periodic", "from_id" ]
 
+  """
+  Type of artificial viscosity regularization. Available options:
+  - `none`
+  - `mono` - monolithic regularization, see [1]
+  - `navierstokes` - regularization based on Navier-Stokes equations, but without heat flux, see [1]
+  - `general` - most general regularization, see [1]; actually implemented version of [2]
+
+  [1] J.-L. Guermond and B. Popov, “Viscous regularization of the euler equations and entropy principles,”
+  SIAM Journal on Applied Mathematics, vol. 74, no. 2, pp. 284–305, 2014.
+  [2] J.-L. Guermond, B. Popov, and V. Tomov, “Entropy–viscosity method for the single material euler equa-
+  tions in lagrangian frame,” Computer Methods in Applied Mechanics and Engineering, vol. 300, pp. 402–
+  426, 2016.
+  """
+  av_regularization = "none"
+  @check av_regularization in ["none", "mono", "navierstokes", "general"]
+
 end
 
 end # module
diff --git a/src/EulerEq/callbacks.jl b/src/EulerEq/callbacks.jl
index 1d2f78fd0225a457545a700526aa88d596a27572..fcda2346035a4904abed020c7b8d580a29c81ad3 100644
--- a/src/EulerEq/callbacks.jl
+++ b/src/EulerEq/callbacks.jl
@@ -45,9 +45,6 @@ function callback_equation(state_u, state_t, env, P, isperiodic, eq::EulerEquati
   t       .= state_t
   Sm1     .= S
   flx_Sm1 .= flx_S
-  dg1d.new_broadcast_volume!(pressure, equation, env.mesh)
-  @unpack rho, q = get_dynamic_variables(cache)
-  @unpack p, eps = get_static_variables(cache)
   dg1d.new_broadcast_volume!(entropy_variables, equation, env.mesh)
 end
 
@@ -74,7 +71,7 @@ function callback_hrsc(state_u, state_t, env, P, isperiodic, hrsc::HRSC.Abstract
   @unpack mu, max_v   = get_cell_variables(cache)
   @unpack v           = get_static_variables(cache)
   @unpack rho         = get_dynamic_variables(cache)
-  dg1d.new_broadcast_volume!(maxspeed, equation, cache)
+  dg1d.new_broadcast_volume!(maxspeed, equation, mesh)
   Npts, K = layout(mesh)
   mat_max_v = reshape(view(v, :), (Npts, K))
   for k = 1:K
diff --git a/src/EulerEq/equation.jl b/src/EulerEq/equation.jl
index df6f37fb68abf8c6dc25053cf52214bd28027e1c..0a843de64149f90542cdd209f4a2964b2c0bff81 100644
--- a/src/EulerEq/equation.jl
+++ b/src/EulerEq/equation.jl
@@ -448,3 +448,71 @@ end
 
   @returns nflx_rho, nflx_qx, nflx_qy, nflx_E
 end
+
+
+# TODO Rename all methods accordingly
+const mono_ldg_nflux   = ldg_nflux
+const mono_ldg_bdryllf = ldg_bdryllf
+const mono_av_flux     = av_flux
+const mono_av_nflux    = av_nflux
+const mono_bdryllf     = bdryllf
+
+
+#######################################################################
+#                    Navier-Stokes regularization                     #
+#######################################################################
+
+
+@with_signature [legacy=false] function navierstokes_u(equation::EulerEquation)
+  @accepts rho, q
+  u = q/rho
+  @returns u
+end
+
+@with_signature [legacy=false] function navierstokes_ldg_nflux(eq::EulerEquation)
+  @accepts        u
+  @accepts [bdry] bdry_u
+  @accepts [bdry] nx
+
+  nflx_u = 0.5 * nx * (u + bdry_u)
+
+  @returns [bdry] nflx_u
+end
+@with_signature [legacy=false] function navierstokes_ldg_bdryllf(eq::EulerEquation)
+  @accepts        u
+  @accepts [bdry] nx
+
+  # set exterior state to 0 according to stability analysis
+  # see dg-notes/tex/evm.pdf for reasoning
+  nflx_u = nx * u
+
+  @returns [bdry] nflx_u
+end
+
+@with_signature [legacy=false] function navierstokes_av_flux(eq::EulerEquation)
+  @accepts flx_q, flx_E, ldg_u, u, smoothed_mu
+
+  # μ = λ, κ = 0
+  G = 4 * smoothed_mu * ldg_u + smoothed_mu * ldg_u
+  flx_q   += G
+  flx_E   += G * u
+
+  @returns flx_q, flx_E
+end
+@with_signature [legacy=false] function navierstokes_av_nflux(eq::EulerEquation)
+  @accepts        ldg_u, u, smoothed_mu
+  @accepts [bdry] bdry_ldg_u, bdry_u, bdry_smoothed_mu
+  @accepts [bdry] nflx_q, nflx_E, nx
+
+  G      = 4 * smoothed_mu * ldg_u + smoothed_mu * ldg_u
+  bdry_G = 4 * bdry_smoothed_mu * bdry_ldg_u + bdry_smoothed_mu * bdry_ldg_u
+  nflx       = nx*G
+  bdry_nflx  = nx*bdry_G
+  nflx_q    += 0.5 * (nflx + bdry_nflx)
+
+  nflx       = nx*G*u
+  bdry_nflx  = nx*bdry_G*bdry_u
+  nflx_E    += 0.5 * (nflx + bdry_nflx)
+
+  @returns [bdry] nflx_q, nflx_E
+end
diff --git a/src/EulerEq/monoreg.jl.bak b/src/EulerEq/monoreg.jl.bak
new file mode 100644
index 0000000000000000000000000000000000000000..70105050b5db04f06a7d52cd186b0ff18215907d
--- /dev/null
+++ b/src/EulerEq/monoreg.jl.bak
@@ -0,0 +1,29 @@
+#######################################################################
+#                            Flux methods                             #
+#######################################################################
+
+
+
+@with_signature [legacy=false] function mono_ldg_nflux(eq::EulerEquation)
+  @accepts        rho, q, E
+  @accepts [bdry] bdry_rho, bdry_q, bdry_E
+  @accepts [bdry] nx
+
+  nflx_rho = 0.5 * nx * (rho + bdry_rho)
+  nflx_q   = 0.5 * nx * (q + bdry_q)
+  nflx_E   = 0.5 * nx * (E + bdry_E)
+
+  @returns [bdry] nflx_rho, nflx_q, nflx_E
+end
+@with_signature [legacy=false] function mono_ldg_bdryllf(eq::EulerEquation)
+  @accepts        rho, q, E
+  @accepts [bdry] nx
+
+  # set exterior state to 0 according to stability analysis
+  # see dg-notes/tex/evm.pdf for reasoning
+  nflx_rho = nx * rho
+  nflx_q   = nx * q
+  nflx_E   = nx * E
+
+  @returns [bdry] nflx_rho, nflx_q, nflx_E
+end
diff --git a/src/EulerEq/rhs.jl b/src/EulerEq/rhs.jl
index 862b6788c07767099ed462928c6a14072f30679e..c0e3e39ba801e1b3ce9c3d212b2eb49a86e247ff 100644
--- a/src/EulerEq/rhs.jl
+++ b/src/EulerEq/rhs.jl
@@ -213,20 +213,33 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mesh1d,
     bdryconds, ldg_bdryconds, av_bdryconds)
 
   @unpack cache, mesh = env
-  @unpack equation    = P
+  @unpack equation, prms = P
+  reg = prms.av_regularization
+
+  if reg === :mono
+    rhs_mono!(cache, mesh, equation, P)
+  elseif reg === :navierstokes
+    rhs_navierstokes!(cache, mesh, equation, P)
+  elseif reg === :general
+    rhs_general!(cache, mesh, equation, P)
+  else
+    TODO(reg)
+  end
+end
 
-  @unpack rho, q, E                            = get_dynamic_variables(cache)
+
+function rhs_mono!(cache, mesh, equation, P)
+
+  @unpack rho, q, E                = get_dynamic_variables(cache)
   @unpack flx_rho, flx_q, flx_E, p, eps,
-          ldg_rho, ldg_q, ldg_E,
-          flx_ldg_rho, flx_ldg_q, flx_ldg_E    = get_static_variables(cache)
-  @unpack rhs_rho, rhs_q, rhs_E                = get_rhs_variables(cache)
-  @unpack nflx_rho, nflx_q, nflx_E,
-          nflx_ldg_rho, nflx_ldg_q, nflx_ldg_E = get_bdry_variables(cache)
+          ldg_rho, ldg_q, ldg_E    = get_static_variables(cache)
+  @unpack rhs_rho, rhs_q, rhs_E    = get_rhs_variables(cache)
+  @unpack nflx_rho, nflx_q, nflx_E = get_bdry_variables(cache)
   @unpack bdry_rho, bdry_q, bdry_E,
           bdry_v,
           bdry_ldg_rho, bdry_ldg_q, bdry_ldg_E,
-          bdry_smoothed_mu                     = get_bdry_variables(mesh.cache)
-  @unpack v, smoothed_mu                       = get_static_variables(mesh.cache)
+          bdry_smoothed_mu         = get_bdry_variables(mesh.cache)
+  @unpack v, smoothed_mu           = get_static_variables(mesh.cache)
 
   dg1d.interpolate_face_data!(mesh, rho, bdry_rho)
   dg1d.interpolate_face_data!(mesh, q, bdry_q)
@@ -234,8 +247,8 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mesh1d,
   dg1d.interpolate_face_data!(mesh, v, bdry_v)
 
   ## solve auxiliary equation: q + ∂x u = 0
-  dg1d.new_broadcast_faces!(ldg_nflux, equation, mesh)
-  dg1d.new_broadcast_bdry!(ldg_bdryllf, equation, mesh)
+  dg1d.new_broadcast_faces!(mono_ldg_nflux, equation, mesh)
+  dg1d.new_broadcast_bdry!(mono_ldg_bdryllf, equation, mesh)
   compute_rhs_weak_form!(ldg_rho, rho, nflx_rho, mesh)
   compute_rhs_weak_form!(ldg_q,   q,   nflx_q,   mesh)
   compute_rhs_weak_form!(ldg_E,   E,   nflx_E,   mesh)
@@ -249,8 +262,51 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mesh1d,
   dg1d.interpolate_face_data!(mesh, ldg_E, bdry_ldg_E)
   dg1d.interpolate_face_data!(mesh, smoothed_mu, bdry_smoothed_mu)
 
-  dg1d.new_broadcast_volume!(av_flux, equation, mesh)
-  dg1d.new_broadcast_faces!(av_nflux, equation, mesh)
+  dg1d.new_broadcast_volume!(mono_av_flux, equation, mesh)
+  dg1d.new_broadcast_faces!(mono_av_nflux, equation, mesh)
+  dg1d.new_broadcast_bdry!(mono_bdryllf, equation, mesh)
+  # TODO don't we also need av_bdrylff here?
+
+  compute_rhs_weak_form!(rhs_rho, flx_rho,  nflx_rho, mesh)
+  compute_rhs_weak_form!(rhs_q,   flx_q,    nflx_q,   mesh)
+  compute_rhs_weak_form!(rhs_E,   flx_E,    nflx_E,   mesh)
+
+  return
+end
+function rhs_navierstokes!(cache, mesh, equation, P)
+
+  @unpack rho, q, E                        = get_dynamic_variables(cache)
+  @unpack flx_rho, flx_q, flx_E, u,
+          ldg_u                            = get_static_variables(cache)
+  @unpack rhs_rho, rhs_q, rhs_E            = get_rhs_variables(cache)
+  @unpack nflx_rho, nflx_q, nflx_E, nflx_u = get_bdry_variables(cache)
+  @unpack bdry_rho, bdry_q, bdry_E,
+          bdry_v, bdry_u,
+          bdry_ldg_u,
+          bdry_smoothed_mu                 = get_bdry_variables(mesh.cache)
+  @unpack v, smoothed_mu                   = get_static_variables(mesh.cache)
+
+  dg1d.new_broadcast_volume!(navierstokes_u, equation, mesh)
+  dg1d.interpolate_face_data!(mesh, rho, bdry_rho)
+  dg1d.interpolate_face_data!(mesh, q,   bdry_q)
+  dg1d.interpolate_face_data!(mesh, E,   bdry_E)
+  dg1d.interpolate_face_data!(mesh, v,   bdry_v)
+  dg1d.interpolate_face_data!(mesh, u,   bdry_u)
+
+  ## solve auxiliary equation: q + ∂x u = 0
+  dg1d.new_broadcast_faces!(navierstokes_ldg_nflux, equation, mesh)
+  dg1d.new_broadcast_bdry!(navierstokes_ldg_bdryllf, equation, mesh)
+  compute_rhs_weak_form!(ldg_u, u, nflx_u, mesh)
+
+  ## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
+  dg1d.new_broadcast_volume!(flux, equation, mesh)
+  dg1d.new_broadcast_faces!(llf, equation, mesh)
+
+  dg1d.interpolate_face_data!(mesh, ldg_u,       bdry_ldg_u)
+  dg1d.interpolate_face_data!(mesh, smoothed_mu, bdry_smoothed_mu)
+
+  dg1d.new_broadcast_volume!(navierstokes_av_flux, equation, mesh)
+  dg1d.new_broadcast_faces!(navierstokes_av_nflux, equation, mesh)
   dg1d.new_broadcast_bdry!(bdryllf, equation, mesh)
   # TODO don't we also need av_bdrylff here?
 
@@ -260,6 +316,56 @@ function rhs!(env, P::Project, hrsc::HRSC.AbstractArtificialViscosity, ::Mesh1d,
 
   return
 end
+function rhs_general!(cache, mesh, equation)
+
+  TODO()
+
+  @unpack rho, q, E                            = get_dynamic_variables(cache)
+  @unpack flx_rho, flx_q, flx_E, p, eps,
+          ldg_rho, ldg_rhoe, ldg_u,
+          flx_ldg_rho, flx_ldg_q, flx_ldg_E    = get_static_variables(cache)
+  @unpack rhs_rho, rhs_q, rhs_E                = get_rhs_variables(cache)
+  @unpack nflx_rho, nflx_q, nflx_E,
+          nflx_ldg_rho, nflx_ldg_q, nflx_ldg_E = get_bdry_variables(cache)
+  @unpack bdry_rho, bdry_q, bdry_E,
+          bdry_v,
+          bdry_ldg_rho, bdry_ldg_rhoe, bdry_ldg_u,
+          bdry_smoothed_mu                     = get_bdry_variables(mesh.cache)
+  @unpack v, smoothed_mu                       = get_static_variables(mesh.cache)
+
+  dg1d.new_broadcast_volume!(mono_rhoeps, equation, mesh)
+  dg1d.interpolate_face_data!(mesh, rho,    bdry_rho)
+  dg1d.interpolate_face_data!(mesh, rhoeps, bdry_rhoeps)
+  dg1d.interpolate_face_data!(mesh, u,      bdry_u)
+  dg1d.interpolate_face_data!(mesh, v,      bdry_v)
+
+  ## solve auxiliary equation: q + ∂x u = 0
+  dg1d.new_broadcast_faces!(general_ldg_nflux, equation, mesh)
+  dg1d.new_broadcast_bdry!(general_ldg_bdryllf, equation, mesh)
+  compute_rhs_weak_form!(ldg_rho, rho,    nflx_rho,    mesh)
+  compute_rhs_weak_form!(ldg_q,   rhoeps, nflx_rhoeps, mesh)
+  compute_rhs_weak_form!(ldg_E,   u,      nflx_u,      mesh)
+
+  ## compute rhs of regularized equation: ∂t u + ∂x f + ∂x μ q = 0
+  dg1d.new_broadcast_volume!(flux, equation, mesh)
+  dg1d.new_broadcast_faces!(llf, equation, mesh)
+
+  dg1d.interpolate_face_data!(mesh, ldg_rho,    bdry_ldg_rho)
+  dg1d.interpolate_face_data!(mesh, ldg_rhoeps, bdry_ldg_rhoeps)
+  dg1d.interpolate_face_data!(mesh, ldg_u,      bdry_ldg_u)
+  dg1d.interpolate_face_data!(mesh, smoothed_mu, bdry_smoothed_mu)
+
+  dg1d.new_broadcast_volume!(general_av_flux, equation, mesh)
+  dg1d.new_broadcast_faces!(general_av_nflux, equation, mesh)
+  dg1d.new_broadcast_bdry!(general_bdryllf, equation, mesh)
+  # TODO don't we also need av_bdrylff here?
+
+  compute_rhs_weak_form!(rhs_rho, flx_rho,  nflx_rho, mesh)
+  compute_rhs_weak_form!(rhs_q,   flx_q,    nflx_q,   mesh)
+  compute_rhs_weak_form!(rhs_E,   flx_E,    nflx_E,   mesh)
+
+  return
+end
 
 
 function rhs!(env, P::Project2d, hrsc::HRSC.AbstractArtificialViscosity, mesh::Mesh2d,
diff --git a/src/EulerEq/setup.jl b/src/EulerEq/setup.jl
index 3879b91e68f1f457c5da758159585df91af735ea..6a44eceb710415f59859393023d16817c0f3194a 100644
--- a/src/EulerEq/setup.jl
+++ b/src/EulerEq/setup.jl
@@ -10,7 +10,9 @@ function Project(env::Environment, mesh::Mesh1d, prms)
   hrsc = HRSC.make_HRSC(env.mesh, prms["HRSC"])
 
   bdrycond = dg1d.DirichletBC2()
-  P = Project(equation, hrsc, tci, bdrycond)
+  fixedprms = (; av_regularization=Symbol(prms["EulerEq"]["av_regularization"]),
+                 bernstein=BernsteinReconstruction(mesh))
+  P = Project(equation, hrsc, tci, bdrycond, fixedprms)
 
   # register variables
   # TODO add a ::Nothing overload for register_variables!
@@ -18,6 +20,8 @@ function Project(env::Environment, mesh::Mesh1d, prms)
   @unpack cache = env
   _register_variables!(mesh, P)
   _register_variables!(mesh, bdrycond)
+  _register_variables!(mesh, P.hrsc, prms["EulerEq"]["av_regularization"])
+  _register_variables!(mesh, P.tci)
   display(cache)
 
   # setup initial data
@@ -57,7 +61,7 @@ function Project(env::Environment, mesh::Mesh2d, prms)
   @unpack cache = env
   _register_variables!(mesh, P, mesh)
   _register_variables!(mesh, bdrycond)
-  _register_variables!(mesh, hrsc)
+  _register_variables!(mesh, hrsc, :mono)
   display(cache)
 
   # setup initial data
@@ -99,37 +103,50 @@ function _register_variables!(mesh, P::Project)
                              :bdry_rho, :bdry_q, :bdry_E, :bdry_v),
     global_variablenames  = (:t, :tm1) # recent time stamps
   )
-  _register_variables!(mesh, P.hrsc)
-  _register_variables!(mesh, P.tci)
 end
 
 
-_register_variables!(mesh, tci_or_hrsc::Nothing) = nothing
+_register_variables!(mesh, tci_or_hrsc::Nothing, args...) = nothing
 
-_register_variables!(mesh, hrsc::AbstractReconstruction) = nothing
+_register_variables!(mesh, hrsc::AbstractReconstruction, regularization) = nothing
 
 
-function _register_variables!(mesh, hrsc::HRSC.AbstractArtificialViscosity)
+function _register_variables!(mesh, hrsc::HRSC.AbstractArtificialViscosity, regularization)
   register_variables!(mesh.cache,
-    static_variablenames  = (:ldg_rho, :ldg_q, :ldg_E,              # local DG variable
-                             :flx_ldg_rho, :flx_ldg_q, :flx_ldg_E), # local DG flux
-    bdry_variablenames    = (:nflx_ldg_rho, :nflx_ldg_q, :nflx_ldg_E), # local DG flux
-    cell_variablenames    = (:mu,) # 'one viscosity to rule them all'
+                      cell_variablenames    = (:mu,) # 'one viscosity to rule them all'
+  )
+  if regularization == "mono"
+    register_variables!(mesh.cache,
+                        static_variablenames  = (:ldg_rho, :ldg_q, :ldg_E,),
+                        bdry_variablenames    = (:nflx_ldg_rho, :nflx_ldg_q, :nflx_ldg_E,
+                                                 :bdry_ldg_rho, :bdry_ldg_q, :bdry_ldg_E),
+    )
+  elseif regularization == "navierstokes"
+    register_variables!(mesh.cache,
+                        static_variablenames  = (:ldg_u, :u),
+                        bdry_variablenames = (:bdry_ldg_u, :bdry_u, :nflx_u),
+    )
+  elseif regularization == "general"
+    TODO()
+    register_variables!(mesh.cache,
+      static_variablenames  = (:ldg_rho, :ldg_q, :ldg_E,              # local DG variable
+                               :flx_ldg_rho, :flx_ldg_q, :flx_ldg_E), # local DG flux
+      bdry_variablenames    = (:nflx_ldg_rho, :nflx_ldg_q, :nflx_ldg_E), # local DG flux
     )
+  end
 end
 
 
-function _register_variables!(mesh::Mesh, hrsc::HRSC.SmoothedArtificialViscosity)
-  _register_variables!(mesh, hrsc.av)
+function _register_variables!(mesh::Mesh, hrsc::HRSC.SmoothedArtificialViscosity, regularization)
+  _register_variables!(mesh, hrsc.av, regularization)
   register_variables!(mesh.cache,
     static_variablenames  = (:smoothed_mu,),
-    bdry_variablenames    = (:bdry_smoothed_mu,
-                            :bdry_ldg_rho, :bdry_ldg_q, :bdry_ldg_E),
+    bdry_variablenames    = (:bdry_smoothed_mu,)
   )
 end
 
 
-function _register_variables!(mesh::Mesh2d, hrsc::HRSC.EntropyViscosity)
+function _register_variables!(mesh::Mesh2d, hrsc::HRSC.EntropyViscosity, regularization)
   register_variables!(mesh.cache,
     static_variablenames = (:S, :Sm1, :flx_S_x, :flx_S_y, :flxm1_S_x, :flxm1_S_y,
                             :ldg_rho_x, :ldg_rho_y, :ldg_qx_x, :ldg_qx_y,
diff --git a/src/EulerEq/types.jl b/src/EulerEq/types.jl
index 2af4af9ac6394bd78e89f798164ad5f0efc5d895..eeb7f938fba186c8bed2c4501d4e513b03399a32 100644
--- a/src/EulerEq/types.jl
+++ b/src/EulerEq/types.jl
@@ -10,12 +10,14 @@ end
 mutable struct Project{T_Equation <:EulerEquation,
                        T_HRSC     <:Maybe{HRSC.AbstractHRSC},
                        T_TCI      <:Maybe{TCI.AbstractTCI},
-                       T_BC}
+                       T_BC,
+                       T_Prms<:NamedTuple}
 
   equation::T_Equation
   hrsc::T_HRSC
   tci::T_TCI
   bdrycond::T_BC
+  prms::T_Prms
 end
 
 mutable struct Project2d{T_Equation <:EulerEquation2d,
diff --git a/src/HRSC/HRSC.jl b/src/HRSC/HRSC.jl
index 20e5995e7f3369b746972f67adb7dc9aa75959f8..d8c7c99fff0f3a725a004e36d37fff66e7e0c3b7 100644
--- a/src/HRSC/HRSC.jl
+++ b/src/HRSC/HRSC.jl
@@ -45,7 +45,7 @@ export AbstractHRSC,
   TODO
   """
   av_smoother_order = 1
-  @check av_smoother_order in [ 1, 2 ]
+  @check av_smoother_order in [ 0, 1, 2 ]
 
   "TODO"
   mda_cmax = 1.0
@@ -138,7 +138,7 @@ function make_HRSC(mesh::Mesh1d, prms)
 
     if av_smoother_order < 0
       error("HRSC: 'av_smoother_order' must be positive, found $av_smoother_order")
-    elseif av_smoother_order > 0
+    elseif av_smoother_order >= 0
       return SmoothedArtificialViscosity(mesh, hrsc, smoother_order=av_smoother_order)
     else
       return hrsc
diff --git a/src/HRSC/bernstein.jl b/src/HRSC/bernstein.jl
index ca29a4400cd703fb09b3c7441bb6df4927cd3536..6d42466202c4936013e9b12dc08e128d85f1d243 100644
--- a/src/HRSC/bernstein.jl
+++ b/src/HRSC/bernstein.jl
@@ -60,7 +60,7 @@ Compute inplace a convex combination of `u!` and its Bernstein reconstruction we
 grid to the bernstein grid and back, respectively.
 `m` and `M` keywords allow to enforce lower and upper bound in reconstruction.
 """
-function reconstruct_cell!(u, w, P, V; m=missing, M=missing)
+@inline function reconstruct_cell!(u, w, P, V; m=missing, M=missing)
 
   uB = P * u
 
diff --git a/src/HRSC/ev.jl b/src/HRSC/ev.jl
index e0142e95d9cdd52842eda10b5ff754d2b4f81edf..e3c4e56e32482de0956277b26da9d89ffb865fa1 100644
--- a/src/HRSC/ev.jl
+++ b/src/HRSC/ev.jl
@@ -62,7 +62,7 @@ function compute_viscosity!(
 
   # upper bound on AV
   maxabs_v = vec(maximum(abs.(v), dims=1))
-  mu_max = @. cmax * maxabs_v * dl / N
+  (mu_max,) = @. cmax * maxabs_v * dl / N
 
   # final AV
   @. mu = min(mu_e, mu_max)
diff --git a/src/HRSC/smoother.jl b/src/HRSC/smoother.jl
index 7fcf16408015d5e231285c328583d37b28b2a03a..b316a5a29cc58130b2fb58af4d7418cc740db1f2 100644
--- a/src/HRSC/smoother.jl
+++ b/src/HRSC/smoother.jl
@@ -15,15 +15,16 @@ end
 
 
 function SmoothedArtificialViscosity(mesh::Mesh1d, av::AbstractArtificialViscosity; smoother_order=2)
-  @toggled_assert smoother_order in ( 1, 2 )
+  @toggled_assert smoother_order in ( 0, 1, 2 )
   quadr = mesh.element.quadr
-  interp_element = dg1d.SpectralElement(smoother_order, quadr)
-  if smoother_order == 1
+  if smoother_order == 0
+    return SmoothedArtificialViscosity(av, smooth_C0_0th_order!)
+  elseif smoother_order == 1
     return SmoothedArtificialViscosity(av, smooth_C0_1st_order!)
   elseif smoother_order == 2
     return SmoothedArtificialViscosity(av, smooth_C0_2nd_order!)
   else
-    error("smoother_order must be 1 or 2")
+    error("smoother_order must be 0, 1 or 2")
   end
 end
 
@@ -31,7 +32,6 @@ end
 function SmoothedArtificialViscosity(mesh::Mesh2d, av::AbstractArtificialViscosity; smoother_order=2)
   @toggled_assert smoother_order in ( 1, 2 )
   quadr = mesh.element.quadr
-  interp_element = dg1d.SpectralElement(smoother_order, quadr)
   if smoother_order == 1
     return SmoothedArtificialViscosity(av, smooth_C0_1st_order!)
   elseif smoother_order == 2
@@ -87,6 +87,24 @@ avg(a,b) = (a+b)/2
 avg(a,b,c,d) = (a+b+c+d)/4
 
 
+"""
+Dummy method which does not smoothing.
+"""
+function smooth_C0_0th_order!(bulk_v, cell_v, mesh::Mesh1d, isperiodic)
+  @unpack Npts, z = mesh.element
+  @unpack cells = mesh.tree
+  Nx = Npts
+  Kx, = mesh.tree.dims
+
+  mat_v = dg1d.vreshape(bulk_v, (Nx,Kx))
+  for k = 1:Kx, i = 1:Nx
+    mat_v[i,k] = cell_v[k]
+  end
+
+  return
+end
+
+
 """
 Smooth cell wise constant values over neighbouring grid cells
 and use a linear polynomial to interpolate nodal values.
diff --git a/src/TCI/entropyproduction.jl b/src/TCI/entropyproduction.jl
index b8717c115da7e2db21f776b9ac8acc2af3c9ce6a..1e3b7901b640c52f87c6835fab038b3857a104e2 100644
--- a/src/TCI/entropyproduction.jl
+++ b/src/TCI/entropyproduction.jl
@@ -36,7 +36,13 @@ function compute_indicator!(
   # entropy residual
   dt = t - tm1
   dt <= 0 && return # occurs on first time step when no previous entropy is available
-  EP .= abs.((E - Em1) / dt + 0.5  * D * (F + Fm1) .* invjac)
+
+  # use inplace operations
+  @. Fm1 -= F
+  mul!(EP, D, Fm1, -1.0, false)
+  @. EP .*= invjac
+  @. EP += (E - Em1) / dt
+  @. EP = abs(EP)
 
   # entropy jumps
   for k = 1:K
@@ -53,11 +59,14 @@ function compute_indicator!(
   # avg_E = dg1d.grid_average(E, mesh) # use as scale for indicator
   # @assert !(isapprox(avg_E, 0.0)) "No scale for entropy production available!"
 
+  # EP_max = @views maximum(EP[:])
+  _, ep_max = extrema(EP[:])
+
   @unpack EP_max = tci
   for k = 1:K
-    cell_EP = maximum( (maximum(EP[:,k]), EPjump_lhs[k], EPjump_rhs[k]) )
-    cell_EP = minimum((cell_EP, EP_max))
-    flag[k] = cell_EP / EP_max
+    # cell_EP = maximum( (maximum(EP[:,k]), EPjump_lhs[k], EPjump_rhs[k]) )
+    # cell_EP = minimum((cell_EP, EP_max))
+    flag[k] = @views maximum(EP[:,k] / ep_max)
   end
 
   return
diff --git a/src/glgl.jl b/src/glgl.jl
index af261d37f826c3084c579402ac311af16a552d73..11198584b2db11996305e5a93a3431157d656929 100644
--- a/src/glgl.jl
+++ b/src/glgl.jl
@@ -53,7 +53,6 @@ f1, f2 ... Vector of function values evaluated on grid x.
 function integrate(w, v, D, f)
 
   N = length(w)
-  @toggled_assert N == length(v)
   @toggled_assert (N,N) == size(D)
   @toggled_assert N == length(f)
 
diff --git a/src/utils.jl b/src/utils.jl
index 29648728fe86352cfe6a3ad9742d80171ba85132..b301b40a56676e1fd12b9fff38507ac52cc0bb46 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -270,7 +270,7 @@ end
 
 
 TODO() = error("Not implemented yet!")
-TODO(msg::String) = error("Not implemented: $msg")
+TODO(msg::Union{AbstractString,Symbol}) = error("Not implemented: $msg")
 TODO(fn::Function) = error("'$fn': Not implemented yet!")
 TODO(fn::Function, msg) = error("""
 '$fn'' Not implemented yet!
diff --git a/test/IntegrationTests/refs/burgers_sine_aventropy/output.h5 b/test/IntegrationTests/refs/burgers_sine_aventropy/output.h5
index 87fe09627e56d0b6195fdaf8444d15059ad2f16d..97849dfebb6b6f33e69570beef627b6a87ed6de2 100644
Binary files a/test/IntegrationTests/refs/burgers_sine_aventropy/output.h5 and b/test/IntegrationTests/refs/burgers_sine_aventropy/output.h5 differ
diff --git a/test/IntegrationTests/refs/euler_isentropic_flow_regulate-ns/euler_isentropic_flow_regulate-ns.toml b/test/IntegrationTests/refs/euler_isentropic_flow_regulate-ns/euler_isentropic_flow_regulate-ns.toml
new file mode 100644
index 0000000000000000000000000000000000000000..9ecd8390cdff27e23f85ff9d4ab12884b58d2f34
--- /dev/null
+++ b/test/IntegrationTests/refs/euler_isentropic_flow_regulate-ns/euler_isentropic_flow_regulate-ns.toml
@@ -0,0 +1,35 @@
+[EulerEq]
+id = "smooth_isentropic_flow"
+av_regularization = "navierstokes"
+# av_regularization = "mono"
+
+[EquationOfState]
+eos = "idealgas"
+idealgas_gamma = 3.0
+
+[Mesh]
+range = [ -1.0, 1.0 ]
+cfl = 0.10
+n = 7
+k = 201
+basis = "lgl"
+
+[Output]
+every_iteration = 1
+# every_dt        = 0.02
+variables       = [ "rho", "q", "E", "smoothed_mu" ]
+
+[Evolution]
+tend = 0.5
+
+[Log]
+progress_stdout = true
+
+[HRSC]
+method = "av"
+av_smoother_order = 2
+av_method = "mda"
+# av_method = "entropy"
+# entropy_cmax = 2.0
+# entropy_ce = 1.0
+# mda_cmax = 1.0
diff --git a/test/IntegrationTests/refs/fv_euler_isentropic_flow/fv_euler_isentropic_flow.toml b/test/IntegrationTests/refs/fv_euler_isentropic_flow/fv_euler_isentropic_flow.toml
new file mode 100644
index 0000000000000000000000000000000000000000..84a38b5b223c2987edcb9f50405e453cc56a1e25
--- /dev/null
+++ b/test/IntegrationTests/refs/fv_euler_isentropic_flow/fv_euler_isentropic_flow.toml
@@ -0,0 +1,22 @@
+[EquationOfState]
+idealgas_gamma = 3.0
+eos = "idealgas"
+
+[EulerEq]
+id = "smooth_isentropic_flow"
+
+[Evolution]
+cfl = 0.9
+tend = 0.5
+
+[Mesh]
+range = [-1.0, 1.0]
+k = 1048
+scheme = "FV"
+
+
+[Output]
+variables = ["rho", "q", "E"]
+# variables = ["rho"]
+every_dt = 0.02
+# aligned_ts = "$(collect(range(0.01,1.00,step=0.01)))"
diff --git a/test/IntegrationTests/refs/fv_euler_isentropic_flow/output.h5 b/test/IntegrationTests/refs/fv_euler_isentropic_flow/output.h5
new file mode 100644
index 0000000000000000000000000000000000000000..fe9558091764ac5876e6f1479cf25b19eb648e36
Binary files /dev/null and b/test/IntegrationTests/refs/fv_euler_isentropic_flow/output.h5 differ
diff --git a/test/IntegrationTests/refs/testconfig.toml b/test/IntegrationTests/refs/testconfig.toml
index ac56d4c658f45e2cd6e88b526f0bf130903e381c..0a389ca844d02a571298f7b766e4d3565baff6af 100644
--- a/test/IntegrationTests/refs/testconfig.toml
+++ b/test/IntegrationTests/refs/testconfig.toml
@@ -44,6 +44,13 @@ variables = [ "u" ]
 # stops before discontinuity appears
 variables = [ "rho", "q", "E" ]
 
+[fv_euler_sod_shock_tube]
+variables = [ "rho", "q", "E" ]
+
+# [euler_isentropic_flow_regulate-ns]
+# # stops before discontinuity appears
+# variables = [ "rho", "q", "E" ]
+
 [euler_sod_shock_tube_bernstein]
 # we did not emit any outputs -- we are happy when it runs
 variables = [ "rho", "q", "E" ]
@@ -52,7 +59,7 @@ variables = [ "rho", "q", "E" ]
 # [euler_sod_shock_tube_entropyav]
 # variables = [ "rho", "q", "E" ]
 
-[fv_euler_sod_shock_tube]
+[fv_euler_isentropic_flow]
 variables = [ "rho", "q", "E" ]
 
 # TODO Disabled due to refactoring
diff --git a/test/IntegrationTests/src/IntegrationTests.jl b/test/IntegrationTests/src/IntegrationTests.jl
index 85ca4aa5f0989d0f9bff8335f1ca0ed824eccd06..1ee11d5796b0e71544953540797342156934d32b 100644
--- a/test/IntegrationTests/src/IntegrationTests.jl
+++ b/test/IntegrationTests/src/IntegrationTests.jl
@@ -10,6 +10,20 @@ const REFDIR = abspath(joinpath(@__DIR__, "..", "refs"))
 const TESTDIR = abspath(joinpath(@__DIR__, "..", "tests"))
 
 
+function my_findmax(f::Function, args...)
+    max, idx = f(first.(args)...), 1
+    for ii in eachindex(args...)
+        ai = ( a[ii] for a in args )
+        fi = f(ai...)
+        if fi > max
+            max = fi
+            idx = ii
+        end
+    end
+    return max, idx
+end
+
+
 """
     discover_output_dir(outdir)
 
@@ -86,9 +100,18 @@ function compare_outputs(io, testname, refoutput, output,
     if length(refts) != length(ts)
         println(io, """
                     FAIL: Mismatch in number of time steps for which outputs were recorded
-                    reference $(@sprintf("%10d",length(refts))) : $refts
-                    test      $(@sprintf("%10d", length(ts))) : $ts
                     """)
+        nref, ntest= length(refts), length(ts)
+        common_n = max(nref,ntest)
+        data = zeros(Union{Float64,Missing}, common_n, 2)
+        data .= missing
+        data[1:nref,1]  .= refts
+        data[1:ntest,2] .= ts
+        pretty_table(io, data;
+            header = [ :ref, :test ],
+            tf = tf_simple,
+            formatters = ft_printf("% .6e"),
+            crop = :none)
         return false
     end
 
@@ -119,7 +142,9 @@ function compare_outputs(io, testname, refoutput, output,
             continue
         end
 
+        check_failed = false
         upload = false
+        diff_ts, diff_xs, diff_is, diff_maxs = Float64[], Float64[], Float64[], Float64[]
         for i in 1:length(ts)
             refvi = read_variable(rout, var, i)
             vi = read_variable(out, var, i)
@@ -127,11 +152,27 @@ function compare_outputs(io, testname, refoutput, output,
             if !all(zip(refvi, vi)) do (rv, v)
                     isapprox(rv, v, rtol=reltol, atol=abstol)
                 end
-                println(io, "FAIL: Mismatch in variable '$var' @ (t,i) = ($(@sprintf("%10.10e",ts[i])),$i), norm(ref-test,1) = $(norm(refvi.-vi,Inf))")
+                max, idx = my_findmax(refvi,vi) do rv,v
+                    abs(rv-v)
+                end
+                push!(diff_is, i)
+                push!(diff_ts, ts[i])
+                push!(diff_xs, ref_x[idx])
+                push!(diff_maxs, max)
                 success = false
                 upload = true
+                check_failed = true
             end
         end
+        if check_failed
+            println(io, "FAIL: Mismatch in variable '$var'")
+            pretty_table(io, hcat(diff_is, diff_ts, diff_xs, diff_maxs);
+                         header = [ "i", "t[i]", "x", "max(abs(ref-test))" ],
+                         tf = tf_simple,
+                         formatters = ft_printf(["%d", "% .6e", "% .6e", "% .6e"],1:4),
+                         crop = :none)
+            println(io)
+        end
 
         if upload
             diffdir = joinpath(TESTDIR, "diffoutputs")