diff --git a/Manifest.toml b/Manifest.toml new file mode 100644 index 0000000..f911c20 --- /dev/null +++ b/Manifest.toml @@ -0,0 +1,1236 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.8.3" +manifest_format = "2.0" +project_hash = "f2973e3b0b2e44777e1e5ff467ac1c36d72da209" + +[[deps.AMD]] +deps = ["Libdl", "LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "00163dc02b882ca5ec032400b919e5f5011dbd31" +uuid = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e" +version = "0.5.0" + +[[deps.AbstractFFTs]] +deps = ["ChainRulesCore", "LinearAlgebra"] +git-tree-sha1 = "69f7020bd72f069c219b5e8c236c1fa90d2cb409" +uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" +version = "1.2.1" + +[[deps.Adapt]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "0310e08cb19f5da31d08341c6120c047598f5b9c" +uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +version = "3.5.0" + +[[deps.ArgCheck]] +git-tree-sha1 = "a3a402a35a2f7e0b87828ccabbd5ebfbebe356b4" +uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197" +version = "2.3.0" + +[[deps.ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +version = "1.1.1" + +[[deps.ArnoldiMethod]] +deps = ["LinearAlgebra", "Random", "StaticArrays"] +git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae" +uuid = "ec485272-7323-5ecc-a04f-4719b315124d" +version = "0.2.0" + +[[deps.ArrayInterface]] +deps = ["ArrayInterfaceCore", "Compat", "IfElse", "LinearAlgebra", "SnoopPrecompile", "Static"] +git-tree-sha1 = "dedc16cbdd1d32bead4617d27572f582216ccf23" +uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +version = "6.0.25" + +[[deps.ArrayInterfaceCore]] +deps = ["LinearAlgebra", "SnoopPrecompile", "SparseArrays", "SuiteSparse"] +git-tree-sha1 = "e5f08b5689b1aad068e01751889f2f615c7db36d" +uuid = "30b0a656-2188-435a-8636-2ec0e6a096e2" +version = "0.1.29" + +[[deps.ArrayInterfaceGPUArrays]] +deps = ["Adapt", "ArrayInterfaceCore", "GPUArraysCore", "LinearAlgebra"] +git-tree-sha1 = "fc114f550b93d4c79632c2ada2924635aabfa5ed" +uuid = "6ba088a2-8465-4c0a-af30-387133b534db" +version = "0.2.2" + +[[deps.ArrayInterfaceOffsetArrays]] +deps = ["ArrayInterface", "OffsetArrays", "Static"] +git-tree-sha1 = "3d1a9a01976971063b3930d1aed1d9c4af0817f8" +uuid = "015c0d05-e682-4f19-8f0a-679ce4c54826" +version = "0.1.7" + +[[deps.ArrayInterfaceStaticArrays]] +deps = ["Adapt", "ArrayInterface", "ArrayInterfaceCore", "ArrayInterfaceStaticArraysCore", "LinearAlgebra", "Static", "StaticArrays"] +git-tree-sha1 = "f12dc65aef03d0a49650b20b2fdaf184928fd886" +uuid = "b0d46f97-bff5-4637-a19a-dd75974142cd" +version = "0.1.5" + +[[deps.ArrayInterfaceStaticArraysCore]] +deps = ["Adapt", "ArrayInterfaceCore", "LinearAlgebra", "StaticArraysCore"] +git-tree-sha1 = "93c8ba53d8d26e124a5a8d4ec914c3a16e6a0970" +uuid = "dd5226c6-a4d4-4bc7-8575-46859f9c95b9" +version = "0.1.3" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.AutoHashEquals]] +git-tree-sha1 = "45bb6705d93be619b81451bb2006b7ee5d4e4453" +uuid = "15f4f7f2-30c1-5605-9d31-71845cf9641f" +version = "0.2.0" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[deps.BitTwiddlingConvenienceFunctions]] +deps = ["Static"] +git-tree-sha1 = "0c5f81f47bbbcf4aea7b2959135713459170798b" +uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b" +version = "0.1.5" + +[[deps.CEnum]] +git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90" +uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" +version = "0.4.2" + +[[deps.CPUSummary]] +deps = ["CpuId", "IfElse", "Static"] +git-tree-sha1 = "2c144ddb46b552f72d7eafe7cc2f50746e41ea21" +uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9" +version = "0.2.2" + +[[deps.Calculus]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad" +uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" +version = "0.5.1" + +[[deps.ChainRulesCore]] +deps = ["Compat", "LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "c6d890a52d2c4d55d326439580c3b8d0875a77d9" +uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +version = "1.15.7" + +[[deps.ChangesOfVariables]] +deps = ["ChainRulesCore", "LinearAlgebra", "Test"] +git-tree-sha1 = "844b061c104c408b24537482469400af6075aae4" +uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" +version = "0.1.5" + +[[deps.CloseOpenIntervals]] +deps = ["ArrayInterface", "Static"] +git-tree-sha1 = "d61300b9895f129f4bd684b2aff97cf319b6c493" +uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9" +version = "0.1.11" + +[[deps.CodeTracking]] +deps = ["InteractiveUtils", "UUIDs"] +git-tree-sha1 = "0e5c14c3bb8a61b3d53b2c0620570c332c8d0663" +uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" +version = "1.2.0" + +[[deps.CodecZlib]] +deps = ["TranscodingStreams", "Zlib_jll"] +git-tree-sha1 = "9c209fb7536406834aa938fb149964b985de6c83" +uuid = "944b1d66-785c-5afd-91f1-9de20f533193" +version = "0.7.1" + +[[deps.CommonSolve]] +git-tree-sha1 = "9441451ee712d1aec22edad62db1a9af3dc8d852" +uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" +version = "0.2.3" + +[[deps.CommonSubexpressions]] +deps = ["MacroTools", "Test"] +git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" +uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" +version = "0.3.0" + +[[deps.Compat]] +deps = ["Dates", "LinearAlgebra", "UUIDs"] +git-tree-sha1 = "61fdd77467a5c3ad071ef8277ac6bd6af7dd4c04" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.6.0" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "0.5.2+0" + +[[deps.ConstructionBase]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "fb21ddd70a051d882a1686a5a550990bbe371a95" +uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +version = "1.4.1" + +[[deps.CovarianceEstimation]] +deps = ["LinearAlgebra", "Statistics", "StatsBase"] +git-tree-sha1 = "3c8de95b4e932d76ec8960e12d681eba580e9674" +uuid = "587fd27a-f159-11e8-2dae-1979310e6154" +version = "0.2.8" + +[[deps.CpuId]] +deps = ["Markdown"] +git-tree-sha1 = "fcbb72b032692610bfbdb15018ac16a36cf2e406" +uuid = "adafc99b-e345-5852-983c-f28acb93d879" +version = "0.3.1" + +[[deps.Crayons]] +git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15" +uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" +version = "4.1.1" + +[[deps.DataAPI]] +git-tree-sha1 = "e8119c1a33d267e16108be441a287a6981ba1630" +uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" +version = "1.14.0" + +[[deps.DataStructures]] +deps = ["Compat", "InteractiveUtils", "OrderedCollections"] +git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0" +uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +version = "0.18.13" + +[[deps.DataValueInterfaces]] +git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" +uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464" +version = "1.0.0" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[deps.DelimitedFiles]] +deps = ["Mmap"] +uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" + +[[deps.DensityInterface]] +deps = ["InverseFunctions", "Test"] +git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b" +uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d" +version = "0.4.0" + +[[deps.DiffEqBase]] +deps = ["ArrayInterfaceCore", "ChainRulesCore", "DataStructures", "Distributions", "DocStringExtensions", "EnumX", "FastBroadcast", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "Printf", "RecursiveArrayTools", "Reexport", "Requires", "SciMLBase", "Setfield", "SparseArrays", "Static", "StaticArraysCore", "Statistics", "Tricks", "ZygoteRules"] +git-tree-sha1 = "81470904b958f3f24173fa013d4a54198842cb8d" +uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" +version = "6.119.0" + +[[deps.DiffResults]] +deps = ["StaticArraysCore"] +git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621" +uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +version = "1.1.0" + +[[deps.DiffRules]] +deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] +git-tree-sha1 = "c5b6685d53f933c11404a3ae9822afe30d522494" +uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" +version = "1.12.2" + +[[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" + +[[deps.Distributions]] +deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"] +git-tree-sha1 = "9258430c176319dc882efa4088e2ff882a0cb1f1" +uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" +version = "0.25.81" + +[[deps.DocStringExtensions]] +deps = ["LibGit2"] +git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.9.3" + +[[deps.Downloads]] +deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +version = "1.6.0" + +[[deps.DualNumbers]] +deps = ["Calculus", "NaNMath", "SpecialFunctions"] +git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566" +uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" +version = "0.6.8" + +[[deps.EllipsisNotation]] +deps = ["ArrayInterface"] +git-tree-sha1 = "03b753748fd193a7f2730c02d880da27c5a24508" +uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" +version = "1.6.0" + +[[deps.EnumX]] +git-tree-sha1 = "bdb1942cd4c45e3c678fd11569d5cccd80976237" +uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56" +version = "1.0.4" + +[[deps.ExponentialUtilities]] +deps = ["Adapt", "ArrayInterfaceCore", "ArrayInterfaceGPUArrays", "GPUArraysCore", "GenericSchur", "LinearAlgebra", "Printf", "SnoopPrecompile", "SparseArrays", "libblastrampoline_jll"] +git-tree-sha1 = "eb58c1e1417a6580b983069f1491ff82c37def2c" +uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18" +version = "1.23.0" + +[[deps.ExprTools]] +git-tree-sha1 = "56559bbef6ca5ea0c0818fa5c90320398a6fbf8d" +uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" +version = "0.1.8" + +[[deps.FFTW]] +deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] +git-tree-sha1 = "90630efff0894f8142308e334473eba54c433549" +uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +version = "1.5.0" + +[[deps.FFTW_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea" +uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" +version = "3.3.10+0" + +[[deps.FastBroadcast]] +deps = ["ArrayInterface", "ArrayInterfaceCore", "LinearAlgebra", "Polyester", "Static", "StrideArraysCore"] +git-tree-sha1 = "4bef892787c972913d4d84e7255400759bb650e5" +uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +version = "0.2.4" + +[[deps.FastClosures]] +git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef" +uuid = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" +version = "0.3.2" + +[[deps.FastGaussQuadrature]] +deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"] +git-tree-sha1 = "1dfc97ea0df8b3689c2608db8bf28b3fc4f4147b" +uuid = "442a2c76-b920-505d-bb47-c5924d526838" +version = "0.5.0" + +[[deps.FastLapackInterface]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "c1293a93193f0ae94be7cf338d33e162c39d8788" +uuid = "29a986be-02c6-4525-aec4-84b980013641" +version = "1.2.9" + +[[deps.FileWatching]] +uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" + +[[deps.FillArrays]] +deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"] +git-tree-sha1 = "d3ba08ab64bdfd27234d3f61956c966266757fe6" +uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" +version = "0.13.7" + +[[deps.FiniteDiff]] +deps = ["ArrayInterfaceCore", "LinearAlgebra", "Requires", "Setfield", "SparseArrays", "StaticArrays"] +git-tree-sha1 = "04ed1f0029b6b3af88343e439b995141cb0d0b8d" +uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" +version = "2.17.0" + +[[deps.Formatting]] +deps = ["Printf"] +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 = "a69dd6db8a809f78846ff259298678f0d6212180" +uuid = "f6369f11-7733-5829-9624-2563aa707210" +version = "0.10.34" + +[[deps.FunctionWrappers]] +git-tree-sha1 = "d62485945ce5ae9c0c48f124a84998d755bae00e" +uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" +version = "1.1.3" + +[[deps.FunctionWrappersWrappers]] +deps = ["FunctionWrappers"] +git-tree-sha1 = "b104d487b34566608f8b4e1c39fb0b10aa279ff8" +uuid = "77dc65aa-8811-40c2-897b-53d922fa7daf" +version = "0.1.3" + +[[deps.Future]] +deps = ["Random"] +uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" + +[[deps.GPUArraysCore]] +deps = ["Adapt"] +git-tree-sha1 = "1cd7f0af1aa58abc02ea1d872953a97359cb87fa" +uuid = "46192b85-c4d5-4398-a991-12ede77f4527" +version = "0.1.4" + +[[deps.GaussQuadrature]] +deps = ["SpecialFunctions"] +git-tree-sha1 = "eb6f1f48aa994f3018cbd029a17863c6535a266d" +uuid = "d54b0c1a-921d-58e0-8e36-89d8069c0969" +version = "0.5.8" + +[[deps.GenericSchur]] +deps = ["LinearAlgebra", "Printf"] +git-tree-sha1 = "fb69b2a645fa69ba5f474af09221b9308b160ce6" +uuid = "c145ed77-6b09-5dd9-b285-bf645a82121e" +version = "0.5.3" + +[[deps.Graphs]] +deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] +git-tree-sha1 = "1cf1d7dcb4bc32d7b4a5add4232db3750c27ecb4" +uuid = "86223c79-3864-5bf0-83f7-82e725a168b6" +version = "1.8.0" + +[[deps.HDF5]] +deps = ["Compat", "HDF5_jll", "Libdl", "Mmap", "Random", "Requires", "UUIDs"] +git-tree-sha1 = "b5df7c3cab3a00c33c2e09c6bd23982a75e2fbb2" +uuid = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" +version = "0.16.13" + +[[deps.HDF5_jll]] +deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenSSL_jll", "Pkg", "Zlib_jll"] +git-tree-sha1 = "4cc2bb72df6ff40b055295fdef6d92955f9dede8" +uuid = "0234f1f7-429e-5d53-9886-15a909be8d59" +version = "1.12.2+2" + +[[deps.HostCPUFeatures]] +deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"] +git-tree-sha1 = "734fd90dd2f920a2f1921d5388dcebe805b262dc" +uuid = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0" +version = "0.1.14" + +[[deps.HypergeometricFunctions]] +deps = ["DualNumbers", "LinearAlgebra", "OpenLibm_jll", "SpecialFunctions", "Test"] +git-tree-sha1 = "709d864e3ed6e3545230601f94e11ebc65994641" +uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" +version = "0.3.11" + +[[deps.IfElse]] +git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" +uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" +version = "0.1.1" + +[[deps.Inflate]] +git-tree-sha1 = "5cd07aab533df5170988219191dfad0519391428" +uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" +version = "0.1.3" + +[[deps.IntelOpenMP_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c" +uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" +version = "2018.0.3+2" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[deps.InverseFunctions]] +deps = ["Test"] +git-tree-sha1 = "49510dfcb407e572524ba94aeae2fced1f3feb0f" +uuid = "3587e190-3f89-42d0-90ee-14403ec27112" +version = "0.1.8" + +[[deps.IrrationalConstants]] +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" +version = "1.0.0" + +[[deps.JLLWrappers]] +deps = ["Preferences"] +git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.4.1" + +[[deps.KLU]] +deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse_jll"] +git-tree-sha1 = "764164ed65c30738750965d55652db9c94c59bfe" +uuid = "ef3ab10e-7fda-4108-b977-705223b18434" +version = "0.4.0" + +[[deps.Kronecker]] +deps = ["LinearAlgebra", "NamedDims", "SparseArrays", "StatsBase"] +git-tree-sha1 = "78d9909daf659c901ae6c7b9de7861ba45a743f4" +uuid = "2c470bb0-bcc8-11e8-3dad-c9649493f05e" +version = "0.5.3" + +[[deps.Krylov]] +deps = ["LinearAlgebra", "Printf", "SparseArrays"] +git-tree-sha1 = "dd90aacbfb622f898a97c2a4411ac49101ebab8a" +uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" +version = "0.9.0" + +[[deps.KrylovKit]] +deps = ["ChainRulesCore", "GPUArraysCore", "LinearAlgebra", "Printf"] +git-tree-sha1 = "1a5e1d9941c783b0119897d29f2eb665d876ecf3" +uuid = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" +version = "0.6.0" + +[[deps.LDLFactorizations]] +deps = ["AMD", "LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "cbf4b646f82bfc58bb48bcca9dcce2eb88da4cd1" +uuid = "40e66cde-538c-5869-a4ad-c39174c6795b" +version = "0.10.0" + +[[deps.LaTeXStrings]] +git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996" +uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +version = "1.3.0" + +[[deps.LayoutPointers]] +deps = ["ArrayInterface", "ArrayInterfaceOffsetArrays", "ArrayInterfaceStaticArrays", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static"] +git-tree-sha1 = "0ad6f0c51ce004dadc24a28a0dfecfb568e52242" +uuid = "10f19ff3-798f-405d-979b-55457f8fc047" +version = "0.1.13" + +[[deps.Lazy]] +deps = ["MacroTools"] +git-tree-sha1 = "1370f8202dac30758f3c345f9909b97f53d87d3f" +uuid = "50d2b5c4-7a5e-59d5-8109-a42b560f39c0" +version = "0.15.1" + +[[deps.LazyArtifacts]] +deps = ["Artifacts", "Pkg"] +uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" + +[[deps.LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.6.3" + +[[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" + +[[deps.LibGit2]] +deps = ["Base64", "NetworkOptions", "Printf", "SHA"] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[deps.LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" +version = "1.10.2+0" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[deps.Libiconv_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c7cb1f5d892775ba13767a87c7ada0b980ea0a71" +uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" +version = "1.16.1+2" + +[[deps.LightXML]] +deps = ["Libdl", "XML2_jll"] +git-tree-sha1 = "e129d9391168c677cd4800f5c0abb1ed8cb3794f" +uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179" +version = "0.9.0" + +[[deps.LineSearches]] +deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"] +git-tree-sha1 = "7bbea35cec17305fc70a0e5b4641477dc0789d9d" +uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +version = "7.2.0" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[deps.LinearMaps]] +deps = ["ChainRulesCore", "LinearAlgebra", "SparseArrays", "Statistics"] +git-tree-sha1 = "42970dad6b0d2515571613010bd32ba37e07f874" +uuid = "7a12625a-238d-50fd-b39a-03d52299707e" +version = "3.9.0" + +[[deps.LinearOperators]] +deps = ["FastClosures", "LDLFactorizations", "LinearAlgebra", "Printf", "SparseArrays", "TimerOutputs"] +git-tree-sha1 = "e445a26217944cc1597e2879c604bd6b47bd7c78" +uuid = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125" +version = "2.5.1" + +[[deps.LinearSolve]] +deps = ["ArrayInterface", "DocStringExtensions", "FastLapackInterface", "GPUArraysCore", "IterativeSolvers", "KLU", "Krylov", "KrylovKit", "LinearAlgebra", "Preferences", "RecursiveFactorization", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "SnoopPrecompile", "SparseArrays", "Sparspak", "SuiteSparse", "UnPack"] +git-tree-sha1 = "d1fce810e9a4213607f0182cf25ffd6ce13e19b6" +uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +version = "1.37.0" + +[[deps.LogExpFunctions]] +deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "680e733c3a0a9cea9e935c8c2184aea6a63fa0b5" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "0.3.21" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[deps.LoopVectorization]] +deps = ["ArrayInterface", "ArrayInterfaceCore", "ArrayInterfaceOffsetArrays", "ArrayInterfaceStaticArrays", "CPUSummary", "ChainRulesCore", "CloseOpenIntervals", "DocStringExtensions", "ForwardDiff", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "SIMDTypes", "SLEEFPirates", "SnoopPrecompile", "SpecialFunctions", "Static", "ThreadingUtilities", "UnPack", "VectorizationBase"] +git-tree-sha1 = "9696a80c21a56b937e3fd89e972f8db5db3186e2" +uuid = "bdcacae8-1622-11e9-2a5c-532679323890" +version = "0.12.150" + +[[deps.MKL_jll]] +deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] +git-tree-sha1 = "2ce8695e1e699b68702c03402672a69f54b8aca9" +uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" +version = "2022.2.0+0" + +[[deps.MPI]] +deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "Requires", "Serialization", "Sockets"] +git-tree-sha1 = "6d72bafd3960f9c119ceb8f034fef28346490fe5" +uuid = "da04e1cc-30fd-572f-bb4f-1f8673147195" +version = "0.20.8" + +[[deps.MPICH_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "Pkg", "TOML"] +git-tree-sha1 = "6d4fa43afab4611d090b11617ecea1a144b21d35" +uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4" +version = "4.0.2+5" + +[[deps.MPIPreferences]] +deps = ["Libdl", "Preferences"] +git-tree-sha1 = "71f937129731a29eabe6969db2c90368a4408933" +uuid = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267" +version = "0.1.7" + +[[deps.MPItrampoline_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "Pkg", "TOML"] +git-tree-sha1 = "b3f9e42685b4ad614eca0b44bd863cd41b1c86ea" +uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748" +version = "5.0.2+1" + +[[deps.MacroTools]] +deps = ["Markdown", "Random"] +git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.10" + +[[deps.ManualMemory]] +git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd" +uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667" +version = "0.1.8" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[deps.MbedTLS_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" +version = "2.28.0+0" + +[[deps.MicrosoftMPI_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "a16aa086d335ed7e0170c5265247db29172af2f9" +uuid = "9237b28f-5490-5468-be7b-bb81f5f5e6cf" +version = "10.1.3+2" + +[[deps.Missings]] +deps = ["DataAPI"] +git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272" +uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" +version = "1.1.0" + +[[deps.Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[deps.MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +version = "2022.2.1" + +[[deps.MuladdMacro]] +git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab" +uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +version = "0.2.4" + +[[deps.NLSolversBase]] +deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] +git-tree-sha1 = "a0b464d183da839699f4c79e7606d9d186ec172c" +uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" +version = "7.8.3" + +[[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]] +deps = ["OpenLibm_jll"] +git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "1.0.2" + +[[deps.NamedDims]] +deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"] +git-tree-sha1 = "cb8ebcee2b4e07b72befb9def593baef8aa12f07" +uuid = "356022a1-0364-5f58-8944-0da4b18d706f" +version = "0.2.50" + +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" + +[[deps.NodesAndModes]] +deps = ["DelimitedFiles", "LinearAlgebra", "SpecialFunctions"] +git-tree-sha1 = "db18627df685c339878b6f63a31addd3806c8e76" +uuid = "7aca2e03-f7e2-4192-9ec8-f4ca66d597fb" +version = "0.9.0" + +[[deps.NonlinearSolve]] +deps = ["ArrayInterface", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "LinearSolve", "RecursiveArrayTools", "Reexport", "SciMLBase", "SimpleNonlinearSolve", "SnoopPrecompile", "SparseArrays", "SparseDiffTools", "StaticArraysCore", "UnPack"] +git-tree-sha1 = "536aa8b33b2c3a10df8ce89bdb0b0affef93d393" +uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +version = "1.4.0" + +[[deps.Octavian]] +deps = ["ArrayInterface", "CPUSummary", "IfElse", "LoopVectorization", "ManualMemory", "PolyesterWeave", "Requires", "SnoopPrecompile", "Static", "ThreadingUtilities", "VectorizationBase"] +git-tree-sha1 = "b6c8c7f574c981546cf7b0f017572e5d019991a3" +uuid = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" +version = "0.3.20" + +[[deps.OffsetArrays]] +deps = ["Adapt"] +git-tree-sha1 = "82d7c9e310fe55aa54996e6f7f94674e2a38fcb4" +uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +version = "1.12.9" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.20+0" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.1+0" + +[[deps.OpenMPI_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "Pkg", "TOML"] +git-tree-sha1 = "346d6b357a480300ed7854dbc70e746ac52e10fd" +uuid = "fe0851c0-eecd-5654-98d4-656369965a5c" +version = "4.1.3+3" + +[[deps.OpenSSL_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "9ff31d101d987eb9d66bd8b176ac7c277beccd09" +uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" +version = "1.1.20+0" + +[[deps.OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.5+0" + +[[deps.OrderedCollections]] +git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.4.1" + +[[deps.OrdinaryDiffEq]] +deps = ["Adapt", "ArrayInterface", "ArrayInterfaceCore", "ArrayInterfaceGPUArrays", "ArrayInterfaceStaticArrays", "ArrayInterfaceStaticArraysCore", "DataStructures", "DiffEqBase", "DocStringExtensions", "ExponentialUtilities", "FastBroadcast", "FastClosures", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "Logging", "LoopVectorization", "MacroTools", "MuladdMacro", "NLsolve", "NonlinearSolve", "Polyester", "PreallocationTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLNLSolve", "SimpleNonlinearSolve", "SnoopPrecompile", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"] +git-tree-sha1 = "295b13e65001f8e08fb6a861bdf61ade8ffc460b" +uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +version = "6.44.1" + +[[deps.P4est]] +deps = ["CEnum", "MPI", "MPIPreferences", "P4est_jll", "Preferences", "Reexport", "UUIDs"] +git-tree-sha1 = "326049aeefa844ba03ab12a57857d04a2f60d8df" +uuid = "7d669430-f675-4ae7-b43e-fab78ec5a902" +version = "0.4.0" + +[[deps.P4est_jll]] +deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "Pkg", "TOML", "Zlib_jll"] +git-tree-sha1 = "d4b48fd3ca75a398916c58c1e4628bf0ce11a7b6" +uuid = "6b5a15aa-cf52-5330-8376-5e5d90283449" +version = "2.8.1+0" + +[[deps.PDMats]] +deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] +git-tree-sha1 = "cf494dca75a69712a72b80bc48f59dcf3dea63ec" +uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" +version = "0.11.16" + +[[deps.Parameters]] +deps = ["OrderedCollections", "UnPack"] +git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe" +uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a" +version = "0.12.3" + +[[deps.PathIntersections]] +deps = ["ForwardDiff", "GaussQuadrature", "LinearAlgebra", "SparseArrays", "StaticArrays", "StructArrays"] +git-tree-sha1 = "5957f6a65b3bb31004a753402d9da24b264244cd" +uuid = "4c1a95c7-462a-4a7e-b284-959c63fbf1dc" +version = "0.1.1" + +[[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" +version = "1.8.0" + +[[deps.Polyester]] +deps = ["ArrayInterface", "BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "ManualMemory", "PolyesterWeave", "Requires", "Static", "StrideArraysCore", "ThreadingUtilities"] +git-tree-sha1 = "e8e0fabcff4df8686c4267503887202a783d498e" +uuid = "f517fe37-dbe3-4b94-8317-1923a5111588" +version = "0.7.2" + +[[deps.PolyesterWeave]] +deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"] +git-tree-sha1 = "240d7170f5ffdb285f9427b92333c3463bf65bf6" +uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad" +version = "0.2.1" + +[[deps.PolynomialBases]] +deps = ["ArgCheck", "AutoHashEquals", "FFTW", "FastGaussQuadrature", "LinearAlgebra", "Requires", "SpecialFunctions", "UnPack"] +git-tree-sha1 = "38629c0a9cace7c6f51c103084589ff8a7a1c02f" +uuid = "c74db56a-226d-5e98-8bb0-a6049094aeea" +version = "0.4.15" + +[[deps.PreallocationTools]] +deps = ["Adapt", "ArrayInterfaceCore", "ForwardDiff", "Requires"] +git-tree-sha1 = "2c7658dd593e3adc118b00429e1048829f1abb8c" +uuid = "d236fae5-4411-538c-8e31-a6e3d9e00b46" +version = "0.4.11" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.3.0" + +[[deps.PrettyTables]] +deps = ["Crayons", "Formatting", "LaTeXStrings", "Markdown", "Reexport", "StringManipulation", "Tables"] +git-tree-sha1 = "96f6db03ab535bdb901300f88335257b0018689d" +uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +version = "2.2.2" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[deps.QuadGK]] +deps = ["DataStructures", "LinearAlgebra"] +git-tree-sha1 = "786efa36b7eff813723c4849c90456609cf06661" +uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +version = "2.8.1" + +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[deps.Random]] +deps = ["SHA", "Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.RecipesBase]] +deps = ["SnoopPrecompile"] +git-tree-sha1 = "261dddd3b862bd2c940cf6ca4d1c8fe593e457c8" +uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +version = "1.3.3" + +[[deps.RecursiveArrayTools]] +deps = ["Adapt", "ArrayInterfaceCore", "ArrayInterfaceStaticArraysCore", "ChainRulesCore", "DocStringExtensions", "FillArrays", "GPUArraysCore", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "Requires", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables", "ZygoteRules"] +git-tree-sha1 = "54e055256bbd41fd10566880bc4baa5316bca6fe" +uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" +version = "2.37.0" + +[[deps.RecursiveFactorization]] +deps = ["LinearAlgebra", "LoopVectorization", "Polyester", "SnoopPrecompile", "StrideArraysCore", "TriangularSolve"] +git-tree-sha1 = "9088515ad915c99026beb5436d0a09cd8c18163e" +uuid = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" +version = "0.2.18" + +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "1.2.2" + +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.0" + +[[deps.Rmath]] +deps = ["Random", "Rmath_jll"] +git-tree-sha1 = "f65dcb5fa46aee0cf9ed6274ccbd597adc49aa7b" +uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa" +version = "0.7.1" + +[[deps.Rmath_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "6ed52fdd3382cf21947b15e8870ac0ddbff736da" +uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f" +version = "0.4.0+0" + +[[deps.RuntimeGeneratedFunctions]] +deps = ["ExprTools", "SHA", "Serialization"] +git-tree-sha1 = "50314d2ef65fce648975a8e80ae6d8409ebbf835" +uuid = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" +version = "0.5.5" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.SIMDTypes]] +git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c" +uuid = "94e857df-77ce-4151-89e5-788b33177be4" +version = "0.1.0" + +[[deps.SLEEFPirates]] +deps = ["IfElse", "Static", "VectorizationBase"] +git-tree-sha1 = "cda0aece8080e992f6370491b08ef3909d1c04e7" +uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa" +version = "0.6.38" + +[[deps.SciMLBase]] +deps = ["ArrayInterfaceCore", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "Preferences", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"] +git-tree-sha1 = "76eec814289c4a249ee3747ceeea0d83defbeb8d" +uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +version = "1.84.1" + +[[deps.SciMLNLSolve]] +deps = ["DiffEqBase", "LineSearches", "NLsolve", "Reexport", "SciMLBase"] +git-tree-sha1 = "66c7f901dbcad51791136e2d90ee67240256ecde" +uuid = "e9a6253c-8580-4d32-9898-8661bb511710" +version = "0.1.3" + +[[deps.SciMLOperators]] +deps = ["ArrayInterfaceCore", "DocStringExtensions", "Lazy", "LinearAlgebra", "Setfield", "SparseArrays", "StaticArraysCore", "Tricks"] +git-tree-sha1 = "c737d575c18bdf9aba0a3c7071d5249d09f45dd8" +uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" +version = "0.1.21" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.Setfield]] +deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"] +git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac" +uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" +version = "1.1.1" + +[[deps.SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + +[[deps.SimpleNonlinearSolve]] +deps = ["ArrayInterface", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "Reexport", "Requires", "SciMLBase", "SnoopPrecompile", "StaticArraysCore"] +git-tree-sha1 = "326789bbaa1b65b809bd4596b74e4fc3be5af6ac" +uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7" +version = "0.1.13" + +[[deps.SimpleTraits]] +deps = ["InteractiveUtils", "MacroTools"] +git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231" +uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d" +version = "0.9.4" + +[[deps.SnoopPrecompile]] +deps = ["Preferences"] +git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c" +uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c" +version = "1.0.3" + +[[deps.Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[deps.SortingAlgorithms]] +deps = ["DataStructures"] +git-tree-sha1 = "a4ada03f999bd01b3a25dcaa30b2d929fe537e00" +uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" +version = "1.1.0" + +[[deps.SparseArrays]] +deps = ["LinearAlgebra", "Random"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.SparseDiffTools]] +deps = ["Adapt", "ArrayInterfaceCore", "ArrayInterfaceStaticArrays", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "Graphs", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays", "VertexSafeGraphs"] +git-tree-sha1 = "4245283bee733122a9cb4545748d64e0c63337c0" +uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804" +version = "1.30.0" + +[[deps.Sparspak]] +deps = ["Libdl", "LinearAlgebra", "Logging", "OffsetArrays", "Printf", "SparseArrays", "Test"] +git-tree-sha1 = "342cf4b449c299d8d1ceaf00b7a49f4fbc7940e7" +uuid = "e56a9233-b9d6-4f03-8d0f-1825330902ac" +version = "0.3.9" + +[[deps.SpecialFunctions]] +deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "d75bda01f8c31ebb72df80a46c88b25d1c79c56d" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "2.1.7" + +[[deps.StartUpDG]] +deps = ["ConstructionBase", "FillArrays", "HDF5", "Kronecker", "LinearAlgebra", "NodesAndModes", "OrderedCollections", "PathIntersections", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "Requires", "Setfield", "SparseArrays", "StaticArrays", "Triangulate", "UnPack", "WriteVTK"] +git-tree-sha1 = "baca9a9af25821f26f34d5489f7dc358732e623f" +uuid = "472ebc20-7c99-4d4b-9470-8fde4e9faa0f" +version = "0.15.9" + +[[deps.Static]] +deps = ["IfElse"] +git-tree-sha1 = "c35b107b61e7f34fa3f124026f2a9be97dea9e1c" +uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" +version = "0.8.3" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"] +git-tree-sha1 = "67d3e75e8af8089ea34ce96974d5468d4a008ca6" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.5.15" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.0" + +[[deps.Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[deps.StatsAPI]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "f9af7f195fb13589dd2e2d57fdb401717d2eb1f6" +uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" +version = "1.5.0" + +[[deps.StatsBase]] +deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] +git-tree-sha1 = "d1bf48bfcc554a3761a133fe3a9bb01488e06916" +uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +version = "0.33.21" + +[[deps.StatsFuns]] +deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] +git-tree-sha1 = "ab6083f09b3e617e34a956b43e9d51b824206932" +uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" +version = "1.1.1" + +[[deps.StrideArrays]] +deps = ["ArrayInterface", "LinearAlgebra", "LoopVectorization", "Octavian", "Random", "SLEEFPirates", "Static", "StaticArraysCore", "Statistics", "StrideArraysCore", "VectorizationBase", "VectorizedRNG", "VectorizedStatistics"] +git-tree-sha1 = "06a592d96244200bfcc31399e9c436a41331a1c2" +uuid = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b" +version = "0.1.23" + +[[deps.StrideArraysCore]] +deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "ManualMemory", "SIMDTypes", "Static", "ThreadingUtilities"] +git-tree-sha1 = "8114ba9c3694827838d45ea3c9f6b9ccb4182cf2" +uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da" +version = "0.4.7" + +[[deps.StringManipulation]] +git-tree-sha1 = "46da2434b41f41ac3594ee9816ce5541c6096123" +uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e" +version = "0.3.0" + +[[deps.StructArrays]] +deps = ["Adapt", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"] +git-tree-sha1 = "b03a3b745aa49b566f128977a7dd1be8711c5e71" +uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" +version = "0.6.14" + +[[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" +version = "5.10.1+0" + +[[deps.SummationByPartsOperators]] +deps = ["ArgCheck", "ArrayInterface", "AutoHashEquals", "FFTW", "InteractiveUtils", "LinearAlgebra", "LoopVectorization", "MuladdMacro", "PolynomialBases", "RecursiveArrayTools", "Reexport", "Requires", "SciMLBase", "SnoopPrecompile", "SparseArrays", "StaticArrays", "UnPack", "Unrolled"] +git-tree-sha1 = "72df7ed25cbc75a0d059f11e9fbc07b2a8fef3e8" +uuid = "9f78cca6-572e-554e-b819-917d2f1cf240" +version = "0.5.28" + +[[deps.SymbolicIndexingInterface]] +deps = ["DocStringExtensions"] +git-tree-sha1 = "6b764c160547240d868be4e961a5037f47ad7379" +uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5" +version = "0.2.1" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.0" + +[[deps.TableTraits]] +deps = ["IteratorInterfaceExtensions"] +git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39" +uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" +version = "1.0.1" + +[[deps.Tables]] +deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"] +git-tree-sha1 = "c79322d36826aa2f4fd8ecfa96ddb47b174ac78d" +uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +version = "1.10.0" + +[[deps.Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" +version = "1.10.1" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[deps.ThreadingUtilities]] +deps = ["ManualMemory"] +git-tree-sha1 = "c97f60dd4f2331e1a495527f80d242501d2f9865" +uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5" +version = "0.5.1" + +[[deps.TimerOutputs]] +deps = ["ExprTools", "Printf"] +git-tree-sha1 = "f2fd3f288dfc6f507b0c3a2eb3bac009251e548b" +uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +version = "0.5.22" + +[[deps.TranscodingStreams]] +deps = ["Random", "Test"] +git-tree-sha1 = "94f38103c984f89cf77c402f2a68dbd870f8165f" +uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" +version = "0.9.11" + +[[deps.Triangle_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "fe28e9a4684f6f54e868b9136afb8fd11f1734a7" +uuid = "5639c1d2-226c-5e70-8d55-b3095415a16a" +version = "1.6.2+0" + +[[deps.TriangularSolve]] +deps = ["CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "LoopVectorization", "Polyester", "Static", "VectorizationBase"] +git-tree-sha1 = "31eedbc0b6d07c08a700e26d31298ac27ef330eb" +uuid = "d5829a12-d9aa-46ab-831f-fb7c9ab06edf" +version = "0.1.19" + +[[deps.Triangulate]] +deps = ["DocStringExtensions", "Libdl", "Printf", "Test", "Triangle_jll"] +git-tree-sha1 = "bbca6ec35426334d615f58859ad40c96d3a4a1f9" +uuid = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" +version = "2.2.0" + +[[deps.Tricks]] +git-tree-sha1 = "6bac775f2d42a611cdfcd1fb217ee719630c4175" +uuid = "410a4b4d-49e4-4fbc-ab6d-cb71b17b3775" +version = "0.1.6" + +[[deps.TriplotBase]] +git-tree-sha1 = "4d4ed7f294cda19382ff7de4c137d24d16adc89b" +uuid = "981d1d27-644d-49a2-9326-4793e63143c3" +version = "0.1.0" + +[[deps.TriplotRecipes]] +deps = ["RecipesBase", "TriplotBase"] +git-tree-sha1 = "fceb3b0f37ff6ccf3c70b9c5198d2eefec46ada0" +uuid = "808ab39a-a642-4abf-81ff-4cb34ebbffa3" +version = "0.1.2" + +[[deps.Trixi]] +deps = ["CodeTracking", "ConstructionBase", "EllipsisNotation", "FillArrays", "ForwardDiff", "HDF5", "IfElse", "LinearAlgebra", "LinearMaps", "LoopVectorization", "MPI", "MuladdMacro", "Octavian", "OffsetArrays", "P4est", "Polyester", "Printf", "RecipesBase", "Reexport", "Requires", "SciMLBase", "Setfield", "SparseArrays", "StartUpDG", "Static", "StaticArrays", "StrideArrays", "StructArrays", "SummationByPartsOperators", "TimerOutputs", "Triangulate", "TriplotBase", "TriplotRecipes", "UnPack"] +git-tree-sha1 = "8ce74d9676f3ce1cee01028d0418e52f53e9b016" +uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" +version = "0.5.12" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[deps.UnPack]] +git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b" +uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" +version = "1.0.2" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[deps.Unrolled]] +deps = ["MacroTools"] +git-tree-sha1 = "69eccecbb9ba198e46b4ddcb164f8dea3c685601" +uuid = "9602ed7d-8fef-5bc8-8597-8f21381861e8" +version = "0.1.4" + +[[deps.VectorizationBase]] +deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static"] +git-tree-sha1 = "4c59c2df8d2676c4691a39fa70495a6db0c5d290" +uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f" +version = "0.21.58" + +[[deps.VectorizedRNG]] +deps = ["Distributed", "Random", "SLEEFPirates", "StaticArraysCore", "StrideArraysCore", "UnPack", "VectorizationBase"] +git-tree-sha1 = "d2032df5b000bfa871e1a6a4112a399540d4f935" +uuid = "33b4df10-0173-11e9-2a0c-851a7edac40e" +version = "0.2.23" + +[[deps.VectorizedStatistics]] +deps = ["LoopVectorization", "SnoopPrecompile", "Static"] +git-tree-sha1 = "d952a2ec061cf1eb79b448fec9cde0a403643213" +uuid = "3b853605-1c98-4422-8364-4bd93ee0529e" +version = "0.5.2" + +[[deps.VertexSafeGraphs]] +deps = ["Graphs"] +git-tree-sha1 = "8351f8d73d7e880bfc042a8b6922684ebeafb35c" +uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" +version = "0.2.0" + +[[deps.WriteVTK]] +deps = ["Base64", "CodecZlib", "FillArrays", "LightXML", "TranscodingStreams"] +git-tree-sha1 = "49353f30da65f377cff0f934bb9f562a2c0441b9" +uuid = "64499a7a-5c06-52f2-abe2-ccb03c286192" +version = "1.17.1" + +[[deps.XML2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"] +git-tree-sha1 = "93c41695bc1c08c46c5899f4fe06d6ead504bb73" +uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" +version = "2.10.3+0" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.12+3" + +[[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" +version = "5.1.1+0" + +[[deps.nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" +version = "1.48.0+0" + +[[deps.p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" +version = "17.4.0+0" diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..c18ff45 --- /dev/null +++ b/Project.toml @@ -0,0 +1,12 @@ +[deps] +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" +Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" +SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" +Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" diff --git a/README.md b/README.md index bc155c4..a7db4d7 100644 --- a/README.md +++ b/README.md @@ -55,10 +55,20 @@ to install [Julia](https://julialang.org/). The numerical experiments presented in this article were performed using Julia v1.8.3. First, you need to download this repository, e.g., by cloning it with `git` -or by downloading an archive via the GitHub interface. +or by downloading an archive via the GitHub interface. Then, you need to start +Julia in this directory and execute the following commands in the Julia REPL. -Then, you need to start Julia in the `code` directory and follow the -instructions in the `README.md` file therein. +```julia +julia> include("code.jl") + +julia> convergence_tests_1d() + +julia> convergence_tests_2d() +``` + +This will show the results in the REPL. If you want to display LaTeX code for +the convergence tables, add the keyword argument `latex = true` to the function +calls shown above. ## Authors diff --git a/code.jl b/code.jl new file mode 100644 index 0000000..b39b841 --- /dev/null +++ b/code.jl @@ -0,0 +1,1315 @@ +# Setup +import Pkg +Pkg.activate(@__DIR__) +Pkg.instantiate() + + +# Load packages +using LinearAlgebra +using Printf +using SparseArrays + +using SpecialFunctions: erf + +using SummationByPartsOperators +using Trixi +using OrdinaryDiffEq: CallbackSet + +using PrettyTables: pretty_table, ft_printf +using FillArrays: Zeros +using LinearOperators: LinearOperator +using Krylov: Krylov + +BLAS.set_num_threads(1) + + +# helper functions creating a (sparse) matrix representation from a matrix-free +# implementation in `f!` +function compute_linear_structure(f!, tmp, params, t) + res = similar(tmp) + + tmp .= zero(eltype(tmp)) + f!(res, tmp, params, t) + b = -vec(res) + + A = zeros(eltype(res), length(res), length(res)) + for j in 1:length(res) + tmp[j] = one(eltype(tmp)) + f!(res, tmp, params, t) + A[:, j] .= vec(res) .+ b + tmp[j] = zero(eltype(tmp)) + end + + return A, b +end + +function compute_sparse_linear_structure(f!, tmp, params, t) + res = similar(tmp) + vec_res = vec(res) + + tmp .= zero(eltype(tmp)) + f!(res, tmp, params, t) + b = -vec(res) + + rowind = Vector{Int}() + nzval = Vector{eltype(tmp)}() + colptr = Vector{Int}(undef, length(tmp) + 1) + + for i in 1:length(res) + tmp[i] = one(eltype(tmp)) + f!(res, tmp, params, t) + @. vec_res = vec_res + b + js = findall(!iszero, vec_res) + colptr[i] = length(nzval) + 1 + if length(js) > 0 + append!(rowind, js) + append!(nzval, vec_res[js]) + end + tmp[i] = zero(eltype(tmp)) + end + colptr[end] = length(nzval) + 1 + + return SparseMatrixCSC(length(tmp), length(tmp), colptr, rowind, nzval), b +end + + +# 1D functions +function compute_q!(q, φ::AbstractArray{<:Real, 2}, params) + (; basis, dl, dr, jacobian, Tᵣ, surface_values, φBC_left, φBC_right) = params + sqrtTᵣ = sqrt(Tᵣ) + + # prolong φ, Dφ to surfaces + for element in axes(φ, 2) + left, right = element, element + 1 + inv_jacobian = inv(jacobian[element]) + surface_values[1, 2, left ] = φ[1, element] + surface_values[2, 2, left ] = dot(dl, view(φ, :, element)) * inv_jacobian + surface_values[1, 1, right] = φ[end, element] + surface_values[2, 1, right] = dot(dr, view(φ, :, element)) * inv_jacobian + end + + # compute q locally + mul!(q, basis.D, φ) + for element in axes(φ, 2) + left, right = element, element + 1 + inv_jacobian = inv(jacobian[element]) + @. q[:, element] *= inv_jacobian + + # left boundary correction + if element == 1 + if φBC_left === :periodic + jump_φ = surface_values[1, 2, left] - surface_values[1, 1, end] + jump_Dφ = surface_values[2, 2, left] - surface_values[2, 1, end] + inv_weight_r = inv(jacobian[element ] * basis.weights[1]) + inv_weight_l = inv(jacobian[end ] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q[1, element] += inv_jacobian / basis.weights[1] * ( + 0.5 * (1 + sqrtTᵣ * c1) * jump_φ - 0.5 * sqrtTᵣ * c2 * jump_Dφ) + else + jump_φ = surface_values[1, 2, left] - φBC_left + q[1, element] += inv_jacobian / basis.weights[1] * jump_φ + end + else + jump_φ = surface_values[1, 2, left] - surface_values[1, 1, left] + jump_Dφ = surface_values[2, 2, left] - surface_values[2, 1, left] + inv_weight_r = inv(jacobian[element ] * basis.weights[1]) + inv_weight_l = inv(jacobian[element-1] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q[1, element] += inv_jacobian / basis.weights[1] * ( + 0.5 * (1 + sqrtTᵣ * c1) * jump_φ - 0.5 * sqrtTᵣ * c2 * jump_Dφ) + end + + # right boundary correction + if element == size(φ, 2) + if φBC_right === :periodic + jump_φ = surface_values[1, 2, 1] - surface_values[1, 1, right] + jump_Dφ = surface_values[2, 2, 1] - surface_values[2, 1, right] + inv_weight_r = inv(jacobian[1 ] * basis.weights[1]) + inv_weight_l = inv(jacobian[element ] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q[end, element] += inv_jacobian / basis.weights[end] * ( + 0.5 * (1 - sqrtTᵣ * c1) * jump_φ + 0.5 * sqrtTᵣ * c2 * jump_Dφ) + else + jump_φ = φBC_right - surface_values[1, 1, right] + q[end, element] += inv_jacobian / basis.weights[end] * jump_φ + end + else + jump_φ = surface_values[1, 2, right] - surface_values[1, 1, right] + jump_Dφ = surface_values[2, 2, right] - surface_values[2, 1, right] + inv_weight_r = inv(jacobian[element+1] * basis.weights[1]) + inv_weight_l = inv(jacobian[element ] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q[end, element] += inv_jacobian / basis.weights[end] * ( + 0.5 * (1 - sqrtTᵣ * c1) * jump_φ + 0.5 * sqrtTᵣ * c2 * jump_Dφ) + end + end + + return nothing +end + +function compute_div_q!(div_q, φ::AbstractArray{<:Real, 2}, q, params) + (; basis, jacobian, Tᵣ, surface_values, φBC_left, φBC_right) = params + sqrtTᵣ = sqrt(Tᵣ) + inv_sqrtTᵣ = inv(sqrtTᵣ) + + # prolong φ, q to surfaces + for element in axes(φ, 2) + left, right = element, element+1 + surface_values[1, 2, left ] = φ[1, element] + surface_values[2, 2, left ] = q[1, element] + surface_values[1, 1, right] = φ[end, element] + surface_values[2, 1, right] = q[end, element] + end + + # compute div q locally + mul!(div_q, basis.D, q) + for element in axes(φ, 2) + left, right = element, element+1 + inv_jacobian = inv(jacobian[element]) + @. div_q[:, element] *= inv_jacobian + + # left boundary correction + if element == 1 + if φBC_left === :periodic + jump_φ = surface_values[1, 2, left] - surface_values[1, 1, end] + jump_q = surface_values[2, 2, left] - surface_values[2, 1, end] + div_q[1, element] -= inv_jacobian / basis.weights[1] * ( + -0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + else + jump_φ = surface_values[1, 2, left] - φBC_left + div_q[1, element] += inv_jacobian / basis.weights[1] * 0.5 * inv_sqrtTᵣ * jump_φ + end + else + jump_φ = surface_values[1, 2, left] - surface_values[1, 1, left] + jump_q = surface_values[2, 2, left] - surface_values[2, 1, left] + div_q[1, element] -= inv_jacobian / basis.weights[1] * ( + -0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + end + + # right boundary correction + if element == size(φ, 2) + if φBC_right === :periodic + jump_φ = surface_values[1, 2, 1] - surface_values[1, 1, right] + jump_q = surface_values[2, 2, 1] - surface_values[2, 1, right] + div_q[end, element] += inv_jacobian / basis.weights[end] * ( + 0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + else + jump_φ = φBC_right - surface_values[1, 1, right] + div_q[end, element] -= inv_jacobian / basis.weights[end] * 0.5 * inv_sqrtTᵣ * jump_φ + end + else + jump_φ = surface_values[1, 2, right] - surface_values[1, 1, right] + jump_q = surface_values[2, 2, right] - surface_values[2, 1, right] + div_q[end, element] += inv_jacobian / basis.weights[end] * ( + 0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + end + end + + return nothing +end + + +# 2D functions +function compute_q!(q, φ::AbstractArray{<:Real, 4}, params) + (; basis, dl, dr, Tᵣ, + jacobian_x, jacobian_y, + surface_values_x, surface_values_y, + φBC_left, φBC_right, φBC_bottom, φBC_top) = params + sqrtTᵣ = sqrt(Tᵣ) + q1, q2 = q + + # prolong φ, Dφ to surfaces + for element_y in axes(φ, 4), element_x in axes(φ, 3) + left, right = element_x, element_x + 1 + bottom, top = element_y, element_y + 1 + inv_jacobian_x = inv(jacobian_x[element_x]) + inv_jacobian_y = inv(jacobian_y[element_y]) + + # x direction + for j in axes(φ, 2) + surface_values_x[1, j, 2, left, element_y] = φ[1, j, element_x, element_y] + surface_values_x[2, j, 2, left, element_y] = dot(dl, view(φ, :, j, element_x, element_y)) * inv_jacobian_x + surface_values_x[1, j, 1, right, element_y] = φ[end, j, element_x, element_y] + surface_values_x[2, j, 1, right, element_y] = dot(dr, view(φ, :, j, element_x, element_y)) * inv_jacobian_x + end + + # y direction + for i in axes(φ, 1) + surface_values_y[1, i, 2, bottom, element_x] = φ[i, 1, element_x, element_y] + surface_values_y[2, i, 2, bottom, element_x] = dot(dl, view(φ, i, :, element_x, element_y)) * inv_jacobian_y + surface_values_y[1, i, 1, top, element_x] = φ[i, end, element_x, element_y] + surface_values_y[2, i, 1, top, element_x] = dot(dr, view(φ, i, :, element_x, element_y)) * inv_jacobian_y + end + end + + # compute q1 locally + for element_y in axes(φ, 4), element_x in axes(φ, 3) + left, right = element_x, element_x + 1 + inv_jacobian = inv(jacobian_x[element_x]) + + for j in axes(φ, 2) + mul!(view(q1, :, j, element_x, element_y), basis.D, + view(φ, :, j, element_x, element_y), inv_jacobian, false) + + # left boundary correction + if element_x == 1 + if φBC_left === :periodic + jump_φ = surface_values_x[1, j, 2, left, element_y] - surface_values_x[1, j, 1, end, element_y] + jump_Dφ = surface_values_x[2, j, 2, left, element_y] - surface_values_x[2, j, 1, end, element_y] + inv_weight_r = inv(jacobian_x[element_x] * basis.weights[1]) + inv_weight_l = inv(jacobian_x[end ] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q1[1, j, element_x, element_y] += inv_jacobian / basis.weights[1] * ( + 0.5 * (1 + sqrtTᵣ * c1) * jump_φ - 0.5 * sqrtTᵣ * c2 * jump_Dφ) + else + jump_φ = surface_values_x[1, j, 2, left, element_y] - φBC_left[j, element_y] + q1[1, j, element_x, element_y] += inv_jacobian / basis.weights[1] * jump_φ + end + else + jump_φ = surface_values_x[1, j, 2, left, element_y] - surface_values_x[1, j, 1, left, element_y] + jump_Dφ = surface_values_x[2, j, 2, left, element_y] - surface_values_x[2, j, 1, left, element_y] + inv_weight_r = inv(jacobian_x[element_x ] * basis.weights[1]) + inv_weight_l = inv(jacobian_x[element_x-1] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q1[1, j, element_x, element_y] += inv_jacobian / basis.weights[1] * ( + 0.5 * (1 + sqrtTᵣ * c1) * jump_φ - 0.5 * sqrtTᵣ * c2 * jump_Dφ) + end + + # right boundary correction + if element_x == size(φ, 3) + if φBC_right === :periodic + jump_φ = surface_values_x[1, j, 2, 1, element_y] - surface_values_x[1, j, 1, right, element_y] + jump_Dφ = surface_values_x[2, j, 2, 1, element_y] - surface_values_x[2, j, 1, right, element_y] + inv_weight_r = inv(jacobian_x[1 ] * basis.weights[1]) + inv_weight_l = inv(jacobian_x[element_x] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q1[end, j, element_x, element_y] += inv_jacobian / basis.weights[end] * ( + 0.5 * (1 - sqrtTᵣ * c1) * jump_φ + 0.5 * sqrtTᵣ * c2 * jump_Dφ) + else + jump_φ = φBC_right[j, element_y] - surface_values_x[1, j, 1, right, element_y] + q1[end, j, element_x, element_y] += inv_jacobian / basis.weights[end] * jump_φ + end + else + jump_φ = surface_values_x[1, j, 2, right, element_y] - surface_values_x[1, j, 1, right, element_y] + jump_Dφ = surface_values_x[2, j, 2, right, element_y] - surface_values_x[2, j, 1, right, element_y] + inv_weight_r = inv(jacobian_x[element_x+1] * basis.weights[1]) + inv_weight_l = inv(jacobian_x[element_x ] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q1[end, j, element_x, element_y] += inv_jacobian / basis.weights[end] * ( + 0.5 * (1 - sqrtTᵣ * c1) * jump_φ + 0.5 * sqrtTᵣ * c2 * jump_Dφ) + end + end + end + + # compute q2 locally + for element_y in axes(φ, 4), element_x in axes(φ, 3) + bottom, top = element_y, element_y + 1 + inv_jacobian = inv(jacobian_y[element_y]) + + for i in axes(φ, 1) + mul!(view(q2, i, :, element_x, element_y), basis.D, + view(φ, i, :, element_x, element_y), inv_jacobian, false) + + # bottom boundary correction + if element_y == 1 + if φBC_bottom === :periodic + jump_φ = surface_values_y[1, i, 2, bottom, element_x] - surface_values_y[1, i, 1, end, element_x] + jump_Dφ = surface_values_y[2, i, 2, bottom, element_x] - surface_values_y[2, i, 1, end, element_x] + inv_weight_r = inv(jacobian_y[element_y] * basis.weights[1]) + inv_weight_l = inv(jacobian_y[end ] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q2[i, 1, element_x, element_y] += inv_jacobian / basis.weights[1] * ( + 0.5 * (1 + sqrtTᵣ * c1) * jump_φ - 0.5 * sqrtTᵣ * c2 * jump_Dφ) + else + jump_φ = surface_values_y[1, i, 2, bottom, element_x] - φBC_bottom[i, element_x] + q2[i, 1, element_x, element_y] += inv_jacobian / basis.weights[1] * jump_φ + end + else + jump_φ = surface_values_y[1, i, 2, bottom, element_x] - surface_values_y[1, i, 1, bottom, element_x] + jump_Dφ = surface_values_y[2, i, 2, bottom, element_x] - surface_values_y[2, i, 1, bottom, element_x] + inv_weight_r = inv(jacobian_y[element_y ] * basis.weights[1]) + inv_weight_l = inv(jacobian_y[element_y-1] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q2[i, 1, element_x, element_y] += inv_jacobian / basis.weights[1] * ( + 0.5 * (1 + sqrtTᵣ * c1) * jump_φ - 0.5 * sqrtTᵣ * c2 * jump_Dφ) + end + + # top boundary correction + if element_y == size(φ, 4) + if φBC_top === :periodic + jump_φ = surface_values_y[1, i, 2, 1, element_x] - surface_values_y[1, i, 1, top, element_x] + jump_Dφ = surface_values_y[2, i, 2, 1, element_x] - surface_values_y[2, i, 1, top, element_x] + inv_weight_r = inv(jacobian_y[1 ] * basis.weights[1]) + inv_weight_l = inv(jacobian_y[element_y] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q2[i, end, element_x, element_y] += inv_jacobian / basis.weights[end] * ( + 0.5 * (1 - sqrtTᵣ * c1) * jump_φ + 0.5 * sqrtTᵣ * c2 * jump_Dφ) + else + jump_φ = φBC_top[i, element_x] - surface_values_y[1, i, 1, top, element_x] + q2[i, end, element_x, element_y] += inv_jacobian / basis.weights[end] * jump_φ + end + else + jump_φ = surface_values_y[1, i, 2, top, element_x] - surface_values_y[1, i, 1, top, element_x] + jump_Dφ = surface_values_y[2, i, 2, top, element_x] - surface_values_y[2, i, 1, top, element_x] + inv_weight_r = inv(jacobian_y[element_y+1] * basis.weights[1]) + inv_weight_l = inv(jacobian_y[element_y ] * basis.weights[end]) + denominator = 1 + 0.5 * sqrtTᵣ * (inv_weight_r + inv_weight_l) + c1 = 0.5 * (inv_weight_r - inv_weight_l) / denominator + c2 = inv(denominator) + q2[i, end, element_x, element_y] += inv_jacobian / basis.weights[end] * ( + 0.5 * (1 - sqrtTᵣ * c1) * jump_φ + 0.5 * sqrtTᵣ * c2 * jump_Dφ) + end + end + end + + return nothing +end + +function compute_div_q!(div_q, φ::AbstractArray{<:Real, 4}, q, params) + (; basis, Tᵣ, + jacobian_x, jacobian_y, + surface_values_x, surface_values_y, + φBC_left, φBC_right, φBC_bottom, φBC_top) = params + sqrtTᵣ = sqrt(Tᵣ) + inv_sqrtTᵣ = inv(sqrtTᵣ) + q1, q2 = q + + # prolong φ, q to surfaces + for element_y in axes(φ, 4), element_x in axes(φ, 3) + left, right = element_x, element_x + 1 + bottom, top = element_y, element_y + 1 + + # x direction + for j in axes(φ, 2) + surface_values_x[1, j, 2, left , element_y] = φ[1, j, element_x, element_y] + surface_values_x[2, j, 2, left , element_y] = q1[1, j, element_x, element_y] + surface_values_x[1, j, 1, right, element_y] = φ[end, j, element_x, element_y] + surface_values_x[2, j, 1, right, element_y] = q1[end, j, element_x, element_y] + end + + # y direction + for i in axes(φ, 1) + surface_values_y[1, i, 2, bottom, element_x] = φ[i, 1, element_x, element_y] + surface_values_y[2, i, 2, bottom, element_x] = q2[i, 1, element_x, element_y] + surface_values_y[1, i, 1, top , element_x] = φ[i, end, element_x, element_y] + surface_values_y[2, i, 1, top , element_x] = q2[i, end, element_x, element_y] + end + end + + # compute div q locally + for element_y in axes(φ, 4), element_x in axes(φ, 3) + left, right = element_x, element_x + 1 + bottom, top = element_y, element_y + 1 + inv_jacobian_x = inv(jacobian_x[element_x]) + inv_jacobian_y = inv(jacobian_y[element_y]) + + # x direction + for j in axes(φ, 2) + mul!(view(div_q, :, j, element_x, element_y), basis.D, + view(q1, :, j, element_x, element_y), inv_jacobian_x, false) + + # left boundary correction + if element_x == 1 + if φBC_left === :periodic + jump_φ = surface_values_x[1, j, 2, left, element_y] - surface_values_x[1, j, 1, end, element_y] + jump_q = surface_values_x[2, j, 2, left, element_y] - surface_values_x[2, j, 1, end, element_y] + div_q[1, j, element_x, element_y] -= inv_jacobian_x / basis.weights[1] * ( + -0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + else + jump_φ = surface_values_x[1, j, 2, left, element_y] - φBC_left[j, element_y] + div_q[1, j, element_x, element_y] += inv_jacobian_x / basis.weights[1] * 0.5 * inv_sqrtTᵣ * jump_φ + end + else + jump_φ = surface_values_x[1, j, 2, left, element_y] - surface_values_x[1, j, 1, left, element_y] + jump_q = surface_values_x[2, j, 2, left, element_y] - surface_values_x[2, j, 1, left, element_y] + div_q[1, j, element_x, element_y] -= inv_jacobian_x / basis.weights[1] * ( + -0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + end + + # right boundary correction + if element_x == size(φ, 3) + if φBC_right === :periodic + jump_φ = surface_values_x[1, j, 2, 1, element_y] - surface_values_x[1, j, 1, right, element_y] + jump_q = surface_values_x[2, j, 2, 1, element_y] - surface_values_x[2, j, 1, right, element_y] + div_q[end, j, element_x, element_y] += inv_jacobian_x / basis.weights[end] * ( + 0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + else + jump_φ = φBC_right[j, element_y] - surface_values_x[1, j, 1, right, element_y] + div_q[end, j, element_x, element_y] -= inv_jacobian_x / basis.weights[end] * 0.5 * inv_sqrtTᵣ * jump_φ + end + else + jump_φ = surface_values_x[1, j, 2, right, element_y] - surface_values_x[1, j, 1, right, element_y] + jump_q = surface_values_x[2, j, 2, right, element_y] - surface_values_x[2, j, 1, right, element_y] + div_q[end, j, element_x, element_y] += inv_jacobian_x / basis.weights[end] * ( + 0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + end + end + + # y direction + for i in axes(φ, 1) + mul!(view(div_q, i, :, element_x, element_y), basis.D, + view(q2, i, :, element_x, element_y), inv_jacobian_y, true) + + # bottom boundary correction + if element_y == 1 + if φBC_bottom === :periodic + jump_φ = surface_values_y[1, i, 2, bottom, element_x] - surface_values_y[1, i, 1, end, element_x] + jump_q = surface_values_y[2, i, 2, bottom, element_x] - surface_values_y[2, i, 1, end, element_x] + div_q[i, 1, element_x, element_y] -= inv_jacobian_y / basis.weights[1] * ( + -0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + else + jump_φ = surface_values_y[1, i, 2, bottom, element_x] - φBC_bottom[i, element_x] + div_q[i, 1, element_x, element_y] += inv_jacobian_y / basis.weights[1] * 0.5 * inv_sqrtTᵣ * jump_φ + end + else + jump_φ = surface_values_y[1, i, 2, bottom, element_x] - surface_values_y[1, i, 1, bottom, element_x] + jump_q = surface_values_y[2, i, 2, bottom, element_x] - surface_values_y[2, i, 1, bottom, element_x] + div_q[i, 1, element_x, element_y] -= inv_jacobian_y / basis.weights[1] * ( + -0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + end + + # top boundary correction + if element_y == size(φ, 4) + if φBC_top === :periodic + jump_φ = surface_values_y[1, i, 2, 1, element_x] - surface_values_y[1, i, 1, top, element_x] + jump_q = surface_values_y[2, i, 2, 1, element_x] - surface_values_y[2, i, 1, top, element_x] + div_q[i, end, element_x, element_y] += inv_jacobian_y / basis.weights[end] * ( + 0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + else + jump_φ = φBC_top[i, element_x] - surface_values_y[1, i, 1, top, element_x] + div_q[i, end, element_x, element_y] -= inv_jacobian_y / basis.weights[end] * 0.5 * inv_sqrtTᵣ * jump_φ + end + else + jump_φ = surface_values_y[1, i, 2, top, element_x] - surface_values_y[1, i, 1, top, element_x] + jump_q = surface_values_y[2, i, 2, top, element_x] - surface_values_y[2, i, 1, top, element_x] + div_q[i, end, element_x, element_y] += inv_jacobian_y / basis.weights[end] * ( + 0.5 * jump_q + 0.5 * inv_sqrtTᵣ * jump_φ) + end + end + end + + return nothing +end + + +function hypdiff_parabolic!(dφ, φ, params, t) + (; f, q) = params + + compute_q!(q, φ, params) + + compute_div_q!(dφ, φ, q, params) + + @. dφ = dφ + f + + return nothing +end + +function hypdiff_parabolic_for_cg!(dφ_vector, φ_vector, params) + (; q) = params + if q isa Array + # 1D + s = size(q) + else + # 2D + s = size(q[1]) + end + dφ = reshape(dφ_vector, s) + φ = reshape(φ_vector, s) + + # need to set Dirichlet boundary conditions to homogeneous since the + # non-homogeneous BCs are contained in the right-hand side vector + params_homogeneous = params + if haskey(params, :φBC_left) + if params.φBC_left isa Number + params_homogeneous = (; params_homogeneous..., φBC_left = zero(params.φBC_left)) + elseif params.φBC_left isa AbstractArray + params_homogeneous = (; params_homogeneous..., φBC_left = Zeros(size(params.φBC_left))) + end + end + if haskey(params, :φBC_right) + if params.φBC_right isa Number + params_homogeneous = (; params_homogeneous..., φBC_right = zero(params.φBC_right)) + elseif params.φBC_right isa AbstractArray + params_homogeneous = (; params_homogeneous..., φBC_right = Zeros(size(params.φBC_right))) + end + end + if haskey(params, :φBC_bottom) + if params.φBC_bottom isa Number + params_homogeneous = (; params_homogeneous..., φBC_bottom = zero(params.φBC_bottom)) + elseif params.φBC_bottom isa AbstractArray + params_homogeneous = (; params_homogeneous..., φBC_bottom = Zeros(size(params.φBC_bottom))) + end + end + if haskey(params, :φBC_top) + if params.φBC_top isa Number + params_homogeneous = (; params_homogeneous..., φBC_top = zero(params.φBC_top)) + elseif params.φBC_top isa AbstractArray + params_homogeneous = (; params_homogeneous..., φBC_top = Zeros(size(params.φBC_top))) + end + end + + compute_q!(q, φ, params_homogeneous) + + compute_div_q!(dφ, φ, q, params_homogeneous) + + mul_by_mass_matrix!(dφ, params_homogeneous) + + return nothing +end + + +# 1D function +function mul_by_mass_matrix!(u::AbstractArray{<:Real, 2}, params) + (; basis, jacobian) = params + weights = basis.weights + + for element in axes(u, 2) + factor = jacobian[element] + for i in axes(u, 1) + u[i, element] *= factor * weights[i] + end + end + + return nothing +end + +# 2D function +function mul_by_mass_matrix!(u::AbstractArray{<:Real, 4}, params) + (; basis, jacobian_x, jacobian_y) = params + weights = basis.weights + + for element_y in axes(u, 4), element_x in axes(u, 3) + factor = jacobian_x[element_x] * jacobian_y[element_y] + for j in axes(u, 2), i in axes(u, 1) + u[i, j, element_x, element_y] *= factor * weights[i] * weights[j] + end + end + + return nothing +end + + +# 1D functions +function l2_norm(u1, u2::AbstractArray{<:Real, 2}, params) + (; basis, jacobian) = params + (; weights) = basis + + res = zero((u1[1] - u2[1]) * weights[1]) + for element in axes(u1, 2) + factor = jacobian[element] + for i in axes(u1, 1) + res += factor * weights[i] * (u1[i, element] - u2[i, element])^2 + end + end + + return sqrt(res) +end + +function integral(u::AbstractArray{<:Real, 2}, params) + (; basis, jacobian) = params + (; weights) = basis + + res = zero(u[1] * weights[1]) + for element in axes(u, 2) + factor = jacobian[element] + for i in axes(u, 1) + res += factor * weights[i] * u[i, element] + end + end + + return res +end + + +# 2D functions +function l2_norm(u1, u2::AbstractArray{<:Real, 4}, params) + (; basis, jacobian_x, jacobian_y) = params + (; weights) = basis + + res = zero((u1[1] - u2[1]) * weights[1]) + for element_y in axes(u1, 4), element_x in axes(u1, 3) + factor = jacobian_x[element_x] * jacobian_y[element_y] + for j in axes(u1, 2), i in axes(u1, 1) + res += factor * weights[i] * weights[j] * + (u1[i, j, element_x, element_y] - u2[i, j, element_x, element_y])^2 + end + end + + return sqrt(res) +end + +function integral(u::AbstractArray{<:Real, 4}, params) + (; basis, jacobian_x, jacobian_y) = params + (; weights) = basis + + res = zero(u[1] * weights[1]) + for element_y in axes(u, 4), element_x in axes(u, 3) + factor = jacobian_x[element_x] * jacobian_y[element_y] + for j in axes(u, 2), i in axes(u, 1) + res += factor * weights[i] * weights[j] * u[i, j, element_x, element_y] + end + end + + return res +end + + +# utility function +function compute_eoc(Ns, errors) + eoc = similar(errors) + eoc[begin] = NaN # no EOC defined for the first grid + for idx in Iterators.drop(eachindex(errors, Ns, eoc), 1) + eoc[idx] = -( log(errors[idx] / errors[idx - 1]) / log(Ns[idx] / Ns[idx - 1]) ) + end + return eoc +end + + +# 1D functions +function setup_1d(xmin, xmax, polydeg, N) + Lᵣ = (xmax - xmin) / (2π) + Tᵣ = Lᵣ^2 # std. choice of Nishikawa, this makes all modes underdamped + # Tᵣ = Lᵣ^2 / 4 # this makes the lowest mode damped ideally + + basis = LobattoLegendre(polydeg) + element_boundaries = range(xmin, stop=xmax, length=N+1) # uniform mesh + + x = zeros(length(basis.weights), length(element_boundaries) - 1) + jacobian = zeros(length(element_boundaries) - 1) + for element in axes(x, 2) + jac = (element_boundaries[element+1] - element_boundaries[element]) / + (basis.nodes[end] - basis.nodes[1]) + @. x[:, element] = element_boundaries[element] + (basis.nodes - basis.nodes[1]) * jac + jacobian[element] = jac + end + + surface_values = similar(x, 2, 2, length(element_boundaries)) # (φ,q), (left,right), surface + + dl = basis.D[1, :] + dr = basis.D[end, :] + + return (; basis, dl, dr, jacobian, x, Tᵣ, surface_values) +end + +function solve_1d_nonperiodic(; polydeg = 3, N = 10) + xmin, xmax = -1.0, 1.0 + + φ_func(x) = exp(-10 * x^2) + q_func(x) = -20 * x * exp(-10 * x^2) + f_func(x) = -20 * (20 * x^2 - 1) * exp(-10 * x^2) + + (; basis, dl, dr, jacobian, x, Tᵣ, surface_values) = setup_1d( + xmin, xmax, polydeg, N) + + φ0 = zero(x) + q = similar(φ0) + f = f_func.(x) + tmp = similar(φ0) + φsol = φ_func.(x) + qsol = q_func.(x) + φBC_left = φ_func(first(x)) + φBC_right = φ_func(last(x)) + params = (; basis, dl, dr, jacobian, x, Tᵣ, surface_values, f, q, φBC_left, φBC_right, tmp, φsol, qsol) + + A, b = compute_sparse_linear_structure(hypdiff_parabolic!, copy(φ0), params, 0.0) + φ = reshape(A \ b, size(φ0)) + q = similar(φ) + compute_q!(q, φ, params) + + error_φ = l2_norm(φ, φsol, params) + error_q = l2_norm(q, qsol, params) + return (; error_φ, error_q, x, φ, q, φsol, qsol, params) +end + +function solve_1d_periodic(; polydeg = 3, N = 10) + xmin, xmax = -2.0, 2.0 + + # needs to have a vanishing mean value for periodic BCs - assume the domain + # is symmetric around the origin + φ_func(x) = exp(-10 * x^2) - sqrt(pi / 10) * erf(sqrt(10) * xmax) / (xmax - xmin) + q_func(x) = -20 * x * exp(-10 * x^2) + f_func(x) = -20 * (20 * x^2 - 1) * exp(-10 * x^2) + + (; basis, dl, dr, jacobian, x, Tᵣ, surface_values) = setup_1d( + xmin, xmax, polydeg, N) + + φ0 = zero(x) + q = similar(φ0) + f = f_func.(x) + tmp = similar(φ0) + φsol = φ_func.(x) + qsol = q_func.(x) + φBC_left = :periodic + φBC_right = :periodic + params = (; basis, dl, dr, jacobian, x, Tᵣ, surface_values, f, q, φBC_left, φBC_right, tmp, φsol, qsol) + + A, b = compute_sparse_linear_structure(hypdiff_parabolic!, copy(φ0), params, 0.0) + φ = reshape(A \ b, size(φ0)) + # subtract mean value since direct solvers may not give the solution with + # vanishing mean value but, e.g., the least norm solution + φ .= φ .- integral(φ, params) / (xmax - xmin) + q = similar(φ) + compute_q!(q, φ, params) + + error_φ = l2_norm(φ, φsol, params) + error_q = l2_norm(q, qsol, params) + return (; error_φ, error_q, x, φ, q, φsol, qsol, params) +end + + +@inline function Trixi.flux_godunov(u_ll, u_rr, orientation::Integer, equations::HyperbolicDiffusionEquations1D) + # Obtain left and right fluxes + phi_ll, q1_ll = u_ll + phi_rr, q1_rr = u_rr + f_ll = flux(u_ll, orientation, equations) + f_rr = flux(u_rr, orientation, equations) + + # this is an optimized version of the application of the upwind dissipation matrix: + # dissipation = 0.5*R_n*|Λ|*inv(R_n)[[u]] + λ_max = sqrt(equations.nu * equations.inv_Tr) + f1 = 1/2 * (f_ll[1] + f_rr[1]) - 1/2 * λ_max * (phi_rr - phi_ll) + f2 = 1/2 * (f_ll[2] + f_rr[2]) - 1/2 * λ_max * (q1_rr - q1_ll) + + return SVector(f1, f2) +end + +function solve_1d_nonperiodic_trixi(; polydeg = 3, N = 10, + surface_flux = flux_godunov) + xmin, xmax = -1.0, 1.0 + + φ_func(x) = exp(-10 * x^2) + q_func(x) = -20 * x * exp(-10 * x^2) + f_func(x) = -20 * (20 * x^2 - 1) * exp(-10 * x^2) + + initial_condition(x, t, equations) = if iszero(t) + return SVector(zero(t), zero(t)) + else + return SVector(φ_func(x[1]), q_func(x[1])) + end + source_terms(u, x, t, equations) = let + harmonic = source_terms_harmonic(u, x, t, equations) + return SVector(f_func(x[1]), 0) + harmonic + end + boundary_condition_nonperiodic(u_inner, orientation, direction, x, t, + surface_flux_function, equations) = let + u_boundary = initial_condition(x, one(t), equations) + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" + flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + end + return flux + end + + equations = HyperbolicDiffusionEquations1D(Lr = (xmax - xmin) / (2π)) + solver = DGSEM(; polydeg, surface_flux) + mesh = StructuredMesh((N,), (xmin,), (xmax,), periodicity=(false,)) + semi = SemidiscretizationHyperbolic( + mesh, equations, initial_condition, solver; + source_terms, boundary_conditions = boundary_condition_nonperiodic) + + ode = semidiscretize(semi, (0.0, 100.0)) + steady_state_callback = SteadyStateCallback(abstol = 5.0e-12, reltol = 0.0) + stepsize_callback = StepsizeCallback(cfl = 1.5) + callbacks = CallbackSet(steady_state_callback, stepsize_callback) + sol = redirect_stderr(devnull) do + # avoid cluttering the REPL with output of the steady_state_callback + Trixi.solve( + ode, Trixi.HypDiffN3Erk3Sstar52(), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks) + end + + usol = Trixi.compute_coefficients(initial_condition, 1.0, semi) + diff = sol.u[end] - usol + error_φ = integrate(abs2 ∘ ((u, equations) -> u[1]), diff, semi; normalize = false) |> sqrt + error_q = integrate(abs2 ∘ ((u, equations) -> u[2]), diff, semi; normalize = false) |> sqrt + return (; error_φ, error_q, sol, usol, semi) +end + +function solve_1d_periodic_trixi(; polydeg = 3, N = 10, + surface_flux = flux_godunov) + xmin, xmax = -2.0, 2.0 + + # needs to have a vanishing mean value for periodic BCs - assume the domain + # is symmetric around the origin + φ_func(x) = exp(-10 * x^2) - sqrt(pi / 10) * erf(sqrt(10) * xmax) / (xmax - xmin) + q_func(x) = -20 * x * exp(-10 * x^2) + f_func_tmp(x) = -20 * (20 * x^2 - 1) * exp(-10 * x^2) + + initial_condition(x, t, equations) = if iszero(t) + return SVector(zero(t), zero(t)) + else + return SVector(φ_func(x[1]), q_func(x[1])) + end + + equations = HyperbolicDiffusionEquations1D(Lr = (xmax - xmin) / (2π)) + # equations = HyperbolicDiffusionEquations1D(Lr = 1.0) + solver = DGSEM(; polydeg, surface_flux) + mesh = StructuredMesh((N,), (xmin,), (xmax,), periodicity=(true,)) + semi = SemidiscretizationHyperbolic( + mesh, equations, initial_condition, solver) + + # the source term must have discretely vanishing mean value for periodic BCs + source = Trixi.compute_coefficients(1.0, semi) do x, t, equations + SVector(f_func_tmp(x[1]), 0) + end + mean_value = integrate(first, source, semi) + f_func(x) = let mean_value = mean_value + return -20 * (20 * x^2 - 1) * exp(-10 * x^2) - mean_value + end + + source_terms(u, x, t, equations) = let + harmonic = source_terms_harmonic(u, x, t, equations) + return SVector(f_func(x[1]), 0) + harmonic + end + + semi = SemidiscretizationHyperbolic( + mesh, equations, initial_condition, solver; + source_terms, boundary_conditions = boundary_condition_periodic) + + ode = semidiscretize(semi, (0.0, 1.0e3)) + steady_state_callback = SteadyStateCallback(abstol = 5.0e-12, reltol = 0.0) + stepsize_callback = StepsizeCallback(cfl = 1.5) + callbacks = CallbackSet(steady_state_callback, stepsize_callback) + sol = redirect_stderr(devnull) do + # avoid cluttering the REPL with output of the steady_state_callback + Trixi.solve( + ode, Trixi.HypDiffN3Erk3Sstar52(), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks) + end + + usol = Trixi.compute_coefficients(initial_condition, 1.0, semi) + diff = sol.u[end] - usol + error_φ = integrate(abs2 ∘ ((u, equations) -> u[1]), diff, semi; normalize = false) |> sqrt + error_q = integrate(abs2 ∘ ((u, equations) -> u[2]), diff, semi; normalize = false) |> sqrt + return (; error_φ, error_q, sol, usol, semi, source) +end + + +# 2D functions +function setup_2d(xmin, xmax, ymin, ymax, polydeg, Nx, Ny) + # equation (26) of Nishikawa and Nakashima (2018) + Lᵣ = (xmax - xmin) * (ymax - ymin) / (2π * sqrt((xmax - xmin)^2 + (ymax - ymin)^2)) + Tᵣ = Lᵣ^2 # std. choice of Nishikawa, this makes all modes underdamped + # Tᵣ = Lᵣ^2 / 4 # this makes the lowest mode damped ideally + + basis = LobattoLegendre(polydeg) + element_boundaries_x = range(xmin, stop = xmax, length = Nx + 1) # uniform mesh + element_boundaries_y = range(ymin, stop = ymax, length = Ny + 1) # uniform mesh + + x = zeros(length(basis.weights), length(basis.weights), + length(element_boundaries_x) - 1, length(element_boundaries_y) - 1) + y = zeros(length(basis.weights), length(basis.weights), + length(element_boundaries_x) - 1, length(element_boundaries_y) - 1) + jacobian_x = zeros(length(element_boundaries_x) - 1) + jacobian_y = zeros(length(element_boundaries_y) - 1) + for element_x in axes(x, 3) + jac = (element_boundaries_x[element_x+1] - element_boundaries_x[element_x]) / + (basis.nodes[end] - basis.nodes[1]) + jacobian_x[element_x] = jac + for element_y in axes(x, 4), j in axes(x, 2) + @. x[:, j, element_x, element_y] = element_boundaries_x[element_x] + (basis.nodes - basis.nodes[1]) * jac + end + end + for element_y in axes(y, 4) + jac = (element_boundaries_y[element_y+1] - element_boundaries_y[element_y]) / + (basis.nodes[end] - basis.nodes[1]) + jacobian_y[element_y] = jac + for element_x in axes(y, 3), i in axes(y, 1) + @. y[i, :, element_x, element_y] = element_boundaries_y[element_y] + (basis.nodes - basis.nodes[1]) * jac + end + end + + # (φ,q), other_nodes, (left,right), surface, eleemnt_in_other_direction + surface_values_x = similar(x, 2, size(x, 2), 2, length(element_boundaries_x), size(x, 4)) + surface_values_y = similar(x, 2, size(x, 1), 2, length(element_boundaries_y), size(x, 3)) + + dl = basis.D[1, :] + dr = basis.D[end, :] + + return (; basis, dl, dr, Tᵣ, + jacobian_x, jacobian_y, + x, y, + surface_values_x, surface_values_y) +end + +function solve_2d_nonperiodic(; polydeg = 3, N = 10, Nx = N, Ny = N) + # setup of Cockburn, Guzman, Wang (2008) + xmin, xmax = -0.5, 0.5 + ymin, ymax = -0.5, 0.5 + φ_func(x, y) = cospi(x) * cospi(y) + q1_func(x, y) = -π * sinpi(x) * cospi(y) + q2_func(x, y) = -π * cospi(x) * sinpi(y) + f_func(x, y) = 2 * π^2 * cospi(x) * cospi(y) + + (; basis, dl, dr, Tᵣ, jacobian_x, jacobian_y, x, y, + surface_values_x, surface_values_y) = setup_2d(xmin, xmax, ymin, ymax, + polydeg, Nx, Ny) + + φ0 = zero(x) + q1 = similar(φ0) + q2 = similar(φ0) + q = (q1, q2) + f = f_func.(x, y) + tmp = similar(φ0) + φsol = φ_func.(x, y) + q1sol = q1_func.(x, y) + q2sol = q2_func.(x, y) + φBC_left = φ_func.(x[1, :, 1, :], y[1, :, 1, :]) + φBC_right = φ_func.(x[end, :, end, :], y[end, :, end, :]) + φBC_bottom = φ_func.(x[:, 1, :, 1, ], y[:, 1, :, 1, ]) + φBC_top = φ_func.(x[:, end, :, end], y[:, end, :, end]) + params = (; basis, dl, dr, Tᵣ, jacobian_x, jacobian_y, x, y, + surface_values_x, surface_values_y, + f, q, φBC_left, φBC_right, φBC_bottom, φBC_top, tmp, + φsol, q1sol, q2sol) + + # # using a sparse direct solver is slow if we compute the sparse matrix + # # via many matrix-vector products + # A, b = compute_sparse_linear_structure(hypdiff_parabolic!, copy(φ0), params, 0.0) + # φ = reshape(A \ b, size(φ0)) + + # using CG is more efficient if we have only cheap matrix-vector products + b = zeros(length(φ0)) + hypdiff_parabolic!(reshape(b, size(φ0)), φ0, params, 0.0) + @. b = -b + mul_by_mass_matrix!(reshape(b, size(φ0)), params) + A = LinearOperator(eltype(b), length(b), length(b), true, true, + (dφ, φ) -> hypdiff_parabolic_for_cg!(dφ, φ, params)) + # use slightly stricter relative tolerance than the default + # sqrt(eps(eltype(b)) to get full accuracy for higher grid resolutions + φ_vector, stats = Krylov.cg(A, b; rtol = 1.0e-12) + # @show stats + φ = reshape(φ_vector, size(φ0)) + + q1 = similar(φ) + q2 = similar(φ) + q = (q1, q2) + compute_q!(q, φ, params) + + error_φ = l2_norm(φ, φsol, params) + error_q1 = l2_norm(q1, q1sol, params) + error_q2 = l2_norm(q2, q2sol, params) + return (; error_φ, error_q1, error_q2, + x, y, φ, q1, q2, φsol, q1sol, q2sol, params) +end + +function solve_2d_periodic(; polydeg = 3, N = 10, Nx = N, Ny = N) + xmin, xmax = -1.0, 1.0 + ymin, ymax = -1.0, 1.0 + + φ_func(x, y) = 2 * cospi(x) * sinpi(2 * y) + q1_func(x, y) = -π * 2 * sinpi(x) * sinpi(2 * y) + q2_func(x, y) = π * 4 * cospi(x) * cospi(2 * y) + f_func(x, y) = π^2*10 * cospi(x) * sinpi(2 * y) + + (; basis, dl, dr, Tᵣ, jacobian_x, jacobian_y, x, y, + surface_values_x, surface_values_y) = setup_2d(xmin, xmax, ymin, ymax, + polydeg, Nx, Ny) + + φ0 = zero(x) + q1 = similar(φ0) + q2 = similar(φ0) + q = (q1, q2) + f = f_func.(x, y) + tmp = similar(φ0) + φsol = φ_func.(x, y) + q1sol = q1_func.(x, y) + q2sol = q2_func.(x, y) + φBC_left = :periodic + φBC_right = :periodic + φBC_bottom = :periodic + φBC_top = :periodic + params = (; basis, dl, dr, Tᵣ, jacobian_x, jacobian_y, x, y, + surface_values_x, surface_values_y, + f, q, φBC_left, φBC_right, φBC_bottom, φBC_top, tmp, + φsol, q1sol, q2sol) + + b = zeros(length(φ0)) + hypdiff_parabolic!(reshape(b, size(φ0)), φ0, params, 0.0) + @. b = -b + mul_by_mass_matrix!(reshape(b, size(φ0)), params) + A = LinearOperator(eltype(b), length(b), length(b), true, true, + (dφ, φ) -> hypdiff_parabolic_for_cg!(dφ, φ, params)) + # use slightly stricter relative tolerance than the default + # sqrt(eps(eltype(b)) to get full accuracy for higher grid resolutions + φ_vector, stats = Krylov.cg(A, b; rtol = 1.0e-12) + # @show stats + φ = reshape(φ_vector, size(φ0)) + + q1 = similar(φ) + q2 = similar(φ) + q = (q1, q2) + compute_q!(q, φ, params) + + error_φ = l2_norm(φ, φsol, params) + error_q1 = l2_norm(q1, q1sol, params) + error_q2 = l2_norm(q2, q2sol, params) + return (; error_φ, error_q1, error_q2, + x, y, φ, q1, q2, φsol, q1sol, q2sol, params) +end + + +@inline function dirichlet_boundary_flux(u_inner, u_boundary, orientation, direction, equations::HyperbolicDiffusionEquations2D) + # Calculate boundary flux by setting the value of the incoming characteristic + # variables to the outgoing characteristic variables plus the boundary value. + # This version is made to impose a BC on phi only. + phi_bc, q1_bc, q2_bc = u_boundary + phi, q1, q2 = u_inner + inv_Tr = equations.inv_Tr + sqrt_inv_Tr = sqrt(inv_Tr) + # Curved version using the `orientation`, which is a vector + if direction == 1 # -x + flux = SVector(-q1 - sqrt_inv_Tr * (phi - phi_bc), + -inv_Tr * phi_bc, + zero(phi_bc)) * orientation[1] + elseif direction == 2 # +x + flux = SVector(-q1 + sqrt_inv_Tr * (phi - phi_bc), + -inv_Tr * phi_bc, + zero(phi_bc)) * orientation[1] + elseif direction == 3 # -y + flux = SVector(-q2 - sqrt_inv_Tr * (phi - phi_bc), + zero(phi_bc), + -inv_Tr * phi_bc) * orientation[2] + else #if direction == 4 # +y + flux = SVector(-q2 + sqrt_inv_Tr * (phi - phi_bc), + zero(phi_bc), + -inv_Tr * phi_bc) * orientation[2] + end + + return flux +end + +function solve_2d_nonperiodic_trixi(; polydeg = 3, N = 10, Nx = N, Ny = N, + surface_flux = flux_godunov) + # setup of Cockburn, Guzman, Wang (2008) + xmin, xmax = -0.5, 0.5 + ymin, ymax = -0.5, 0.5 + φ_func(x, y) = cospi(x) * cospi(y) + q1_func(x, y) = -π * sinpi(x) * cospi(y) + q2_func(x, y) = -π * cospi(x) * sinpi(y) + f_func(x, y) = 2 * π^2 * cospi(x) * cospi(y) + + initial_condition(x, t, equations) = if iszero(t) + return SVector(zero(t), zero(t), zero(t)) + else + return SVector(φ_func(x[1], x[2]), q1_func(x[1], x[2]), q2_func(x[1], x[2])) + end + source_terms(u, x, t, equations) = let + harmonic = source_terms_harmonic(u, x, t, equations) + return SVector(f_func(x[1], x[2]), 0, 0) + harmonic + end + boundary_condition_nonperiodic(u_inner, orientation, direction, x, t, + surface_flux_function, equations) = let + u_boundary = initial_condition(x, one(t), equations) + flux = dirichlet_boundary_flux(u_inner, u_boundary, orientation, direction, equations) + return flux + end + + equations = HyperbolicDiffusionEquations2D(Lr = + (xmax - xmin) * (ymax - ymin) / (2π * sqrt((xmax - xmin)^2 + (ymax - ymin)^2))) + solver = DGSEM(; polydeg, surface_flux) + mesh = StructuredMesh((Nx, Ny), (xmin, ymin), (xmax, ymax), periodicity=(false, false)) + semi = SemidiscretizationHyperbolic( + mesh, equations, initial_condition, solver; + source_terms, boundary_conditions = ( + x_neg = boundary_condition_nonperiodic, + x_pos = boundary_condition_nonperiodic, + y_neg = boundary_condition_nonperiodic, + y_pos = boundary_condition_nonperiodic,)) + + ode = semidiscretize(semi, (0.0, 1.0e2)) + steady_state_callback = SteadyStateCallback(abstol = 5.0e-12, reltol = 0.0) + stepsize_callback = StepsizeCallback(cfl = 1.5) + callbacks = CallbackSet(steady_state_callback, stepsize_callback) + sol = redirect_stderr(devnull) do + # avoid cluttering the REPL with output of the steady_state_callback + Trixi.solve( + ode, Trixi.HypDiffN3Erk3Sstar52(), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks) + end + + usol = Trixi.compute_coefficients(initial_condition, 1.0, semi) + diff = sol.u[end] - usol + error_φ = integrate(abs2 ∘ ((u, equations) -> u[1]), diff, semi; normalize = false) |> sqrt + error_q1 = integrate(abs2 ∘ ((u, equations) -> u[2]), diff, semi; normalize = false) |> sqrt + error_q2 = integrate(abs2 ∘ ((u, equations) -> u[3]), diff, semi; normalize = false) |> sqrt + return (; error_φ, error_q1, error_q2, sol, usol, semi) +end + +function solve_2d_periodic_trixi(; polydeg = 3, N = 10, Nx = N, Ny = N, + surface_flux = flux_godunov) + xmin, xmax = -1.0, 1.0 + ymin, ymax = -1.0, 1.0 + + φ_func(x, y) = 2 * cospi(x) * sinpi(2 * y) + q1_func(x, y) = -π * 2 * sinpi(x) * sinpi(2 * y) + q2_func(x, y) = π * 4 * cospi(x) * cospi(2 * y) + f_func(x, y) = π^2*10 * cospi(x) * sinpi(2 * y) + + initial_condition(x, t, equations) = if iszero(t) + return SVector(zero(t), zero(t), zero(t)) + else + return SVector(φ_func(x[1], x[2]), q1_func(x[1], x[2]), q2_func(x[1], x[2])) + end + source_terms(u, x, t, equations) = let + harmonic = source_terms_harmonic(u, x, t, equations) + return SVector(f_func(x[1], x[2]), 0, 0) + harmonic + end + + equations = HyperbolicDiffusionEquations2D(Lr = + (xmax - xmin) * (ymax - ymin) / (2π * sqrt((xmax - xmin)^2 + (ymax - ymin)^2))) + solver = DGSEM(; polydeg, surface_flux) + mesh = StructuredMesh((Nx, Ny), (xmin, ymin), (xmax, ymax), periodicity=(true, true)) + semi = SemidiscretizationHyperbolic( + mesh, equations, initial_condition, solver; + source_terms, boundary_conditions = boundary_condition_periodic) + + ode = semidiscretize(semi, (0.0, 1.0e2)) + steady_state_callback = SteadyStateCallback(abstol = 5.0e-12, reltol = 0.0) + stepsize_callback = StepsizeCallback(cfl = 1.5) + callbacks = CallbackSet(steady_state_callback, stepsize_callback) + sol = redirect_stderr(devnull) do + # avoid cluttering the REPL with output of the steady_state_callback + Trixi.solve( + ode, Trixi.HypDiffN3Erk3Sstar52(), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks) + end + + usol = Trixi.compute_coefficients(initial_condition, 1.0, semi) + diff = sol.u[end] - usol + error_φ = integrate(abs2 ∘ ((u, equations) -> u[1]), diff, semi; normalize = false) |> sqrt + error_q1 = integrate(abs2 ∘ ((u, equations) -> u[2]), diff, semi; normalize = false) |> sqrt + error_q2 = integrate(abs2 ∘ ((u, equations) -> u[3]), diff, semi; normalize = false) |> sqrt + return (; error_φ, error_q1, error_q2, sol, usol, semi) +end + + +# 1D functions +function convergence_tests_1d(solve, polydeg, Ns; latex = false) + errors_φ = Vector{Float64}() + errors_q = Vector{Float64}() + + for N in Ns + res = solve(; polydeg, N) + push!(errors_φ, res.error_φ) + push!(errors_q, res.error_q) + end + eoc_φ = compute_eoc(Ns, errors_φ) + eoc_q = compute_eoc(Ns, errors_q) + + # print results + data = hcat(Ns, errors_φ, eoc_φ, errors_q, eoc_q) + header = ["N", "L2 error φ", "L2 EOC φ", "L2 error q", "L2 EOC q"] + kwargs = (; header, formatters=(ft_printf("%3d", [1]), + ft_printf("%.2e", [2, 4]), + ft_printf("%.2f", [3, 5]))) + pretty_table(data; kwargs...) + if latex + pretty_table(data; kwargs..., backend=Val(:latex)) + end + + return nothing +end + +function convergence_tests_1d(; kwargs...) + Ns = 10 .* 2 .^ (0:6) + + for polydeg in 2:3 + println() + @info "Convergence tests 1D nonperiodic" polydeg + convergence_tests_1d(solve_1d_nonperiodic, polydeg, Ns; kwargs...) + + println() + @info "Convergence tests 1D nonperiodic with Trixi.jl" polydeg + convergence_tests_1d(solve_1d_nonperiodic_trixi, polydeg, Ns[1:end-2]; kwargs...) + end + + for polydeg in 2:3 + println() + @info "Convergence tests 1D periodic" polydeg + convergence_tests_1d(solve_1d_periodic, polydeg, Ns; kwargs...) + + println() + @info "Convergence tests 1D periodic with Trixi.jl" polydeg + convergence_tests_1d(solve_1d_periodic_trixi, polydeg, Ns[1:end-2]; kwargs...) + end + + return nothing +end + + +# 2D functions +function convergence_tests_2d(solve, polydeg, Ns; latex = false) + errors_φ = Vector{Float64}() + errors_q1 = Vector{Float64}() + errors_q2 = Vector{Float64}() + + for N in Ns + res = solve(; polydeg, N) + push!(errors_φ, res.error_φ) + push!(errors_q1, res.error_q1) + push!(errors_q2, res.error_q2) + end + eoc_φ = compute_eoc(Ns, errors_φ) + eoc_q1 = compute_eoc(Ns, errors_q1) + eoc_q2 = compute_eoc(Ns, errors_q2) + + # print results + data = hcat(Ns, errors_φ, eoc_φ, errors_q1, eoc_q1, errors_q2, eoc_q2) + header = ["N", "L2 error φ", "L2 EOC φ", "L2 error q1", "L2 EOC q1", + "L2 error q2", "L2 EOC q2"] + kwargs = (; header, formatters=(ft_printf("%3d", [1]), + ft_printf("%.2e", [2, 4, 6]), + ft_printf("%.2f", [3, 5, 7]))) + pretty_table(data; kwargs...) + if latex + pretty_table(data; kwargs..., backend=Val(:latex)) + end + + return nothing +end + +function convergence_tests_2d(; kwargs...) + Ns = 4 .* 2 .^ (0:4) + + for polydeg in 2:3 + println() + @info "Convergence tests 2D nonperiodic" polydeg + convergence_tests_2d(solve_2d_nonperiodic, polydeg, Ns; kwargs...) + end + + for polydeg in 2:3 + println() + @info "Convergence tests 2D periodic" polydeg + convergence_tests_2d(solve_2d_periodic, polydeg, Ns; kwargs...) + end + + return nothing +end