diff --git a/NEWS.md b/NEWS.md index 3a3a504a911..02a723fca45 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,7 @@ for human readability. - `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D` - Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, can now be digested by Trixi in 2D and 3D. +- Subcell (positivity) limiting support for nonlinear variables in 2D for `TreeMesh` ## Changes when updating to v0.6 from v0.5.x diff --git a/Project.toml b/Project.toml index 0bbdec206d8..e99b08e0e81 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -57,6 +58,7 @@ ConstructionBase = "1.3" DataStructures = "0.18.15" DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" +Downloads = "1.6" EllipsisNotation = "1.0" FillArrays = "0.13.2, 1" ForwardDiff = "0.10.18" diff --git a/benchmark/elixir_2d_euler_vortex_unstructured.jl b/benchmark/elixir_2d_euler_vortex_unstructured.jl index 082b6648abf..43e4b6559de 100644 --- a/benchmark/elixir_2d_euler_vortex_unstructured.jl +++ b/benchmark/elixir_2d_euler_vortex_unstructured.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -49,11 +48,9 @@ end initial_condition = initial_condition_isentropic_vortex solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) -default_mesh_file = joinpath(@__DIR__, "mesh_uniform_cartesian.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/ranocha/f4ea19ba3b62348968c971db43d7798b/raw/a506abb9479c020920cf6068c142670fc1a9aadc/mesh_uniform_cartesian.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/ranocha/f4ea19ba3b62348968c971db43d7798b/raw/a506abb9479c020920cf6068c142670fc1a9aadc/mesh_uniform_cartesian.mesh", + joinpath(@__DIR__, "mesh_uniform_cartesian.mesh")) + mesh = UnstructuredMesh2D(mesh_file, periodicity = true) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/docs/make.jl b/docs/make.jl index dee87371bd1..7fce3b31e24 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -108,6 +108,7 @@ makedocs( ], "Time integration" => "time_integration.md", "Callbacks" => "callbacks.md", + "Coupling" => "multi-physics_coupling.md" ], "Advanced topics & developers" => [ "Conventions" =>"conventions.md", diff --git a/docs/src/multi-physics_coupling.md b/docs/src/multi-physics_coupling.md new file mode 100644 index 00000000000..eec92bc21de --- /dev/null +++ b/docs/src/multi-physics_coupling.md @@ -0,0 +1,46 @@ +# [Multi-physics coupling](@id multi-physics-coupling) +A complex simulation can consist of different spatial domains in which +different equations are being solved, different numerical methods being used +or the grid structure is different. +One example would be a fluid in a tank and an extended hot plate attached to it. +We would then like to solve the Navier-Stokes equations in the fluid domain +and the heat conduction equations in the plate. +The coupling would happen at the interface through the exchange of thermal energy. + + +## Converter coupling +It may happen that the two systems to be coupled do not share any variables, but +share some of the physics. +In such a situation, the same physics is just represented in a different form and with +a different set of variables. +This is the case, for instance assuming two domains, if there is a fluid system in one domain +and a Vlasov system in the other domain. +In that case we would have variables representing distribution functions of +the Vlasov system on one side and variables representing the mechanical quantities, like density, +of the fluid system. +To translate the fields from one description to the other one needs to use +converter functions. +These functions need to be hand tailored by the user in the elixir file where each +pair of coupled systems requires two coupling functions, one for each direction. + +In the general case, we have a system $A$ with $m$ variables +$u_{A,i}, \: i = 1, \dots, m$ and another +system $B$ with $n$ variables $u_{B,j}, \: j = 1, \dots, n$. +We then define two coupling functions, one that transforms $u_A$ into $u_B$ +and one that goes the other way. + +In their minimal form they take the position vector $x$, state vector $u$ +and the equations of the two coupled systems +and return the transformed variables. +By passing the equations we can make use of their parameters, if they are required. +Examples can be seen in `examples/structured_2d_dgsem/elixir_advection_coupled.jl`. + + +## Warning about binary compatibility +Currently the coordinate values on the nodes can differ by machine precision when +simulating the mesh and when splitting the mesh in multiple domains. +This is an issue coming from the coordinate interpolation on the nodes. +As a result, running a simulation in a single system and in two coupled domains +may result in a difference of the order of the machine precision. +While this is not an issue for most practical problems, it is best to keep this in mind when comparing test runs. + diff --git a/examples/dgmulti_2d/elixir_euler_hohqmesh.jl b/examples/dgmulti_2d/elixir_euler_hohqmesh.jl index f534b5bc8ad..9b14a5c6827 100644 --- a/examples/dgmulti_2d/elixir_euler_hohqmesh.jl +++ b/examples/dgmulti_2d/elixir_euler_hohqmesh.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -30,12 +29,8 @@ dg = DGMulti(polydeg = 8, element_type = Quad(), approximation_type = SBP(), ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) - -default_mesh_file = joinpath(@__DIR__, "mesh_trixi_unstructured_mesh_docs.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/52056f1487853fab63b7f4ed7f171c80/raw/9d573387dfdbb8bce2a55db7246f4207663ac07f/mesh_trixi_unstructured_mesh_docs.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/52056f1487853fab63b7f4ed7f171c80/raw/9d573387dfdbb8bce2a55db7246f4207663ac07f/mesh_trixi_unstructured_mesh_docs.mesh", + joinpath(@__DIR__, "mesh_trixi_unstructured_mesh_docs.mesh")) mesh = DGMultiMesh(dg, mesh_file) diff --git a/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl index 0a50b3644f0..4bfb2d3e375 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_amr_unstructured_flag.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -31,10 +30,8 @@ Trixi.validate_faces(faces) mapping_flag = Trixi.transfinite_mapping(faces) # Unstructured mesh with 24 cells of the square domain [-1, 1]^n -mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + joinpath(@__DIR__, "square_unstructured_2.inp")) mesh = P4estMesh{2}(mesh_file, polydeg = 3, mapping = mapping_flag, diff --git a/examples/p4est_2d_dgsem/elixir_advection_unstructured_flag.jl b/examples/p4est_2d_dgsem/elixir_advection_unstructured_flag.jl index 37fcc547f60..1ab96925fe6 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_unstructured_flag.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_unstructured_flag.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -29,10 +28,8 @@ Trixi.validate_faces(faces) mapping_flag = Trixi.transfinite_mapping(faces) # Unstructured mesh with 24 cells of the square domain [-1, 1]^n -mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + joinpath(@__DIR__, "square_unstructured_2.inp")) mesh = P4estMesh{2}(mesh_file, polydeg = 3, mapping = mapping_flag, diff --git a/examples/p4est_2d_dgsem/elixir_euler_blast_wave_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_blast_wave_amr.jl index 0ca4fdc2eb7..5db5f74a686 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_blast_wave_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_blast_wave_amr.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -50,10 +49,8 @@ volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; solver = DGSEM(basis, surface_flux, volume_integral) # Unstructured mesh with 48 cells of the square domain [-1, 1]^n -mesh_file = joinpath(@__DIR__, "square_unstructured_1.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp", + joinpath(@__DIR__, "square_unstructured_1.inp")) mesh = P4estMesh{2}(mesh_file, polydeg = 3, initial_refinement_level = 1) diff --git a/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr.jl index 92928146d7b..fbc11e89185 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_double_mach_amr.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -99,11 +98,8 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) # Get the unstructured quad mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "abaqus_double_mach.inp") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/a0806ef0d03cf5ea221af523167b6e32/raw/61ed0eb017eb432d996ed119a52fb041fe363e8c/abaqus_double_mach.inp", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/a0806ef0d03cf5ea221af523167b6e32/raw/61ed0eb017eb432d996ed119a52fb041fe363e8c/abaqus_double_mach.inp", + joinpath(@__DIR__, "abaqus_double_mach.inp")) mesh = P4estMesh{2}(mesh_file) diff --git a/examples/p4est_2d_dgsem/elixir_euler_forward_step_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_forward_step_amr.jl index 0ec9fc222f2..654efd5e209 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_forward_step_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_forward_step_amr.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -104,11 +103,8 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) # Get the unstructured quad mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "abaqus_forward_step.inp") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/b346ee6aa5446687f128eab8b37d52a7/raw/cd1e1d43bebd8d2631a07caec45585ec8456ca4c/abaqus_forward_step.inp", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/b346ee6aa5446687f128eab8b37d52a7/raw/cd1e1d43bebd8d2631a07caec45585ec8456ca4c/abaqus_forward_step.inp", + joinpath(@__DIR__, "abaqus_forward_step.inp")) mesh = P4estMesh{2}(mesh_file) diff --git a/examples/p4est_2d_dgsem/elixir_euler_free_stream.jl b/examples/p4est_2d_dgsem/elixir_euler_free_stream.jl index 38307a7d781..ab11dc11567 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_free_stream.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_free_stream.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -31,10 +30,8 @@ end # Get the uncurved mesh from a file (downloads the file if not available locally) # Unstructured mesh with 48 cells of the square domain [-1, 1]^n -mesh_file = joinpath(@__DIR__, "square_unstructured_1.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp", + joinpath(@__DIR__, "square_unstructured_1.inp")) # Map the unstructured mesh with the mapping above mesh = P4estMesh{2}(mesh_file, polydeg = 3, mapping = mapping, initial_refinement_level = 1) diff --git a/examples/p4est_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl b/examples/p4est_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl index 09d018309a6..084fd699b8e 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -31,10 +30,8 @@ mapping_flag = Trixi.transfinite_mapping(faces) # Get the uncurved mesh from a file (downloads the file if not available locally) # Unstructured mesh with 24 cells of the square domain [-1, 1]^n -mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + joinpath(@__DIR__, "square_unstructured_2.inp")) mesh = P4estMesh{2}(mesh_file, polydeg = 3, mapping = mapping_flag, diff --git a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl index 36c5624ba97..76ee96d4766 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder.jl @@ -13,7 +13,6 @@ # # Keywords: supersonic flow, shock capturing, AMR, unstructured curved mesh, positivity preservation, compressible Euler, 2D -using Downloads: download using OrdinaryDiffEq using Trixi @@ -82,11 +81,8 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) # Get the unstructured quad mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "abaqus_cylinder_in_channel.inp") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/a08f78f6b185b63c3baeff911a63f628/raw/addac716ea0541f588b9d2bd3f92f643eb27b88f/abaqus_cylinder_in_channel.inp", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/a08f78f6b185b63c3baeff911a63f628/raw/addac716ea0541f588b9d2bd3f92f643eb27b88f/abaqus_cylinder_in_channel.inp", + joinpath(@__DIR__, "abaqus_cylinder_in_channel.inp")) mesh = P4estMesh{2}(mesh_file) diff --git a/examples/p4est_2d_dgsem/elixir_euler_wall_bc_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_wall_bc_amr.jl index 8b8d05bade8..75e60d0c78b 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_wall_bc_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_wall_bc_amr.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -39,11 +38,8 @@ solver = DGSEM(polydeg = 5, surface_flux = flux_lax_friedrichs, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # Get the unstructured quad mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "abaqus_gingerbread_man.inp") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/0e9e990a04b5105d1d2e3096a6e41272/raw/0d924b1d7e7d3cc1070a6cc22fe1d501687aa6dd/abaqus_gingerbread_man.inp", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/0e9e990a04b5105d1d2e3096a6e41272/raw/0d924b1d7e7d3cc1070a6cc22fe1d501687aa6dd/abaqus_gingerbread_man.inp", + joinpath(@__DIR__, "abaqus_gingerbread_man.inp")) mesh = P4estMesh{2}(mesh_file) diff --git a/examples/p4est_2d_dgsem/elixir_mhd_rotor.jl b/examples/p4est_2d_dgsem/elixir_mhd_rotor.jl index 380db487356..089e82580c9 100644 --- a/examples/p4est_2d_dgsem/elixir_mhd_rotor.jl +++ b/examples/p4est_2d_dgsem/elixir_mhd_rotor.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -69,10 +68,8 @@ function mapping_twist(xi, eta) return SVector(x, y) end -mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + joinpath(@__DIR__, "square_unstructured_2.inp")) mesh = P4estMesh{2}(mesh_file, polydeg = 4, diff --git a/examples/p4est_3d_dgsem/elixir_advection_amr_unstructured_curved.jl b/examples/p4est_3d_dgsem/elixir_advection_amr_unstructured_curved.jl index cd280cf5bf6..33afd2e030e 100644 --- a/examples/p4est_3d_dgsem/elixir_advection_amr_unstructured_curved.jl +++ b/examples/p4est_3d_dgsem/elixir_advection_amr_unstructured_curved.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -49,10 +48,8 @@ function mapping(xi, eta, zeta) end # Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + joinpath(@__DIR__, "cube_unstructured_2.inp")) # Mesh polydeg of 2 (half the solver polydeg) to ensure FSP (see above). mesh = P4estMesh{3}(mesh_file, polydeg = 2, diff --git a/examples/p4est_3d_dgsem/elixir_advection_unstructured_curved.jl b/examples/p4est_3d_dgsem/elixir_advection_unstructured_curved.jl index 6df9ac0b16a..83adcbf6a63 100644 --- a/examples/p4est_3d_dgsem/elixir_advection_unstructured_curved.jl +++ b/examples/p4est_3d_dgsem/elixir_advection_unstructured_curved.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -46,10 +45,8 @@ function mapping(xi, eta, zeta) end # Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", + joinpath(@__DIR__, "cube_unstructured_1.inp")) mesh = P4estMesh{3}(mesh_file, polydeg = 3, mapping = mapping, diff --git a/examples/p4est_3d_dgsem/elixir_euler_ec.jl b/examples/p4est_3d_dgsem/elixir_euler_ec.jl index d9d774a7ffc..91698545052 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_ec.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_ec.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -46,10 +45,8 @@ function mapping(xi_, eta_, zeta_) end # Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + joinpath(@__DIR__, "cube_unstructured_2.inp")) mesh = P4estMesh{3}(mesh_file, polydeg = 5, mapping = mapping, diff --git a/examples/p4est_3d_dgsem/elixir_euler_free_stream.jl b/examples/p4est_3d_dgsem/elixir_euler_free_stream.jl index 24a781ca59e..6406a38186b 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_free_stream.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_free_stream.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -47,10 +46,8 @@ function mapping(xi_, eta_, zeta_) end # Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", + joinpath(@__DIR__, "cube_unstructured_1.inp")) # Mesh polydeg of 2 (half the solver polydeg) to ensure FSP (see above). mesh = P4estMesh{3}(mesh_file, polydeg = 2, diff --git a/examples/p4est_3d_dgsem/elixir_euler_free_stream_extruded.jl b/examples/p4est_3d_dgsem/elixir_euler_free_stream_extruded.jl index f56fe3a429d..08307a449a7 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_free_stream_extruded.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_free_stream_extruded.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -36,10 +35,8 @@ function mapping(xi, eta_, zeta_) end # Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + joinpath(@__DIR__, "cube_unstructured_2.inp")) mesh = P4estMesh{3}(mesh_file, polydeg = 3, mapping = mapping, diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl index 0de22eaea40..e7ca0cad4ba 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -49,10 +48,8 @@ function mapping(xi, eta, zeta) end # Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", + joinpath(@__DIR__, "cube_unstructured_1.inp")) # Mesh polydeg of 2 (half the solver polydeg) to ensure FSP (see above). mesh = P4estMesh{3}(mesh_file, polydeg = 2, diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh.jl index 0fa3a28fe8b..7d81d6739bf 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -19,11 +18,8 @@ boundary_conditions = Dict(:Bottom => boundary_condition, solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs) # Unstructured 3D half circle mesh from HOHQMesh -default_mesh_file = joinpath(@__DIR__, "abaqus_half_circle_3d.inp") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/11461efbfb02c42e06aca338b3d0b645/raw/81deeb1ebc4945952c30af5bb75fe222a18d975c/abaqus_half_circle_3d.inp", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/11461efbfb02c42e06aca338b3d0b645/raw/81deeb1ebc4945952c30af5bb75fe222a18d975c/abaqus_half_circle_3d.inp", + joinpath(@__DIR__, "abaqus_half_circle_3d.inp")) mesh = P4estMesh{3}(mesh_file, initial_refinement_level = 0) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index 2a56d23f4c0..43b68f21b03 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -2,31 +2,38 @@ using OrdinaryDiffEq using Trixi ############################################################################### -# Coupled semidiscretization of two linear advection systems, which are connected periodically +# Coupled semidiscretization of four linear advection systems using converter functions such that +# they are also coupled across the domain boundaries to generate a periodic system. # -# In this elixir, we have a square domain that is divided into a left half and a right half. On each -# half of the domain, a completely independent SemidiscretizationHyperbolic is created for the -# linear advection equations. The two systems are coupled in the x-direction and have periodic -# boundaries in the y-direction. For a high-level overview, see also the figure below: +# In this elixir, we have a square domain that is divided into a upper-left, lower-left, +# upper-right and lower-right quarter. On each quarter +# of the domain, a completely independent SemidiscretizationHyperbolic is created for the +# linear advection equations. The four systems are coupled in the x and y-direction. +# For a high-level overview, see also the figure below: # # (-1, 1) ( 1, 1) # ┌────────────────────┬────────────────────┐ -# │ ↑ periodic ↑ │ ↑ periodic ↑ │ -# │ │ │ +# │ ↑ coupled ↑ │ ↑ coupled ↑ │ # │ │ │ # │ ========= │ ========= │ # │ system #1 │ system #2 │ # │ ========= │ ========= │ # │ │ │ +# │<-- coupled │<-- coupled │ +# │ coupled -->│ coupled -->│ # │ │ │ +# │ ↓ coupled ↓ │ ↓ coupled ↓ │ +# ├────────────────────┼────────────────────┤ +# │ ↑ coupled ↑ │ ↑ coupled ↑ │ # │ │ │ +# │ ========= │ ========= │ +# │ system #3 │ system #4 │ +# │ ========= │ ========= │ # │ │ │ -# │ coupled -->│<-- coupled │ -# │ │ │ -# │<-- coupled │ coupled -->│ -# │ │ │ +# │<-- coupled │<-- coupled │ +# │ coupled -->│ coupled -->│ # │ │ │ -# │ ↓ periodic ↓ │ ↓ periodic ↓ │ +# │ ↓ coupled ↓ │ ↓ coupled ↓ │ # └────────────────────┴────────────────────┘ # (-1, -1) ( 1, -1) @@ -36,60 +43,135 @@ equations = LinearScalarAdvectionEquation2D(advection_velocity) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) -# First mesh is the left half of a [-1,1]^2 square -coordinates_min1 = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +# This will be the number of elements for each quarter/semidiscretization. +cells_per_dimension = (8, 8) + +########### +# system #1 +########### + +coordinates_min1 = (-1.0, 0.0) # minimum coordinates (min(x), min(y)) coordinates_max1 = (0.0, 1.0) # maximum coordinates (max(x), max(y)) -# Define identical resolution as a variable such that it is easier to change from `trixi_include` -cells_per_dimension = (8, 16) +mesh1 = StructuredMesh(cells_per_dimension, coordinates_min1, coordinates_max1) -cells_per_dimension1 = cells_per_dimension +# Define the coupling functions +coupling_function12 = (x, u, equations_other, equations_own) -> u +coupling_function13 = (x, u, equations_other, equations_own) -> u -mesh1 = StructuredMesh(cells_per_dimension1, coordinates_min1, coordinates_max1) +# Define the coupling boundary conditions and the system it is coupled to. +boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, + coupling_function12) +boundary_conditions_x_pos1 = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, + coupling_function12) +boundary_conditions_y_neg1 = BoundaryConditionCoupled(3, (:i_forward, :end), Float64, + coupling_function13) +boundary_conditions_y_pos1 = BoundaryConditionCoupled(3, (:i_forward, :begin), Float64, + coupling_function13) # A semidiscretization collects data structures and functions for the spatial discretization semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, solver, - boundary_conditions = ( - # Connect left boundary with right boundary of right mesh - x_neg = BoundaryConditionCoupled(2, - (:end, - :i_forward), - Float64), - # Connect right boundary with left boundary of right mesh - x_pos = BoundaryConditionCoupled(2, - (:begin, - :i_forward), - Float64), - y_neg = boundary_condition_periodic, - y_pos = boundary_condition_periodic)) - -# Second mesh is the right half of a [-1,1]^2 square -coordinates_min2 = (0.0, -1.0) # minimum coordinates (min(x), min(y)) + boundary_conditions = (x_neg = boundary_conditions_x_neg1, + x_pos = boundary_conditions_x_pos1, + y_neg = boundary_conditions_y_neg1, + y_pos = boundary_conditions_y_pos1)) + +########### +# system #2 +########### + +coordinates_min2 = (0.0, 0.0) # minimum coordinates (min(x), min(y)) coordinates_max2 = (1.0, 1.0) # maximum coordinates (max(x), max(y)) -cells_per_dimension2 = cells_per_dimension +mesh2 = StructuredMesh(cells_per_dimension, coordinates_min2, coordinates_max2) -mesh2 = StructuredMesh(cells_per_dimension2, coordinates_min2, coordinates_max2) +# Define the coupling functions +coupling_function21 = (x, u, equations_other, equations_own) -> u +coupling_function24 = (x, u, equations_other, equations_own) -> u +# Define the coupling boundary conditions and the system it is coupled to. +boundary_conditions_x_neg2 = BoundaryConditionCoupled(1, (:end, :i_forward), Float64, + coupling_function21) +boundary_conditions_x_pos2 = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, + coupling_function21) +boundary_conditions_y_neg2 = BoundaryConditionCoupled(4, (:i_forward, :end), Float64, + coupling_function24) +boundary_conditions_y_pos2 = BoundaryConditionCoupled(4, (:i_forward, :begin), Float64, + coupling_function24) + +# A semidiscretization collects data structures and functions for the spatial discretization semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, solver, - boundary_conditions = ( - # Connect left boundary with right boundary of left mesh - x_neg = BoundaryConditionCoupled(1, - (:end, - :i_forward), - Float64), - # Connect right boundary with left boundary of left mesh - x_pos = BoundaryConditionCoupled(1, - (:begin, - :i_forward), - Float64), - y_neg = boundary_condition_periodic, - y_pos = boundary_condition_periodic)) - -# Create a semidiscretization that bundles semi1 and semi2 -semi = SemidiscretizationCoupled(semi1, semi2) + boundary_conditions = (x_neg = boundary_conditions_x_neg2, + x_pos = boundary_conditions_x_pos2, + y_neg = boundary_conditions_y_neg2, + y_pos = boundary_conditions_y_pos2)) + +########### +# system #3 +########### + +coordinates_min3 = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max3 = (0.0, 0.0) # maximum coordinates (max(x), max(y)) + +mesh3 = StructuredMesh(cells_per_dimension, coordinates_min3, coordinates_max3) + +# Define the coupling functions +coupling_function34 = (x, u, equations_other, equations_own) -> u +coupling_function31 = (x, u, equations_other, equations_own) -> u + +# Define the coupling boundary conditions and the system it is coupled to. +boundary_conditions_x_neg3 = BoundaryConditionCoupled(4, (:end, :i_forward), Float64, + coupling_function34) +boundary_conditions_x_pos3 = BoundaryConditionCoupled(4, (:begin, :i_forward), Float64, + coupling_function34) +boundary_conditions_y_neg3 = BoundaryConditionCoupled(1, (:i_forward, :end), Float64, + coupling_function31) +boundary_conditions_y_pos3 = BoundaryConditionCoupled(1, (:i_forward, :begin), Float64, + coupling_function31) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi3 = SemidiscretizationHyperbolic(mesh3, equations, initial_condition_convergence_test, + solver, + boundary_conditions = (x_neg = boundary_conditions_x_neg3, + x_pos = boundary_conditions_x_pos3, + y_neg = boundary_conditions_y_neg3, + y_pos = boundary_conditions_y_pos3)) + +########### +# system #4 +########### + +coordinates_min4 = (0.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max4 = (1.0, 0.0) # maximum coordinates (max(x), max(y)) + +mesh4 = StructuredMesh(cells_per_dimension, coordinates_min4, coordinates_max4) + +# Define the coupling functions +coupling_function43 = (x, u, equations_other, equations_own) -> u +coupling_function42 = (x, u, equations_other, equations_own) -> u + +# Define the coupling boundary conditions and the system it is coupled to. +boundary_conditions_x_neg4 = BoundaryConditionCoupled(3, (:end, :i_forward), Float64, + coupling_function43) +boundary_conditions_x_pos4 = BoundaryConditionCoupled(3, (:begin, :i_forward), Float64, + coupling_function43) +boundary_conditions_y_neg4 = BoundaryConditionCoupled(2, (:i_forward, :end), Float64, + coupling_function42) +boundary_conditions_y_pos4 = BoundaryConditionCoupled(2, (:i_forward, :begin), Float64, + coupling_function42) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi4 = SemidiscretizationHyperbolic(mesh4, equations, initial_condition_convergence_test, + solver, + boundary_conditions = (x_neg = boundary_conditions_x_neg4, + x_pos = boundary_conditions_x_pos4, + y_neg = boundary_conditions_y_neg4, + y_pos = boundary_conditions_y_pos4)) + +# Create a semidiscretization that bundles all the semidiscretizations. +semi = SemidiscretizationCoupled(semi1, semi2, semi3, semi4) ############################################################################### # ODE solvers, callbacks etc. @@ -104,7 +186,10 @@ summary_callback = SummaryCallback() # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results analysis_callback1 = AnalysisCallback(semi1, interval = 100) analysis_callback2 = AnalysisCallback(semi2, interval = 100) -analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2) +analysis_callback3 = AnalysisCallback(semi3, interval = 100) +analysis_callback4 = AnalysisCallback(semi4, interval = 100) +analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2, + analysis_callback3, analysis_callback4) # The SaveSolutionCallback allows to save the solution to a file in regular intervals save_solution = SaveSolutionCallback(interval = 100, diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl index 61dd252fd83..a6a56aa807c 100644 --- a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl index f285d24fc6c..0923e328487 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl @@ -1,4 +1,3 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -31,10 +30,8 @@ Trixi.validate_faces(faces) mapping_flag = Trixi.transfinite_mapping(faces) # Unstructured mesh with 24 cells of the square domain [-1, 1]^n -mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + joinpath(@__DIR__, "square_unstructured_2.inp")) # INP mesh files are only support by p4est. Hence, we # create a p4est connecvity object first from which diff --git a/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl index 5ba1ab15489..ba8f1b59b80 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl @@ -1,4 +1,3 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -28,10 +27,8 @@ Trixi.validate_faces(faces) mapping_flag = Trixi.transfinite_mapping(faces) # Unstructured mesh with 24 cells of the square domain [-1, 1]^n. -mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + joinpath(@__DIR__, "square_unstructured_2.inp")) # INP mesh files are only support by p4est. Hence, we # create a p4est connecvity object first from which diff --git a/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl b/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl index 37d15f38566..5e6c4193c50 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl @@ -1,4 +1,3 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -30,10 +29,8 @@ end # Get the uncurved mesh from a file (downloads the file if not available locally) # Unstructured mesh with 48 cells of the square domain [-1, 1]^n -mesh_file = joinpath(@__DIR__, "square_unstructured_1.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp", + joinpath(@__DIR__, "square_unstructured_1.inp")) # INP mesh files are only support by p4est. Hence, we # create a p4est connecvity object first from which diff --git a/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl index bcc1abc560e..e496eb76729 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl @@ -1,4 +1,3 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -30,10 +29,8 @@ mapping_flag = Trixi.transfinite_mapping(faces) # Get the uncurved mesh from a file (downloads the file if not available locally) # Unstructured mesh with 24 cells of the square domain [-1, 1]^n -mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + joinpath(@__DIR__, "square_unstructured_2.inp")) # INP mesh files are only support by p4est. Hence, we # create a p4est connecvity object first from which diff --git a/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl b/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl index adb154948fb..ff2e40ae607 100644 --- a/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl +++ b/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl @@ -1,4 +1,3 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -68,10 +67,8 @@ function mapping_twist(xi, eta) return SVector(x, y) end -mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + joinpath(@__DIR__, "square_unstructured_2.inp")) # INP mesh files are only support by p4est. Hence, we # create a p4est connecvity object first from which diff --git a/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl index 617736afbdd..e7c0f4b7318 100644 --- a/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl +++ b/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl @@ -1,4 +1,3 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -47,11 +46,9 @@ function mapping(xi, eta, zeta) return SVector(5 * x, 5 * y, 5 * z) end -# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", - mesh_file) +# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3. +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + joinpath(@__DIR__, "cube_unstructured_2.inp")) # INP mesh files are only support by p4est. Hence, we # create a p4est connectivity object first from which diff --git a/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl index df358435c9a..ee27ee117fe 100644 --- a/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl +++ b/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl @@ -1,4 +1,3 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -45,10 +44,8 @@ function mapping(xi, eta, zeta) end # Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", + joinpath(@__DIR__, "cube_unstructured_1.inp")) # INP mesh files are only support by p4est. Hence, we # create a p4est connectivity object first from which diff --git a/examples/t8code_3d_dgsem/elixir_euler_ec.jl b/examples/t8code_3d_dgsem/elixir_euler_ec.jl index 07745c3ac56..b720bfcd375 100644 --- a/examples/t8code_3d_dgsem/elixir_euler_ec.jl +++ b/examples/t8code_3d_dgsem/elixir_euler_ec.jl @@ -1,4 +1,3 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -45,10 +44,8 @@ function mapping(xi_, eta_, zeta_) end # Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + joinpath(@__DIR__, "cube_unstructured_2.inp")) # INP mesh files are only support by p4est. Hence, we # create a p4est connectivity object first from which diff --git a/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl b/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl index e135d464810..b70a6091adf 100644 --- a/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl +++ b/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl @@ -1,4 +1,3 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -46,10 +45,8 @@ function mapping(xi_, eta_, zeta_) end # Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", + joinpath(@__DIR__, "cube_unstructured_1.inp")) # INP mesh files are only support by p4est. Hence, we # create a p4est connectivity object first from which diff --git a/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl b/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl index d129b59826e..6ae38d20b5a 100644 --- a/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl +++ b/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl @@ -1,4 +1,3 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -35,10 +34,8 @@ function mapping(xi, eta_, zeta_) end # Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + joinpath(@__DIR__, "cube_unstructured_2.inp")) # INP mesh files are only support by p4est. Hence, we # create a p4est connecvity object first from which diff --git a/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl index d4664522bea..6856be36ea1 100644 --- a/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl +++ b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl @@ -1,4 +1,3 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -48,10 +47,8 @@ function mapping(xi, eta, zeta) end # Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 -mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") -isfile(mesh_file) || - download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", - mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", + joinpath(@__DIR__, "cube_unstructured_1.inp")) # INP mesh files are only support by p4est. Hence, we # create a p4est connecvity object first from which diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl new file mode 100644 index 00000000000..1817672778a --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl @@ -0,0 +1,91 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +""" + initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D) + +A version of the classical Kelvin-Helmholtz instability based on +- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021) + A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations + of the Euler Equations + [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017) +""" +function initial_condition_kelvin_helmholtz_instability(x, t, + equations::CompressibleEulerEquations2D) + # change discontinuity to tanh + # typical resolution 128^2, 256^2 + # domain size is [-1,+1]^2 + slope = 15 + amplitude = 0.02 + B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5) + rho = 0.5 + 0.75 * B + v1 = 0.5 * (B - 1) + v2 = 0.1 * sin(2 * pi * x[1]) + p = 1.0 + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_kelvin_helmholtz_instability + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) + +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons = ["rho"], + positivity_variables_nonlinear = [pressure]) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 5, + n_cells_max = 100_000) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3.7) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +save_restart = SaveRestartCallback(interval = 1000, + save_final_restart = true) + +stepsize_callback = StepsizeCallback(cfl = 0.7) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback, + save_restart, save_solution) + +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false)) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl index fe9ad92467f..74d0370647a 100644 --- a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl @@ -22,7 +22,7 @@ function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D) r = sqrt(x[1]^2 + x[2]^2) pmax = 10.0 - pmin = 1.0 + pmin = 0.01 rhomax = 1.0 rhomin = 0.01 if r <= 0.09 @@ -52,7 +52,8 @@ basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; positivity_variables_cons = ["rho"], - positivity_correction_factor = 0.5) + positivity_variables_nonlinear = [pressure], + positivity_correction_factor = 0.1) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) @@ -84,7 +85,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -cfl = 0.5 +cfl = 0.4 stepsize_callback = StepsizeCallback(cfl = cfl) glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl) diff --git a/examples/unstructured_2d_dgsem/elixir_acoustics_gauss_wall.jl b/examples/unstructured_2d_dgsem/elixir_acoustics_gauss_wall.jl index 8f8e867dca8..9741430d11c 100644 --- a/examples/unstructured_2d_dgsem/elixir_acoustics_gauss_wall.jl +++ b/examples/unstructured_2d_dgsem/elixir_acoustics_gauss_wall.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -14,11 +13,8 @@ equations = AcousticPerturbationEquations2D(v_mean_global = (0.0, -0.5), solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs) # Create unstructured quadrilateral mesh from a file -default_mesh_file = joinpath(@__DIR__, "mesh_five_circles_in_circle.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/3c79baad6b4d73bb26ec6420b5d16f45/raw/22aefc4ec2107cf0bffc40e81dfbc52240c625b1/mesh_five_circles_in_circle.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/3c79baad6b4d73bb26ec6420b5d16f45/raw/22aefc4ec2107cf0bffc40e81dfbc52240c625b1/mesh_five_circles_in_circle.mesh", + joinpath(@__DIR__, "mesh_five_circles_in_circle.mesh")) mesh = UnstructuredMesh2D(mesh_file) diff --git a/examples/unstructured_2d_dgsem/elixir_advection_basic.jl b/examples/unstructured_2d_dgsem/elixir_advection_basic.jl index afef6c2c38f..c0ee453344d 100644 --- a/examples/unstructured_2d_dgsem/elixir_advection_basic.jl +++ b/examples/unstructured_2d_dgsem/elixir_advection_basic.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -16,12 +15,8 @@ solver = DGSEM(polydeg = 6, surface_flux = flux_lax_friedrichs) ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) - -default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", + joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl index cd6a1995757..f8976120d53 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_basic.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_basic.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -25,12 +24,8 @@ solver = DGSEM(polydeg = 8, surface_flux = flux_lax_friedrichs) ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) - -default_mesh_file = joinpath(@__DIR__, "mesh_trixi_unstructured_mesh_docs.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/52056f1487853fab63b7f4ed7f171c80/raw/9d573387dfdbb8bce2a55db7246f4207663ac07f/mesh_trixi_unstructured_mesh_docs.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/52056f1487853fab63b7f4ed7f171c80/raw/9d573387dfdbb8bce2a55db7246f4207663ac07f/mesh_trixi_unstructured_mesh_docs.mesh", + joinpath(@__DIR__, "mesh_trixi_unstructured_mesh_docs.mesh")) mesh = UnstructuredMesh2D(mesh_file) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_ec.jl b/examples/unstructured_2d_dgsem/elixir_euler_ec.jl index 0f53aa62a18..58b4d9a1dd2 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_ec.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_ec.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -19,12 +18,8 @@ solver = DGSEM(polydeg = 6, surface_flux = flux_ranocha, ############################################################################### # Get the curved quad mesh from a file - -default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", + joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_free_stream.jl b/examples/unstructured_2d_dgsem/elixir_euler_free_stream.jl index a2fec1a320a..f266a3de0b2 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_free_stream.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_free_stream.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -26,12 +25,8 @@ solver = DGSEM(polydeg = 6, surface_flux = flux_hll) ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) - -default_mesh_file = joinpath(@__DIR__, "mesh_gingerbread_man.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/2c6440b5f8a57db131061ad7aa78ee2b/raw/1f89fdf2c874ff678c78afb6fe8dc784bdfd421f/mesh_gingerbread_man.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/2c6440b5f8a57db131061ad7aa78ee2b/raw/1f89fdf2c874ff678c78afb6fe8dc784bdfd421f/mesh_gingerbread_man.mesh", + joinpath(@__DIR__, "mesh_gingerbread_man.mesh")) mesh = UnstructuredMesh2D(mesh_file) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_periodic.jl b/examples/unstructured_2d_dgsem/elixir_euler_periodic.jl index afd177f0740..e640001ad7f 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_periodic.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_periodic.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -19,12 +18,8 @@ solver = DGSEM(polydeg = 6, surface_flux = FluxRotated(flux_hll)) ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) - -default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", + joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl b/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl index e1cc0932969..06053273b74 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl @@ -55,11 +55,8 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) # Get the curved quad mesh from a file -default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", + joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_wall_bc.jl b/examples/unstructured_2d_dgsem/elixir_euler_wall_bc.jl index b2abefe7eeb..65e5eb51ce6 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_wall_bc.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_wall_bc.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -40,12 +39,8 @@ solver = DGSEM(polydeg = 4, surface_flux = flux_hll) ############################################################################### # Get the curved quad mesh from a file - -default_mesh_file = joinpath(@__DIR__, "mesh_box_around_circle.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/8b9b11a1eedfa54b215c122c3d17b271/raw/0d2b5d98c87e67a6f384693a8b8e54b4c9fcbf3d/mesh_box_around_circle.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8b9b11a1eedfa54b215c122c3d17b271/raw/0d2b5d98c87e67a6f384693a8b8e54b4c9fcbf3d/mesh_box_around_circle.mesh", + joinpath(@__DIR__, "mesh_box_around_circle.mesh")) mesh = UnstructuredMesh2D(mesh_file) diff --git a/examples/unstructured_2d_dgsem/elixir_mhd_alfven_wave.jl b/examples/unstructured_2d_dgsem/elixir_mhd_alfven_wave.jl index 3ed3e828ca8..0c7152a6ea0 100644 --- a/examples/unstructured_2d_dgsem/elixir_mhd_alfven_wave.jl +++ b/examples/unstructured_2d_dgsem/elixir_mhd_alfven_wave.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -17,12 +16,9 @@ solver = DGSEM(polydeg = 7, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # Get the unstructured quad mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", - default_mesh_file) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) -mesh_file = default_mesh_file mesh = UnstructuredMesh2D(mesh_file, periodicity = true) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/unstructured_2d_dgsem/elixir_mhd_ec.jl b/examples/unstructured_2d_dgsem/elixir_mhd_ec.jl index a40f92cac02..805934e305d 100644 --- a/examples/unstructured_2d_dgsem/elixir_mhd_ec.jl +++ b/examples/unstructured_2d_dgsem/elixir_mhd_ec.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -35,11 +34,8 @@ solver = DGSEM(polydeg = 6, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # Get the unstructured quad mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_dirichlet.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_dirichlet.jl index 1148f25fae3..df1a69192ce 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_dirichlet.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_dirichlet.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -38,11 +37,8 @@ solver = DGSEM(polydeg = 4, surface_flux = (flux_hll, flux_nonconservative_fjord # This setup is for the curved, split form well-balancedness testing # Get the unstructured quad mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "mesh_outer_circle.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/9beddd9cd00e2a0a15865129eeb24928/raw/be71e67fa48bc4e1e97f5f6cd77c3ed34c6ba9be/mesh_outer_circle.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/9beddd9cd00e2a0a15865129eeb24928/raw/be71e67fa48bc4e1e97f5f6cd77c3ed34c6ba9be/mesh_outer_circle.mesh", + joinpath(@__DIR__, "mesh_outer_circle.mesh")) mesh = UnstructuredMesh2D(mesh_file) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl index 8e9d396d826..9122fb8287d 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -25,11 +24,8 @@ solver = DGSEM(polydeg = 6, # This setup is for the curved, split form entropy conservation testing (needs periodic BCs) # Get the unstructured quad mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl index 94202b81df0..98408db5a78 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -34,11 +33,8 @@ solver = DGSEM(basis, surface_flux, volume_integral) # This setup is for the curved, split form entropy conservation testing (needs periodic BCs) # Get the unstructured quad mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_source_terms.jl index 07668688406..a7aa5808955 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_source_terms.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_source_terms.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -23,11 +22,8 @@ solver = DGSEM(polydeg = 6, surface_flux = surface_flux, # This setup is for the curved, split form convergence test on a periodic domain # Get the unstructured quad mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl new file mode 100644 index 00000000000..6164f9d4a55 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl @@ -0,0 +1,140 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the shallow water equations +# +# TODO: TrixiShallowWater: wet/dry example elixir + +equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 1.875, + threshold_limiter = 1e-12, threshold_wet = 1e-14) + +""" + initial_condition_three_mounds(x, t, equations::ShallowWaterEquations2D) + +Initial condition simulating a dam break. The bottom topography is given by one large and two smaller +mounds. The mounds are flooded by the water for t > 0. To smooth the discontinuity, a logistic function +is applied. + +The initial conditions is taken from Section 6.3 of the paper: +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018) + An entropy stable discontinuous Galerkin method for the shallow water equations on + curvilinear meshes with wet/dry fronts accelerated by GPUs\n + [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038) +""" +function initial_condition_three_mounds(x, t, equations::ShallowWaterEquations2D) + + # Set the background values + v1 = 0.0 + v2 = 0.0 + + x1, x2 = x + M_1 = 1 - 0.1 * sqrt((x1 - 30.0)^2 + (x2 - 22.5)^2) + M_2 = 1 - 0.1 * sqrt((x1 - 30.0)^2 + (x2 - 7.5)^2) + M_3 = 2.8 - 0.28 * sqrt((x1 - 47.5)^2 + (x2 - 15.0)^2) + + b = max(0.0, M_1, M_2, M_3) + + # use a logistic function to transfer water height value smoothly + L = equations.H0 # maximum of function + x0 = 8 # center point of function + k = -75.0 # sharpness of transfer + + H = max(b, L / (1.0 + exp(-k * (x1 - x0)))) + + # Avoid division by zero by adjusting the initial condition with a small dry state threshold + # that defaults to 500*eps() ≈ 1e-13 in double precision and is set in the constructor above + # for the ShallowWaterEquations struct. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_three_mounds + +function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t, + surface_flux_function, + equations::ShallowWaterEquations2D) + # Impulse and bottom from inside, height from external state + u_outer = SVector(equations.threshold_wet, u_inner[2], u_inner[3], u_inner[4]) + + # calculate the boundary flux + flux = surface_flux_function(u_inner, u_outer, normal_direction, equations) + + return flux +end + +boundary_conditions = Dict(:Bottom => boundary_condition_slip_wall, + :Top => boundary_condition_slip_wall, + :Right => boundary_condition_outflow, + :Left => boundary_condition_slip_wall) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, + hydrostatic_reconstruction_chen_noelle), + flux_nonconservative_chen_noelle) + +basis = LobattoLegendreBasis(4) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Get the unstructured quad mesh from a file (downloads the file if not available locally) + +default_meshfile = joinpath(@__DIR__, "mesh_three_mound.mesh") + +isfile(default_meshfile) || + download("https://gist.githubusercontent.com/svengoldberg/c3c87fecb3fc6e46be7f0d1c7cb35f83/raw/e817ecd9e6c4686581d63c46128f9b6468d396d3/mesh_three_mound.mesh", + default_meshfile) + +meshfile = default_meshfile + +mesh = UnstructuredMesh2D(meshfile) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solver + +tspan = (0.0, 20.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl new file mode 100644 index 00000000000..e69de29bb2d diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl new file mode 100644 index 00000000000..e69de29bb2d diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl new file mode 100644 index 00000000000..e69de29bb2d diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_wall_bc_shockcapturing.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_wall_bc_shockcapturing.jl index 76b9642d595..f115113ed27 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_wall_bc_shockcapturing.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_wall_bc_shockcapturing.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -55,12 +54,8 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, ############################################################################### # Get the unstructured quad mesh from a file (downloads the file if not available locally) - -default_mesh_file = joinpath(@__DIR__, "mesh_outer_circle.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/9beddd9cd00e2a0a15865129eeb24928/raw/be71e67fa48bc4e1e97f5f6cd77c3ed34c6ba9be/mesh_outer_circle.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/9beddd9cd00e2a0a15865129eeb24928/raw/be71e67fa48bc4e1e97f5f6cd77c3ed34c6ba9be/mesh_outer_circle.mesh", + joinpath(@__DIR__, "mesh_outer_circle.mesh")) mesh = UnstructuredMesh2D(mesh_file) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl index bf4d0be682a..6bad3a77f03 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -36,13 +35,9 @@ solver = DGSEM(polydeg = 6, surface_flux = surface_flux, ############################################################################### # This setup is for the curved, split form well-balancedness testing - # Get the unstructured quad mesh from a file (downloads the file if not available locally) -default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) diff --git a/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl b/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl index c181203e7a4..fe7e708f3b3 100644 --- a/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl +++ b/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -21,12 +20,8 @@ solver = FDSBP(D_SBP, ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) - -default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", + joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl index 7ada50c0c65..25a81c16bf9 100644 --- a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl +++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -33,12 +32,8 @@ solver = FDSBP(D_SBP, ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) - -default_mesh_file = joinpath(@__DIR__, "mesh_gingerbread_man.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/2c6440b5f8a57db131061ad7aa78ee2b/raw/1f89fdf2c874ff678c78afb6fe8dc784bdfd421f/mesh_gingerbread_man.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/2c6440b5f8a57db131061ad7aa78ee2b/raw/1f89fdf2c874ff678c78afb6fe8dc784bdfd421f/mesh_gingerbread_man.mesh", + joinpath(@__DIR__, "mesh_gingerbread_man.mesh")) mesh = UnstructuredMesh2D(mesh_file) diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl index edcd221bf59..5f11d41ad5c 100644 --- a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl +++ b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -22,12 +21,8 @@ solver = FDSBP(D_SBP, ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) - -default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh") -isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", - default_mesh_file) -mesh_file = default_mesh_file +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", + joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")) mesh = UnstructuredMesh2D(mesh_file, periodicity = true) diff --git a/src/Trixi.jl b/src/Trixi.jl index e6373380d39..65ace2f238e 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -39,6 +39,7 @@ import SciMLBase: get_du, get_tmp_cache, u_modified!, get_proposed_dt, set_proposed_dt!, terminate!, remake, add_tstop!, has_tstop, first_tstop +using Downloads: Downloads using CodeTracking: CodeTracking using ConstructionBase: ConstructionBase using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 1f7d30d6aa8..92da9a5ba8b 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -345,4 +345,27 @@ function register_error_hints() return nothing end + +""" + Trixi.download(src_url, file_path) + +Download a file from given `src_url` to given `file_path` if +`file_path` is not already a file. This function just returns +`file_path`. +This is a small wrapper of `Downloads.download(src_url, file_path)` +that avoids race conditions when multiple MPI ranks are used. +""" +function download(src_url, file_path) + # Note that `mpi_isroot()` is also `true` if running + # in serial (without MPI). + if mpi_isroot() + isfile(file_path) || Downloads.download(src_url, file_path) + end + + if mpi_isparallel() + MPI.Barrier(mpi_comm()) + end + + return file_path +end end # @muladd diff --git a/src/auxiliary/t8code.jl b/src/auxiliary/t8code.jl index db01476bb86..7c1399fc803 100644 --- a/src/auxiliary/t8code.jl +++ b/src/auxiliary/t8code.jl @@ -46,230 +46,6 @@ function init_t8code() return nothing end -function trixi_t8_unref_forest(forest) - t8_forest_unref(Ref(forest)) -end - -function t8_free(ptr) - T8code.Libt8.sc_free(t8_get_package_id(), ptr) -end - -function trixi_t8_count_interfaces(forest) - # Check that forest is a committed, that is valid and usable, forest. - @assert t8_forest_is_committed(forest) != 0 - - # Get the number of local elements of forest. - num_local_elements = t8_forest_get_local_num_elements(forest) - # Get the number of ghost elements of forest. - num_ghost_elements = t8_forest_get_num_ghosts(forest) - # Get the number of trees that have elements of this process. - num_local_trees = t8_forest_get_num_local_trees(forest) - - current_index = t8_locidx_t(0) - - local_num_conform = 0 - local_num_mortars = 0 - local_num_boundary = 0 - - for itree in 0:(num_local_trees - 1) - tree_class = t8_forest_get_tree_class(forest, itree) - eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) - - # Get the number of elements of this tree. - num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) - - for ielement in 0:(num_elements_in_tree - 1) - element = t8_forest_get_element_in_tree(forest, itree, ielement) - - level = t8_element_level(eclass_scheme, element) - - num_faces = t8_element_num_faces(eclass_scheme, element) - - for iface in 0:(num_faces - 1) - pelement_indices_ref = Ref{Ptr{t8_locidx_t}}() - pneighbor_leafs_ref = Ref{Ptr{Ptr{t8_element}}}() - pneigh_scheme_ref = Ref{Ptr{t8_eclass_scheme}}() - - dual_faces_ref = Ref{Ptr{Cint}}() - num_neighbors_ref = Ref{Cint}() - - forest_is_balanced = Cint(1) - - t8_forest_leaf_face_neighbors(forest, itree, element, - pneighbor_leafs_ref, iface, dual_faces_ref, - num_neighbors_ref, - pelement_indices_ref, pneigh_scheme_ref, - forest_is_balanced) - - num_neighbors = num_neighbors_ref[] - neighbor_ielements = unsafe_wrap(Array, pelement_indices_ref[], - num_neighbors) - neighbor_leafs = unsafe_wrap(Array, pneighbor_leafs_ref[], num_neighbors) - neighbor_scheme = pneigh_scheme_ref[] - - if num_neighbors > 0 - neighbor_level = t8_element_level(neighbor_scheme, neighbor_leafs[1]) - - # Conforming interface: The second condition ensures we only visit the interface once. - if level == neighbor_level && current_index <= neighbor_ielements[1] - local_num_conform += 1 - elseif level < neighbor_level - local_num_mortars += 1 - end - else - local_num_boundary += 1 - end - - t8_free(dual_faces_ref[]) - t8_free(pneighbor_leafs_ref[]) - t8_free(pelement_indices_ref[]) - end # for - - current_index += 1 - end # for - end # for - - return (interfaces = local_num_conform, - mortars = local_num_mortars, - boundaries = local_num_boundary) -end - -function trixi_t8_fill_mesh_info(forest, elements, interfaces, mortars, boundaries, - boundary_names) - # Check that forest is a committed, that is valid and usable, forest. - @assert t8_forest_is_committed(forest) != 0 - - # Get the number of local elements of forest. - num_local_elements = t8_forest_get_local_num_elements(forest) - # Get the number of ghost elements of forest. - num_ghost_elements = t8_forest_get_num_ghosts(forest) - # Get the number of trees that have elements of this process. - num_local_trees = t8_forest_get_num_local_trees(forest) - - current_index = t8_locidx_t(0) - - local_num_conform = 0 - local_num_mortars = 0 - local_num_boundary = 0 - - for itree in 0:(num_local_trees - 1) - tree_class = t8_forest_get_tree_class(forest, itree) - eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) - - # Get the number of elements of this tree. - num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) - - for ielement in 0:(num_elements_in_tree - 1) - element = t8_forest_get_element_in_tree(forest, itree, ielement) - - level = t8_element_level(eclass_scheme, element) - - num_faces = t8_element_num_faces(eclass_scheme, element) - - for iface in 0:(num_faces - 1) - - # Compute the `orientation` of the touching faces. - if t8_element_is_root_boundary(eclass_scheme, element, iface) == 1 - cmesh = t8_forest_get_cmesh(forest) - itree_in_cmesh = t8_forest_ltreeid_to_cmesh_ltreeid(forest, itree) - iface_in_tree = t8_element_tree_face(eclass_scheme, element, iface) - orientation_ref = Ref{Cint}() - - t8_cmesh_get_face_neighbor(cmesh, itree_in_cmesh, iface_in_tree, C_NULL, - orientation_ref) - orientation = orientation_ref[] - else - orientation = zero(Cint) - end - - pelement_indices_ref = Ref{Ptr{t8_locidx_t}}() - pneighbor_leafs_ref = Ref{Ptr{Ptr{t8_element}}}() - pneigh_scheme_ref = Ref{Ptr{t8_eclass_scheme}}() - - dual_faces_ref = Ref{Ptr{Cint}}() - num_neighbors_ref = Ref{Cint}() - - forest_is_balanced = Cint(1) - - t8_forest_leaf_face_neighbors(forest, itree, element, - pneighbor_leafs_ref, iface, dual_faces_ref, - num_neighbors_ref, - pelement_indices_ref, pneigh_scheme_ref, - forest_is_balanced) - - num_neighbors = num_neighbors_ref[] - dual_faces = unsafe_wrap(Array, dual_faces_ref[], num_neighbors) - neighbor_ielements = unsafe_wrap(Array, pelement_indices_ref[], - num_neighbors) - neighbor_leafs = unsafe_wrap(Array, pneighbor_leafs_ref[], num_neighbors) - neighbor_scheme = pneigh_scheme_ref[] - - if num_neighbors > 0 - neighbor_level = t8_element_level(neighbor_scheme, neighbor_leafs[1]) - - # Conforming interface: The second condition ensures we only visit the interface once. - if level == neighbor_level && current_index <= neighbor_ielements[1] - local_num_conform += 1 - - faces = (iface, dual_faces[1]) - interface_id = local_num_conform - - # Write data to interfaces container. - interfaces.neighbor_ids[1, interface_id] = current_index + 1 - interfaces.neighbor_ids[2, interface_id] = neighbor_ielements[1] + 1 - - # Save interfaces.node_indices dimension specific in containers_3d.jl. - init_interface_node_indices!(interfaces, faces, orientation, - interface_id) - # Non-conforming interface. - elseif level < neighbor_level - local_num_mortars += 1 - - faces = (dual_faces[1], iface) - - mortar_id = local_num_mortars - - # Last entry is the large element. - mortars.neighbor_ids[end, mortar_id] = current_index + 1 - - # Fill in the `mortars.neighbor_ids` array and reorder if necessary. - init_mortar_neighbor_ids!(mortars, faces[2], faces[1], - orientation, neighbor_ielements, - mortar_id) - - # Fill in the `mortars.node_indices` array. - init_mortar_node_indices!(mortars, faces, orientation, mortar_id) - - # else: "level > neighbor_level" is skipped since we visit the mortar interface only once. - end - - # Domain boundary. - else - local_num_boundary += 1 - boundary_id = local_num_boundary - - boundaries.neighbor_ids[boundary_id] = current_index + 1 - - init_boundary_node_indices!(boundaries, iface, boundary_id) - - # One-based indexing. - boundaries.name[boundary_id] = boundary_names[iface + 1, itree + 1] - end - - t8_free(dual_faces_ref[]) - t8_free(pneighbor_leafs_ref[]) - t8_free(pelement_indices_ref[]) - end # for iface = ... - - current_index += 1 - end # for - end # for - - return (interfaces = local_num_conform, - mortars = local_num_mortars, - boundaries = local_num_boundary) -end - function trixi_t8_get_local_element_levels(forest) # Check that forest is a committed, that is valid and usable, forest. @assert t8_forest_is_committed(forest) != 0 @@ -341,23 +117,16 @@ function adapt_callback(forest, end function trixi_t8_adapt_new(old_forest, indicators) - # Check that forest is a committed, that is valid and usable, forest. - @assert t8_forest_is_committed(old_forest) != 0 - - # Init new forest. new_forest_ref = Ref{t8_forest_t}() t8_forest_init(new_forest_ref) new_forest = new_forest_ref[] - let set_from = C_NULL, recursive = 0, set_for_coarsening = 0, no_repartition = 0, - do_ghost = 1 - + let set_from = C_NULL, recursive = 0, no_repartition = 1, do_ghost = 1 t8_forest_set_user_data(new_forest, pointer(indicators)) t8_forest_set_adapt(new_forest, old_forest, @t8_adapt_callback(adapt_callback), recursive) t8_forest_set_balance(new_forest, set_from, no_repartition) - t8_forest_set_partition(new_forest, set_from, set_for_coarsening) - t8_forest_set_ghost(new_forest, do_ghost, T8_GHOST_FACES) # Note: MPI support not available yet so it is a dummy call. + t8_forest_set_ghost(new_forest, do_ghost, T8_GHOST_FACES) t8_forest_commit(new_forest) end diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 9f34a6b3b4b..4dbf44d29c4 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -97,6 +97,9 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi end print(f, ", " * string(variables[v]) * "_min") end + for variable in limiter.positivity_variables_nonlinear + print(f, ", " * string(variable) * "_min") + end end println(f) end @@ -142,6 +145,11 @@ end println(string(variables[v]) * ":\n- positivity: ", idp_bounds_delta_global[Symbol(string(v), "_min")]) end + for variable in limiter.positivity_variables_nonlinear + variable_string = string(variable) + println(variable_string * ":\n- positivity: ", + idp_bounds_delta_global[Symbol(variable_string, "_min")]) + end end println("─"^100 * "\n") diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 545d19b5136..19d73968c9a 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -60,6 +60,20 @@ deviation_threaded[stride_size * Threads.threadid()] = deviation end end + for variable in limiter.positivity_variables_nonlinear + key = Symbol(string(variable), "_min") + deviation_threaded = idp_bounds_delta_local[key] + @threaded for element in eachelement(solver, cache) + deviation = deviation_threaded[stride_size * Threads.threadid()] + for j in eachnode(solver), i in eachnode(solver) + var = variable(get_node_vars(u, equations, solver, i, j, element), + equations) + deviation = max(deviation, + variable_bounds[key][i, j, element] - var) + end + deviation_threaded[stride_size * Threads.threadid()] = deviation + end + end end for (key, _) in idp_bounds_delta_local @@ -92,6 +106,10 @@ print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")][stride_size]) end + for variable in limiter.positivity_variables_nonlinear + print(f, ", ", + idp_bounds_delta_local[Symbol(string(variable), "_min")][stride_size]) + end end println(f) end diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index 5854c8617c3..6f57d6647fc 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -726,7 +726,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::P4estMesh, return has_changed end -function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::SerialT8codeMesh, +function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::T8codeMesh, equations, dg::DG, cache, semi, t, iter; only_refine = false, only_coarsen = false, @@ -754,29 +754,29 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::SerialT8codeMe @trixi_timeit timer() "adapt" begin difference = @trixi_timeit timer() "mesh" trixi_t8_adapt!(mesh, indicators) - @trixi_timeit timer() "solver" adapt!(u_ode, adaptor, mesh, equations, dg, - cache, difference) - end + # Store whether there were any cells coarsened or refined and perform load balancing. + has_changed = any(difference .!= 0) - # Store whether there were any cells coarsened or refined and perform load balancing. - has_changed = any(difference .!= 0) + # Check if mesh changed on other processes + if mpi_isparallel() + has_changed = MPI.Allreduce!(Ref(has_changed), |, mpi_comm())[] + end - # TODO: T8codeMesh for MPI not implemented yet. - # Check if mesh changed on other processes - # if mpi_isparallel() - # has_changed = MPI.Allreduce!(Ref(has_changed), |, mpi_comm())[] - # end + if has_changed + @trixi_timeit timer() "solver" adapt!(u_ode, adaptor, mesh, equations, dg, + cache, difference) + end + end if has_changed - # TODO: T8codeMesh for MPI not implemented yet. - # if mpi_isparallel() && amr_callback.dynamic_load_balancing - # @trixi_timeit timer() "dynamic load balancing" begin - # global_first_quadrant = unsafe_wrap(Array, mesh.p4est.global_first_quadrant, mpi_nranks() + 1) - # old_global_first_quadrant = copy(global_first_quadrant) - # partition!(mesh) - # rebalance_solver!(u_ode, mesh, equations, dg, cache, old_global_first_quadrant) - # end - # end + if mpi_isparallel() && amr_callback.dynamic_load_balancing + @trixi_timeit timer() "dynamic load balancing" begin + old_global_first_element_ids = get_global_first_element_ids(mesh) + partition!(mesh) + rebalance_solver!(u_ode, mesh, equations, dg, cache, + old_global_first_element_ids) + end + end reinitialize_boundaries!(semi.boundary_conditions, cache) end diff --git a/src/callbacks_step/amr_dg.jl b/src/callbacks_step/amr_dg.jl index 1dcfdccdea8..0a7055af409 100644 --- a/src/callbacks_step/amr_dg.jl +++ b/src/callbacks_step/amr_dg.jl @@ -6,11 +6,14 @@ #! format: noindent # Redistribute data for load balancing after partitioning the mesh -function rebalance_solver!(u_ode::AbstractVector, mesh::ParallelP4estMesh, equations, +function rebalance_solver!(u_ode::AbstractVector, + mesh::Union{ParallelP4estMesh, ParallelT8codeMesh}, + equations, dg::DGSEM, cache, old_global_first_quadrant) - # mpi ranks are 0-based, this array uses 1-based indices - global_first_quadrant = unsafe_wrap(Array, mesh.p4est.global_first_quadrant, - mpi_nranks() + 1) + + # MPI ranks are 0-based. This array uses 1-based indices. + global_first_quadrant = get_global_first_element_ids(mesh) + if global_first_quadrant[mpi_rank() + 1] == old_global_first_quadrant[mpi_rank() + 1] && global_first_quadrant[mpi_rank() + 2] == diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index b816bc06e65..94524b23a3a 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -385,7 +385,12 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, # Return early if there is nothing to do. if !any(difference .!= 0) - return nothing + if mpi_isparallel() + # MPICache init uses all-to-all communication -> reinitialize even if there is nothing to do + # locally (there still might be other MPI ranks that have refined elements) + reinitialize_containers!(mesh, equations, dg, cache) + end + return end # Number of (local) cells/elements. diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 392cbba9e28..3f67951bafe 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -316,7 +316,12 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, # Return early if there is nothing to do. if !any(difference .!= 0) - return nothing + if mpi_isparallel() + # MPICache init uses all-to-all communication -> reinitialize even if there is nothing to do + # locally (there still might be other MPI ranks that have refined elements) + reinitialize_containers!(mesh, equations, dg, cache) + end + return end # Number of (local) cells/elements. diff --git a/src/callbacks_step/analysis_dg2d_parallel.jl b/src/callbacks_step/analysis_dg2d_parallel.jl index a04bf732604..000daa015dc 100644 --- a/src/callbacks_step/analysis_dg2d_parallel.jl +++ b/src/callbacks_step/analysis_dg2d_parallel.jl @@ -91,7 +91,8 @@ function calc_error_norms_per_element(func, u, t, analyzer, end function calc_error_norms(func, u, t, analyzer, - mesh::ParallelP4estMesh{2}, equations, + mesh::Union{ParallelP4estMesh{2}, ParallelT8codeMesh{2}}, + equations, initial_condition, dg::DGSEM, cache, cache_analysis) @unpack vandermonde, weights = analyzer @unpack node_coordinates, inverse_jacobian = cache.elements @@ -171,7 +172,8 @@ function integrate_via_indices(func::Func, u, end function integrate_via_indices(func::Func, u, - mesh::ParallelP4estMesh{2}, equations, + mesh::Union{ParallelP4estMesh{2}, ParallelT8codeMesh{2}}, + equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @unpack weights = dg.basis diff --git a/src/callbacks_step/analysis_dg3d_parallel.jl b/src/callbacks_step/analysis_dg3d_parallel.jl index d8756d91c9d..de777be406d 100644 --- a/src/callbacks_step/analysis_dg3d_parallel.jl +++ b/src/callbacks_step/analysis_dg3d_parallel.jl @@ -6,7 +6,8 @@ #! format: noindent function calc_error_norms(func, u, t, analyzer, - mesh::ParallelP4estMesh{3}, equations, + mesh::Union{ParallelP4estMesh{3}, ParallelT8codeMesh{3}}, + equations, initial_condition, dg::DGSEM, cache, cache_analysis) @unpack vandermonde, weights = analyzer @unpack node_coordinates, inverse_jacobian = cache.elements @@ -64,7 +65,8 @@ function calc_error_norms(func, u, t, analyzer, end function integrate_via_indices(func::Func, u, - mesh::ParallelP4estMesh{3}, equations, + mesh::Union{ParallelP4estMesh{3}, ParallelT8codeMesh{3}}, + equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @unpack weights = dg.basis diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 673c3ba6aa6..c6d32c0f6dc 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -174,4 +174,36 @@ function max_dt(u, t, mesh::ParallelP4estMesh{2}, return dt end + +function max_dt(u, t, mesh::ParallelT8codeMesh{2}, + constant_speed::False, equations, dg::DG, cache) + # call the method accepting a general `mesh::T8codeMesh{2}` + # TODO: MPI, we should improve this; maybe we should dispatch on `u` + # and create some MPI array type, overloading broadcasting and mapreduce etc. + # Then, this specific array type should also work well with DiffEq etc. + dt = invoke(max_dt, + Tuple{typeof(u), typeof(t), T8codeMesh{2}, + typeof(constant_speed), typeof(equations), typeof(dg), + typeof(cache)}, + u, t, mesh, constant_speed, equations, dg, cache) + dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[] + + return dt +end + +function max_dt(u, t, mesh::ParallelT8codeMesh{2}, + constant_speed::True, equations, dg::DG, cache) + # call the method accepting a general `mesh::T8codeMesh{2}` + # TODO: MPI, we should improve this; maybe we should dispatch on `u` + # and create some MPI array type, overloading broadcasting and mapreduce etc. + # Then, this specific array type should also work well with DiffEq etc. + dt = invoke(max_dt, + Tuple{typeof(u), typeof(t), T8codeMesh{2}, + typeof(constant_speed), typeof(equations), typeof(dg), + typeof(cache)}, + u, t, mesh, constant_speed, equations, dg, cache) + dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[] + + return dt +end end # @muladd diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl index 822ab2f87ec..664596f989e 100644 --- a/src/callbacks_step/stepsize_dg3d.jl +++ b/src/callbacks_step/stepsize_dg3d.jl @@ -150,4 +150,36 @@ function max_dt(u, t, mesh::ParallelP4estMesh{3}, return dt end + +function max_dt(u, t, mesh::ParallelT8codeMesh{3}, + constant_speed::False, equations, dg::DG, cache) + # call the method accepting a general `mesh::T8codeMesh{3}` + # TODO: MPI, we should improve this; maybe we should dispatch on `u` + # and create some MPI array type, overloading broadcasting and mapreduce etc. + # Then, this specific array type should also work well with DiffEq etc. + dt = invoke(max_dt, + Tuple{typeof(u), typeof(t), T8codeMesh{3}, + typeof(constant_speed), typeof(equations), typeof(dg), + typeof(cache)}, + u, t, mesh, constant_speed, equations, dg, cache) + dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[] + + return dt +end + +function max_dt(u, t, mesh::ParallelT8codeMesh{3}, + constant_speed::True, equations, dg::DG, cache) + # call the method accepting a general `mesh::T8codeMesh{3}` + # TODO: MPI, we should improve this; maybe we should dispatch on `u` + # and create some MPI array type, overloading broadcasting and mapreduce etc. + # Then, this specific array type should also work well with DiffEq etc. + dt = invoke(max_dt, + Tuple{typeof(u), typeof(t), T8codeMesh{3}, + typeof(constant_speed), typeof(equations), typeof(dg), + typeof(cache)}, + u, t, mesh, constant_speed, equations, dg, cache) + dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[] + + return dt +end end # @muladd diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 3c6f759db2b..f5a632723cf 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1632,6 +1632,18 @@ end return p end +# Transformation from conservative variables u to d(p)/d(u) +@inline function gradient_conservative(::typeof(pressure), + u, equations::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + + return (equations.gamma - 1.0) * SVector(0.5 * v_square, -v1, -v2, 1.0) +end + @inline function density_pressure(u, equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u rho_times_p = (equations.gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2)) @@ -1699,4 +1711,13 @@ end @inline function energy_internal(cons, equations::CompressibleEulerEquations2D) return energy_total(cons, equations) - energy_kinetic(cons, equations) end + +# State validation for Newton-bisection method of subcell IDP limiting +@inline function Base.isvalid(u, equations::CompressibleEulerEquations2D) + p = pressure(u, equations) + if u[1] <= 0.0 || p <= 0.0 + return false + end + return true +end end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index bfd051c0106..d11ead8a2b8 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -376,6 +376,12 @@ of the correct length `nvariables(equations)`. """ function energy_internal end +# Default implementation of gradient for `variable`. Used for subcell limiting. +# Implementing a gradient function for a specific variable improves the performance. +@inline function gradient_conservative(variable, u, equations) + return ForwardDiff.gradient(x -> variable(x, equations), u) +end + #################################################################################################### # Include files with actual implementations for different systems of equations. diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 43d1991e34b..4366cd32f08 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -1118,6 +1118,20 @@ end return p end +# Transformation from conservative variables u to d(p)/d(u) +@inline function gradient_conservative(::typeof(pressure), + u, equations::IdealGlmMhdEquations2D) + rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + v_square = v1^2 + v2^2 + v3^2 + + return (equations.gamma - 1.0) * + SVector(0.5 * v_square, -v1, -v2, -v3, 1.0, -B1, -B2, -B3, -psi) +end + @inline function density_pressure(u, equations::IdealGlmMhdEquations2D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho @@ -1384,6 +1398,15 @@ end cons[9]^2 / 2) end +# State validation for Newton-bisection method of subcell IDP limiting +@inline function Base.isvalid(u, equations::IdealGlmMhdEquations2D) + p = pressure(u, equations) + if u[1] <= 0.0 || p <= 0.0 + return false + end + return true +end + # Calculate the cross helicity (\vec{v}⋅\vec{B}) for a conservative state `cons' @inline function cross_helicity(cons, ::IdealGlmMhdEquations2D) return (cons[2] * cons[6] + cons[3] * cons[7] + cons[4] * cons[8]) / cons[1] diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index c5d39ef00c0..abe9d9345b5 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -1700,6 +1700,10 @@ function bilinear_interpolation!(coordinate, face_vertices, u, v) end end +function get_global_first_element_ids(mesh::P4estMesh) + return unsafe_wrap(Array, mesh.p4est.global_first_quadrant, mpi_nranks() + 1) +end + function balance!(mesh::P4estMesh{2}, init_fn = C_NULL) p4est_balance(mesh.p4est, P4EST_CONNECT_FACE, init_fn) # Due to a bug in `p4est`, the forest needs to be rebalanced twice sometimes diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index c9665a22af9..6fb4d861d10 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -7,8 +7,6 @@ to manage trees and mesh refinement. """ mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <: AbstractMesh{NDIMS} - cmesh :: Ptr{t8_cmesh} # cpointer to coarse mesh - scheme :: Ptr{t8_eclass_scheme} # cpointer to element scheme forest :: Ptr{t8_forest} # cpointer to forest is_parallel :: IsParallel @@ -25,14 +23,15 @@ mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <: nmortars :: Int nboundaries :: Int - function T8codeMesh{NDIMS}(cmesh, scheme, forest, tree_node_coordinates, nodes, + nmpiinterfaces :: Int + nmpimortars :: Int + + function T8codeMesh{NDIMS}(forest, tree_node_coordinates, nodes, boundary_names, current_filename) where {NDIMS} - is_parallel = False() + is_parallel = mpi_isparallel() ? True() : False() - mesh = new{NDIMS, Float64, typeof(is_parallel), NDIMS + 2, length(nodes)}(cmesh, - scheme, - forest, + mesh = new{NDIMS, Float64, typeof(is_parallel), NDIMS + 2, length(nodes)}(forest, is_parallel) mesh.nodes = nodes @@ -52,7 +51,7 @@ mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <: # further down. However, this might cause a pile-up of `mesh` # objects during long-running sessions. if !MPI.Finalized() - trixi_t8_unref_forest(mesh.forest) + t8_forest_unref(Ref(mesh.forest)) end end @@ -63,7 +62,7 @@ mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <: # more information. if haskey(ENV, "TRIXI_T8CODE_SC_FINALIZE") MPI.add_finalize_hook!() do - trixi_t8_unref_forest(mesh.forest) + t8_forest_unref(Ref(mesh.forest)) end end @@ -72,16 +71,15 @@ mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <: end const SerialT8codeMesh{NDIMS} = T8codeMesh{NDIMS, <:Real, <:False} +const ParallelT8codeMesh{NDIMS} = T8codeMesh{NDIMS, <:Real, <:True} @inline mpi_parallel(mesh::SerialT8codeMesh) = False() +@inline mpi_parallel(mesh::ParallelT8codeMesh) = True() @inline Base.ndims(::T8codeMesh{NDIMS}) where {NDIMS} = NDIMS @inline Base.real(::T8codeMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT -@inline ntrees(mesh::T8codeMesh) = Int(t8_forest_get_num_local_trees(mesh.forest)) +@inline ntrees(mesh::T8codeMesh) = size(mesh.tree_node_coordinates)[end] @inline ncells(mesh::T8codeMesh) = Int(t8_forest_get_local_num_elements(mesh.forest)) -@inline ninterfaces(mesh::T8codeMesh) = mesh.ninterfaces -@inline nmortars(mesh::T8codeMesh) = mesh.nmortars -@inline nboundaries(mesh::T8codeMesh) = mesh.nboundaries function Base.show(io::IO, mesh::T8codeMesh) print(io, "T8codeMesh{", ndims(mesh), ", ", real(mesh), "}") @@ -184,21 +182,23 @@ function T8codeMesh(trees_per_dimension; polydeg = 1, T8code.Libt8.p8est_connectivity_destroy(conn) end + do_face_ghost = mpi_isparallel() scheme = t8_scheme_new_default_cxx() - forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, 0, mpi_comm()) + forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, do_face_ghost, + mpi_comm()) basis = LobattoLegendreBasis(RealT, polydeg) nodes = basis.nodes + num_trees = t8_cmesh_get_num_trees(cmesh) + tree_node_coordinates = Array{RealT, NDIMS + 2}(undef, NDIMS, ntuple(_ -> length(nodes), NDIMS)..., - prod(trees_per_dimension)) + num_trees) # Get cell length in reference mesh: Omega_ref = [-1,1]^NDIMS. dx = [2 / n for n in trees_per_dimension] - num_local_trees = t8_cmesh_get_num_local_trees(cmesh) - # Non-periodic boundaries. boundary_names = fill(Symbol("---"), 2 * NDIMS, prod(trees_per_dimension)) @@ -208,7 +208,7 @@ function T8codeMesh(trees_per_dimension; polydeg = 1, mapping_ = mapping end - for itree in 1:num_local_trees + for itree in 1:num_trees veptr = t8_cmesh_get_tree_vertices(cmesh, itree - 1) verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS)) @@ -256,7 +256,7 @@ function T8codeMesh(trees_per_dimension; polydeg = 1, end end - return T8codeMesh{NDIMS}(cmesh, scheme, forest, tree_node_coordinates, nodes, + return T8codeMesh{NDIMS}(forest, tree_node_coordinates, nodes, boundary_names, "") end @@ -290,21 +290,25 @@ function T8codeMesh(cmesh::Ptr{t8_cmesh}; @assert (NDIMS == 2||NDIMS == 3) "NDIMS should be 2 or 3." + do_face_ghost = mpi_isparallel() scheme = t8_scheme_new_default_cxx() - forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, 0, mpi_comm()) + forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, do_face_ghost, + mpi_comm()) basis = LobattoLegendreBasis(RealT, polydeg) nodes = basis.nodes - num_local_trees = t8_cmesh_get_num_local_trees(cmesh) + num_trees = t8_cmesh_get_num_trees(cmesh) tree_node_coordinates = Array{RealT, NDIMS + 2}(undef, NDIMS, ntuple(_ -> length(nodes), NDIMS)..., - num_local_trees) + num_trees) nodes_in = [-1.0, 1.0] matrix = polynomial_interpolation_matrix(nodes_in, nodes) + num_local_trees = t8_cmesh_get_num_local_trees(cmesh) + if NDIMS == 2 data_in = Array{RealT, 3}(undef, 2, 2, 2) tmp1 = zeros(RealT, 2, length(nodes), length(nodes_in)) @@ -353,7 +357,7 @@ function T8codeMesh(cmesh::Ptr{t8_cmesh}; tmp1 = zeros(RealT, 3, length(nodes), length(nodes_in), length(nodes_in)) verts = zeros(3, 8) - for itree in 0:(num_local_trees - 1) + for itree in 0:(num_trees - 1) veptr = t8_cmesh_get_tree_vertices(cmesh, itree) # Note, `verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS))` @@ -387,9 +391,9 @@ function T8codeMesh(cmesh::Ptr{t8_cmesh}; map_node_coordinates!(tree_node_coordinates, mapping) # There's no simple and generic way to distinguish boundaries. Name all of them :all. - boundary_names = fill(:all, 2 * NDIMS, num_local_trees) + boundary_names = fill(:all, 2 * NDIMS, num_trees) - return T8codeMesh{NDIMS}(cmesh, scheme, forest, tree_node_coordinates, nodes, + return T8codeMesh{NDIMS}(forest, tree_node_coordinates, nodes, boundary_names, "") end @@ -442,7 +446,7 @@ function T8codeMesh(conn::Ptr{p8est_connectivity}; kwargs...) end """ - T8codeMesh{NDIMS}(meshfile::String; kwargs...) + T8codeMesh(meshfile::String, ndims; kwargs...) Main mesh constructor for the `T8codeMesh` that imports an unstructured, conforming mesh from a Gmsh mesh file (`.msh`). @@ -461,7 +465,6 @@ mesh from a Gmsh mesh file (`.msh`). - `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. """ function T8codeMesh(meshfile::String, ndims; kwargs...) - # Prevent `t8code` from crashing Julia if the file doesn't exist. @assert isfile(meshfile) @@ -586,13 +589,525 @@ function adapt!(mesh::T8codeMesh, adapt_callback; recursive = true, balance = tr return nothing end -# TODO: Just a placeholder. Will be implemented later when MPI is supported. -function balance!(mesh::T8codeMesh, init_fn = C_NULL) +""" + Trixi.balance!(mesh::T8codeMesh) + +Balance a `T8codeMesh` to ensure 2^(NDIMS-1):1 face neighbors. +""" +function balance!(mesh::T8codeMesh) + new_forest_ref = Ref{t8_forest_t}() + t8_forest_init(new_forest_ref) + new_forest = new_forest_ref[] + + let set_from = mesh.forest, no_repartition = 1, do_ghost = 1 + t8_forest_set_balance(new_forest, set_from, no_repartition) + t8_forest_set_ghost(new_forest, do_ghost, T8_GHOST_FACES) + t8_forest_commit(new_forest) + end + + mesh.forest = new_forest + + return nothing +end + +""" + Trixi.partition!(mesh::T8codeMesh) + +Partition a `T8codeMesh` in order to redistribute elements evenly among MPI ranks. + +# Arguments +- `mesh::T8codeMesh`: Initialized mesh object. +""" +function partition!(mesh::T8codeMesh) + new_forest_ref = Ref{t8_forest_t}() + t8_forest_init(new_forest_ref) + new_forest = new_forest_ref[] + + let set_from = mesh.forest, do_ghost = 1, allow_for_coarsening = 1 + t8_forest_set_partition(new_forest, set_from, allow_for_coarsening) + t8_forest_set_ghost(new_forest, do_ghost, T8_GHOST_FACES) + t8_forest_commit(new_forest) + end + + mesh.forest = new_forest + return nothing end -# TODO: Just a placeholder. Will be implemented later when MPI is supported. -function partition!(mesh::T8codeMesh; allow_coarsening = true, weight_fn = C_NULL) +# Compute the global ids (zero-indexed) of first element in each MPI rank. +function get_global_first_element_ids(mesh::T8codeMesh) + n_elements_local = Int(t8_forest_get_local_num_elements(mesh.forest)) + n_elements_by_rank = Vector{Int}(undef, mpi_nranks()) + n_elements_by_rank[mpi_rank() + 1] = n_elements_local + MPI.Allgather!(MPI.UBuffer(n_elements_by_rank, 1), mpi_comm()) + return [sum(n_elements_by_rank[1:(rank - 1)]) for rank in 1:(mpi_nranks() + 1)] +end + +function count_interfaces(mesh::T8codeMesh) + @assert t8_forest_is_committed(mesh.forest) != 0 + + num_local_elements = t8_forest_get_local_num_elements(mesh.forest) + num_local_trees = t8_forest_get_num_local_trees(mesh.forest) + + current_index = t8_locidx_t(0) + + local_num_conform = 0 + local_num_mortars = 0 + local_num_boundary = 0 + + local_num_mpi_conform = 0 + local_num_mpi_mortars = 0 + + visited_global_mortar_ids = Set{UInt64}([]) + + max_level = t8_forest_get_maxlevel(mesh.forest) #UInt64 + max_tree_num_elements = UInt64(2^ndims(mesh))^max_level + + if mpi_isparallel() + ghost_num_trees = t8_forest_ghost_num_trees(mesh.forest) + + ghost_tree_element_offsets = [num_local_elements + + t8_forest_ghost_get_tree_element_offset(mesh.forest, + itree) + for itree in 0:(ghost_num_trees - 1)] + ghost_global_treeids = [t8_forest_ghost_get_global_treeid(mesh.forest, itree) + for itree in 0:(ghost_num_trees - 1)] + end + + for itree in 0:(num_local_trees - 1) + tree_class = t8_forest_get_tree_class(mesh.forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(mesh.forest, tree_class) + + num_elements_in_tree = t8_forest_get_tree_num_elements(mesh.forest, itree) + + global_itree = t8_forest_global_tree_id(mesh.forest, itree) + + for ielement in 0:(num_elements_in_tree - 1) + element = t8_forest_get_element_in_tree(mesh.forest, itree, ielement) + + level = t8_element_level(eclass_scheme, element) + + num_faces = t8_element_num_faces(eclass_scheme, element) + + # Note: This works only for forests of one element class. + current_linear_id = global_itree * max_tree_num_elements + + t8_element_get_linear_id(eclass_scheme, element, max_level) + + for iface in 0:(num_faces - 1) + pelement_indices_ref = Ref{Ptr{t8_locidx_t}}() + pneighbor_leafs_ref = Ref{Ptr{Ptr{t8_element}}}() + pneigh_scheme_ref = Ref{Ptr{t8_eclass_scheme}}() + + dual_faces_ref = Ref{Ptr{Cint}}() + num_neighbors_ref = Ref{Cint}() + + forest_is_balanced = Cint(1) + + t8_forest_leaf_face_neighbors(mesh.forest, itree, element, + pneighbor_leafs_ref, iface, dual_faces_ref, + num_neighbors_ref, + pelement_indices_ref, pneigh_scheme_ref, + forest_is_balanced) + + num_neighbors = num_neighbors_ref[] + dual_faces = unsafe_wrap(Array, dual_faces_ref[], num_neighbors) + neighbor_ielements = unsafe_wrap(Array, pelement_indices_ref[], + num_neighbors) + neighbor_leafs = unsafe_wrap(Array, pneighbor_leafs_ref[], num_neighbors) + neighbor_scheme = pneigh_scheme_ref[] + + if num_neighbors == 0 + local_num_boundary += 1 + else + neighbor_level = t8_element_level(neighbor_scheme, neighbor_leafs[1]) + + if all(neighbor_ielements .< num_local_elements) + # Conforming interface: The second condition ensures we + # only visit the interface once. + if level == neighbor_level && current_index <= neighbor_ielements[1] + local_num_conform += 1 + elseif level < neighbor_level + local_num_mortars += 1 + # `else level > neighbor_level` is ignored since we + # only want to count the mortar interface once. + end + else + if level == neighbor_level + local_num_mpi_conform += 1 + elseif level < neighbor_level + local_num_mpi_mortars += 1 + + global_mortar_id = 2 * ndims(mesh) * current_linear_id + iface + + else # level > neighbor_level + neighbor_global_ghost_itree = ghost_global_treeids[findlast(ghost_tree_element_offsets .<= + neighbor_ielements[1])] + neighbor_linear_id = neighbor_global_ghost_itree * + max_tree_num_elements + + t8_element_get_linear_id(neighbor_scheme, + neighbor_leafs[1], + max_level) + global_mortar_id = 2 * ndims(mesh) * neighbor_linear_id + + dual_faces[1] + + if !(global_mortar_id in visited_global_mortar_ids) + push!(visited_global_mortar_ids, global_mortar_id) + local_num_mpi_mortars += 1 + end + end + end + end + + t8_free(dual_faces_ref[]) + t8_free(pneighbor_leafs_ref[]) + t8_free(pelement_indices_ref[]) + end # for + + current_index += 1 + end # for + end # for + + return (interfaces = local_num_conform, + mortars = local_num_mortars, + boundaries = local_num_boundary, + mpi_interfaces = local_num_mpi_conform, + mpi_mortars = local_num_mpi_mortars) +end + +# I know this routine is an unmaintainable behemoth. However, I see no real +# and elegant way to refactor this into, for example, smaller parts. The +# `t8_forest_leaf_face_neighbors` routine is as of now rather costly and it +# makes sense to query it only once per face per element and extract all the +# information needed at once in order to fill the connectivity information. +# Instead, I opted for good documentation. +function fill_mesh_info!(mesh::T8codeMesh, interfaces, mortars, boundaries, + boundary_names; mpi_mesh_info = nothing) + @assert t8_forest_is_committed(mesh.forest) != 0 + + num_local_elements = t8_forest_get_local_num_elements(mesh.forest) + num_local_trees = t8_forest_get_num_local_trees(mesh.forest) + + if !isnothing(mpi_mesh_info) + #! format: off + remotes = t8_forest_ghost_get_remotes(mesh.forest) + ghost_num_trees = t8_forest_ghost_num_trees(mesh.forest) + + ghost_remote_first_elem = [num_local_elements + + t8_forest_ghost_remote_first_elem(mesh.forest, remote) + for remote in remotes] + + ghost_tree_element_offsets = [num_local_elements + + t8_forest_ghost_get_tree_element_offset(mesh.forest, itree) + for itree in 0:(ghost_num_trees - 1)] + + ghost_global_treeids = [t8_forest_ghost_get_global_treeid(mesh.forest, itree) + for itree in 0:(ghost_num_trees - 1)] + #! format: on + end + + # Process-local index of the current element in the space-filling curve. + current_index = t8_locidx_t(0) + + # Increment counters for the different interface/mortar/boundary types. + local_num_conform = 0 + local_num_mortars = 0 + local_num_boundary = 0 + + local_num_mpi_conform = 0 + local_num_mpi_mortars = 0 + + # Works for quads and hexs only. This mapping is needed in the MPI mortar + # sections below. + map_iface_to_ichild_to_position = [ + # 0 1 2 3 4 5 6 7 ichild/iface + [1, 0, 2, 0, 3, 0, 4, 0], # 0 + [0, 1, 0, 2, 0, 3, 0, 4], # 1 + [1, 2, 0, 0, 3, 4, 0, 0], # 2 + [0, 0, 1, 2, 0, 0, 3, 4], # 3 + [1, 2, 3, 4, 0, 0, 0, 0], # 4 + [0, 0, 0, 0, 1, 2, 3, 4], # 5 + ] + + # Helper variables to compute unique global MPI interface/mortar ids. + max_level = t8_forest_get_maxlevel(mesh.forest) #UInt64 + max_tree_num_elements = UInt64(2^ndims(mesh))^max_level + + # These two variables help to ensure that we count MPI mortars from smaller + # elements point of view only once. + visited_global_mortar_ids = Set{UInt64}([]) + global_mortar_id_to_local = Dict{UInt64, Int}([]) + + # Loop over all local trees. + for itree in 0:(num_local_trees - 1) + tree_class = t8_forest_get_tree_class(mesh.forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(mesh.forest, tree_class) + + num_elements_in_tree = t8_forest_get_tree_num_elements(mesh.forest, itree) + + global_itree = t8_forest_global_tree_id(mesh.forest, itree) + + # Loop over all local elements of the current local tree. + for ielement in 0:(num_elements_in_tree - 1) + element = t8_forest_get_element_in_tree(mesh.forest, itree, ielement) + + level = t8_element_level(eclass_scheme, element) + + num_faces = t8_element_num_faces(eclass_scheme, element) + + # Note: This works only for forests of one element class. + current_linear_id = global_itree * max_tree_num_elements + + t8_element_get_linear_id(eclass_scheme, element, max_level) + + # Loop over all faces of the current local element. + for iface in 0:(num_faces - 1) + # Compute the `orientation` of the touching faces. + if t8_element_is_root_boundary(eclass_scheme, element, iface) == 1 + cmesh = t8_forest_get_cmesh(mesh.forest) + itree_in_cmesh = t8_forest_ltreeid_to_cmesh_ltreeid(mesh.forest, itree) + iface_in_tree = t8_element_tree_face(eclass_scheme, element, iface) + orientation_ref = Ref{Cint}() + + t8_cmesh_get_face_neighbor(cmesh, itree_in_cmesh, iface_in_tree, C_NULL, + orientation_ref) + orientation = orientation_ref[] + else + orientation = zero(Cint) + end + + pelement_indices_ref = Ref{Ptr{t8_locidx_t}}() + pneighbor_leafs_ref = Ref{Ptr{Ptr{t8_element}}}() + pneigh_scheme_ref = Ref{Ptr{t8_eclass_scheme}}() + + dual_faces_ref = Ref{Ptr{Cint}}() + num_neighbors_ref = Ref{Cint}() + + forest_is_balanced = Cint(1) + + # Query neighbor information from t8code. + t8_forest_leaf_face_neighbors(mesh.forest, itree, element, + pneighbor_leafs_ref, iface, dual_faces_ref, + num_neighbors_ref, + pelement_indices_ref, pneigh_scheme_ref, + forest_is_balanced) + + num_neighbors = num_neighbors_ref[] + dual_faces = unsafe_wrap(Array, dual_faces_ref[], num_neighbors) + neighbor_ielements = unsafe_wrap(Array, pelement_indices_ref[], + num_neighbors) + neighbor_leafs = unsafe_wrap(Array, pneighbor_leafs_ref[], num_neighbors) + neighbor_scheme = pneigh_scheme_ref[] + + # Now we check for the different cases. The nested if-structure is as follows: + # + # if `boundary`: + # + # + # else: // It must be an interface or mortar. + # + # if `all neighbors are local elements`: + # + # if `local interface`: + # + # elseif `local mortar from larger element point of view`: + # + # else: // `local mortar from smaller elements point of view` + # // We only count local mortars once. + # + # else: // It must be either a MPI interface or a MPI mortar. + # + # if `MPI interface`: + # + # elseif `MPI mortar from larger element point of view`: + # + # else: // `MPI mortar from smaller elements point of view` + # + # + # // end + + # Domain boundary. + if num_neighbors == 0 + local_num_boundary += 1 + boundary_id = local_num_boundary + + boundaries.neighbor_ids[boundary_id] = current_index + 1 + + init_boundary_node_indices!(boundaries, iface, boundary_id) + + # One-based indexing. + boundaries.name[boundary_id] = boundary_names[iface + 1, itree + 1] + + # Interface or mortar. + else + neighbor_level = t8_element_level(neighbor_scheme, neighbor_leafs[1]) + + # Local interface or mortar. + if all(neighbor_ielements .< num_local_elements) + + # Local interface: The second condition ensures we only visit the interface once. + if level == neighbor_level && current_index <= neighbor_ielements[1] + local_num_conform += 1 + + interfaces.neighbor_ids[1, local_num_conform] = current_index + + 1 + interfaces.neighbor_ids[2, local_num_conform] = neighbor_ielements[1] + + 1 + + init_interface_node_indices!(interfaces, (iface, dual_faces[1]), + orientation, + local_num_conform) + # Local mortar. + elseif level < neighbor_level + local_num_mortars += 1 + + # Last entry is the large element. + mortars.neighbor_ids[end, local_num_mortars] = current_index + 1 + + init_mortar_neighbor_ids!(mortars, iface, dual_faces[1], + orientation, neighbor_ielements, + local_num_mortars) + + init_mortar_node_indices!(mortars, (dual_faces[1], iface), + orientation, local_num_mortars) + + # else: `level > neighbor_level` is skipped since we visit the mortar interface only once. + end + + # MPI interface or MPI mortar. + else + + # MPI interface. + if level == neighbor_level + local_num_mpi_conform += 1 + + neighbor_global_ghost_itree = ghost_global_treeids[findlast(ghost_tree_element_offsets .<= + neighbor_ielements[1])] + + neighbor_linear_id = neighbor_global_ghost_itree * + max_tree_num_elements + + t8_element_get_linear_id(neighbor_scheme, + neighbor_leafs[1], + max_level) + + if current_linear_id < neighbor_linear_id + local_side = 1 + smaller_iface = iface + smaller_linear_id = current_linear_id + faces = (iface, dual_faces[1]) + else + local_side = 2 + smaller_iface = dual_faces[1] + smaller_linear_id = neighbor_linear_id + faces = (dual_faces[1], iface) + end + + global_interface_id = 2 * ndims(mesh) * smaller_linear_id + + smaller_iface + + mpi_mesh_info.mpi_interfaces.local_neighbor_ids[local_num_mpi_conform] = current_index + + 1 + mpi_mesh_info.mpi_interfaces.local_sides[local_num_mpi_conform] = local_side + + init_mpi_interface_node_indices!(mpi_mesh_info.mpi_interfaces, + faces, local_side, orientation, + local_num_mpi_conform) + + neighbor_rank = remotes[findlast(ghost_remote_first_elem .<= + neighbor_ielements[1])] + mpi_mesh_info.neighbor_ranks_interface[local_num_mpi_conform] = neighbor_rank + + mpi_mesh_info.global_interface_ids[local_num_mpi_conform] = global_interface_id + + # MPI Mortar: from larger element point of view + elseif level < neighbor_level + local_num_mpi_mortars += 1 + + global_mortar_id = 2 * ndims(mesh) * current_linear_id + iface + + neighbor_ids = neighbor_ielements .+ 1 + + local_neighbor_positions = findall(neighbor_ids .<= + num_local_elements) + local_neighbor_ids = [neighbor_ids[i] + for i in local_neighbor_positions] + local_neighbor_positions = [map_iface_to_ichild_to_position[dual_faces[1] + 1][t8_element_child_id(neighbor_scheme, neighbor_leafs[i]) + 1] + for i in local_neighbor_positions] + + # Last entry is the large element. + push!(local_neighbor_ids, current_index + 1) + push!(local_neighbor_positions, 2^(ndims(mesh) - 1) + 1) + + mpi_mesh_info.mpi_mortars.local_neighbor_ids[local_num_mpi_mortars] = local_neighbor_ids + mpi_mesh_info.mpi_mortars.local_neighbor_positions[local_num_mpi_mortars] = local_neighbor_positions + + init_mortar_node_indices!(mpi_mesh_info.mpi_mortars, + (dual_faces[1], iface), orientation, + local_num_mpi_mortars) + + neighbor_ranks = [remotes[findlast(ghost_remote_first_elem .<= + ineighbor_ghost)] + for ineighbor_ghost in filter(x -> x >= + num_local_elements, + neighbor_ielements)] + mpi_mesh_info.neighbor_ranks_mortar[local_num_mpi_mortars] = neighbor_ranks + + mpi_mesh_info.global_mortar_ids[local_num_mpi_mortars] = global_mortar_id + + # MPI Mortar: from smaller elements point of view + else + neighbor_global_ghost_itree = ghost_global_treeids[findlast(ghost_tree_element_offsets .<= + neighbor_ielements[1])] + neighbor_linear_id = neighbor_global_ghost_itree * + max_tree_num_elements + + t8_element_get_linear_id(neighbor_scheme, + neighbor_leafs[1], + max_level) + global_mortar_id = 2 * ndims(mesh) * neighbor_linear_id + + dual_faces[1] + + if global_mortar_id in visited_global_mortar_ids + local_mpi_mortar_id = global_mortar_id_to_local[global_mortar_id] + + push!(mpi_mesh_info.mpi_mortars.local_neighbor_ids[local_mpi_mortar_id], + current_index + 1) + push!(mpi_mesh_info.mpi_mortars.local_neighbor_positions[local_mpi_mortar_id], + map_iface_to_ichild_to_position[iface + 1][t8_element_child_id(eclass_scheme, element) + 1]) + else + local_num_mpi_mortars += 1 + local_mpi_mortar_id = local_num_mpi_mortars + push!(visited_global_mortar_ids, global_mortar_id) + global_mortar_id_to_local[global_mortar_id] = local_mpi_mortar_id + + mpi_mesh_info.mpi_mortars.local_neighbor_ids[local_mpi_mortar_id] = [ + current_index + 1, + ] + mpi_mesh_info.mpi_mortars.local_neighbor_positions[local_mpi_mortar_id] = [ + map_iface_to_ichild_to_position[iface + 1][t8_element_child_id(eclass_scheme, element) + 1], + ] + init_mortar_node_indices!(mpi_mesh_info.mpi_mortars, + (iface, dual_faces[1]), + orientation, local_mpi_mortar_id) + + neighbor_ranks = [ + remotes[findlast(ghost_remote_first_elem .<= + neighbor_ielements[1])], + ] + mpi_mesh_info.neighbor_ranks_mortar[local_mpi_mortar_id] = neighbor_ranks + + mpi_mesh_info.global_mortar_ids[local_mpi_mortar_id] = global_mortar_id + end + end + end + end + + t8_free(dual_faces_ref[]) + t8_free(pneighbor_leafs_ref[]) + t8_free(pelement_indices_ref[]) + end # for iface + + current_index += 1 + end # for ielement + end # for itree + return nothing end diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 0941ae6a8ca..dc21dbe9a1e 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -1,3 +1,10 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + """ SemidiscretizationCoupled @@ -65,11 +72,13 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationCoupled) summary_line(io, "system", i) mesh, equations, solver, _ = mesh_equations_solver_cache(semi.semis[i]) summary_line(increment_indent(io), "mesh", mesh |> typeof |> nameof) - summary_line(increment_indent(io), "equations", equations |> typeof |> nameof) + summary_line(increment_indent(io), "equations", + equations |> typeof |> nameof) summary_line(increment_indent(io), "initial condition", semi.semis[i].initial_condition) # no boundary conditions since that could be too much - summary_line(increment_indent(io), "source terms", semi.semis[i].source_terms) + summary_line(increment_indent(io), "source terms", + semi.semis[i].source_terms) summary_line(increment_indent(io), "solver", solver |> typeof |> nameof) end summary_line(io, "total #DOFs per field", ndofs(semi)) @@ -106,20 +115,14 @@ end @inline Base.real(semi::SemidiscretizationCoupled) = promote_type(real.(semi.semis)...) -@inline Base.eltype(semi::SemidiscretizationCoupled) = promote_type(eltype.(semi.semis)...) +@inline function Base.eltype(semi::SemidiscretizationCoupled) + promote_type(eltype.(semi.semis)...) +end @inline function ndofs(semi::SemidiscretizationCoupled) sum(ndofs, semi.semis) end -@inline function nelements(semi::SemidiscretizationCoupled) - return sum(semi.semis) do semi_ - mesh, equations, solver, cache = mesh_equations_solver_cache(semi_) - - nelements(mesh, solver, cache) - end -end - function compute_coefficients(t, semi::SemidiscretizationCoupled) @unpack u_indices = semi @@ -137,23 +140,40 @@ end @view u_ode[semi.u_indices[index]] end +# Same as `foreach(enumerate(something))`, but without allocations. +# +# Note that compile times may increase if this is used with big tuples. +@inline foreach_enumerate(func, collection) = foreach_enumerate(func, collection, 1) +@inline foreach_enumerate(func, collection::Tuple{}, index) = nothing + +@inline function foreach_enumerate(func, collection, index) + element = first(collection) + remaining_collection = Base.tail(collection) + + func((index, element)) + + # Process remaining collection + foreach_enumerate(func, remaining_collection, index + 1) +end + function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) @unpack u_indices = semi time_start = time_ns() @trixi_timeit timer() "copy to coupled boundaries" begin - for semi_ in semi.semis - copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi) + foreach(semi.semis) do semi_ + copy_to_coupled_boundary!(semi_.boundary_conditions, u_ode, semi, semi_) end end # Call rhs! for each semidiscretization - for i in eachsystem(semi) - u_loc = get_system_u_ode(u_ode, i, semi) - du_loc = get_system_u_ode(du_ode, i, semi) - - @trixi_timeit timer() "system #$i" rhs!(du_loc, u_loc, semi.semis[i], t) + @trixi_timeit timer() "copy to coupled boundaries" begin + foreach_enumerate(semi.semis) do (i, semi_) + u_loc = get_system_u_ode(u_ode, i, semi) + du_loc = get_system_u_ode(du_ode, i, semi) + rhs!(du_loc, u_loc, semi_, t) + end end runtime = time_ns() - time_start @@ -309,7 +329,8 @@ end for i in eachsystem(semi) u_ode_slice = get_system_u_ode(u_ode, i, semi) - save_solution_file(semis[i], u_ode_slice, solution_callback, integrator, system = i) + save_solution_file(semis[i], u_ode_slice, solution_callback, integrator, + system = i) end end @@ -332,7 +353,7 @@ end ################################################################################ """ - BoundaryConditionCoupled(other_semi_index, indices, uEltype) + BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter) Boundary condition to glue two meshes together. Solution values at the boundary of another mesh will be used as boundary values. This requires the use @@ -348,32 +369,37 @@ This is currently only implemented for [`StructuredMesh`](@ref). - `indices::Tuple`: node/cell indices at the boundary of the mesh in the other semidiscretization. See examples below. - `uEltype::Type`: element type of solution +- `coupling_converter::CouplingConverter`: function to call for converting the solution + state of one system to the other system # Examples ```julia # Connect the left boundary of mesh 2 to our boundary such that our positive # boundary direction will match the positive y direction of the other boundary -BoundaryConditionCoupled(2, (:begin, :i), Float64) +BoundaryConditionCoupled(2, (:begin, :i), Float64, fun) # Connect the same two boundaries oppositely oriented -BoundaryConditionCoupled(2, (:begin, :i_backwards), Float64) +BoundaryConditionCoupled(2, (:begin, :i_backwards), Float64, fun) # Using this as y_neg boundary will connect `our_cells[i, 1, j]` to `other_cells[j, end-i, end]` -BoundaryConditionCoupled(2, (:j, :i_backwards, :end), Float64) +BoundaryConditionCoupled(2, (:j, :i_backwards, :end), Float64, fun) ``` !!! warning "Experimental code" This is an experimental feature and can change any time. """ -mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indices} +mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indices, + CouplingConverter} # NDIMST2M1 == NDIMS * 2 - 1 # Buffer for boundary values: [variable, nodes_i, nodes_j, cell_i, cell_j] - u_boundary :: Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1 - other_semi_index :: Int - other_orientation :: Int - indices :: Indices - - function BoundaryConditionCoupled(other_semi_index, indices, uEltype) + u_boundary :: Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1 + other_semi_index :: Int + other_orientation :: Int + indices :: Indices + coupling_converter :: CouplingConverter + + function BoundaryConditionCoupled(other_semi_index, indices, uEltype, + coupling_converter) NDIMS = length(indices) u_boundary = Array{uEltype, NDIMS * 2 - 1}(undef, ntuple(_ -> 0, NDIMS * 2 - 1)) @@ -385,8 +411,10 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic other_orientation = 3 end - new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices)}(u_boundary, other_semi_index, - other_orientation, indices) + new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices), + typeof(coupling_converter)}(u_boundary, + other_semi_index, other_orientation, + indices, coupling_converter) end end @@ -395,8 +423,10 @@ function Base.eltype(boundary_condition::BoundaryConditionCoupled) end function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, direction, - cell_indices, surface_node_indices, - surface_flux_function, equations) + cell_indices, + surface_node_indices, + surface_flux_function, + equations) # get_node_vars(boundary_condition.u_boundary, equations, solver, surface_node_indices..., cell_indices...), # but we don't have a solver here u_boundary = SVector(ntuple(v -> boundary_condition.u_boundary[v, @@ -421,13 +451,15 @@ function allocate_coupled_boundary_conditions(semi::AbstractSemidiscretization) for direction in 1:n_boundaries boundary_condition = semi.boundary_conditions[direction] - allocate_coupled_boundary_condition(boundary_condition, direction, mesh, equations, + allocate_coupled_boundary_condition(boundary_condition, direction, mesh, + equations, solver) end end # Don't do anything for other BCs than BoundaryConditionCoupled -function allocate_coupled_boundary_condition(boundary_condition, direction, mesh, equations, +function allocate_coupled_boundary_condition(boundary_condition, direction, mesh, + equations, solver) return nothing end @@ -448,43 +480,69 @@ function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditi end # Don't do anything for other BCs than BoundaryConditionCoupled -function copy_to_coupled_boundary!(boundary_condition, u_ode, semi) +function copy_to_coupled_boundary!(boundary_condition, u_ode, semi_coupled, semi) return nothing end +function copy_to_coupled_boundary!(u_ode, semi_coupled, semi, i, n_boundaries, + boundary_condition, boundary_conditions...) + copy_to_coupled_boundary!(boundary_condition, u_ode, semi_coupled, semi) + if i < n_boundaries + copy_to_coupled_boundary!(u_ode, semi_coupled, semi, i + 1, n_boundaries, + boundary_conditions...) + end +end + function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple}, u_ode, - semi) - for boundary_condition in boundary_conditions - copy_to_coupled_boundary!(boundary_condition, u_ode, semi) + semi_coupled, semi) + copy_to_coupled_boundary!(u_ode, semi_coupled, semi, 1, length(boundary_conditions), + boundary_conditions...) +end + +function mesh_equations_solver_cache(other_semi_index, i, semi_, semi_tuple...) + if i == other_semi_index + return mesh_equations_solver_cache(semi_) + else + # Walk through semidiscretizations until we find `i` + mesh_equations_solver_cache(other_semi_index, i + 1, semi_tuple...) end end # In 2D -function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, u_ode, - semi) - @unpack u_indices = semi +function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, + u_ode, + semi_coupled, semi) + @unpack u_indices = semi_coupled @unpack other_semi_index, other_orientation, indices = boundary_condition + @unpack coupling_converter, u_boundary = boundary_condition + + mesh_own, equations_own, solver_own, cache_own = mesh_equations_solver_cache(semi) + + mesh_other, equations_other, solver_other, cache_other = mesh_equations_solver_cache(other_semi_index, + 1, + semi_coupled.semis...) - mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semis[other_semi_index]) - u = wrap_array(get_system_u_ode(u_ode, other_semi_index, semi), mesh, equations, solver, - cache) + node_coordinates_other = cache_other.elements.node_coordinates + u_ode_other = get_system_u_ode(u_ode, other_semi_index, semi_coupled) + u_other = wrap_array(u_ode_other, mesh_other, equations_other, solver_other, + cache_other) - linear_indices = LinearIndices(size(mesh)) + linear_indices = LinearIndices(size(mesh_other)) if other_orientation == 1 - cells = axes(mesh, 2) + cells = axes(mesh_other, 2) else # other_orientation == 2 - cells = axes(mesh, 1) + cells = axes(mesh_other, 1) end # Copy solution data to the coupled boundary using "delayed indexing" with # a start value and a step size to get the correct face and orientation. - node_index_range = eachnode(solver) + node_index_range = eachnode(solver_other) i_node_start, i_node_step = index_to_start_step_2d(indices[1], node_index_range) j_node_start, j_node_step = index_to_start_step_2d(indices[2], node_index_range) - i_cell_start, i_cell_step = index_to_start_step_2d(indices[1], axes(mesh, 1)) - j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], axes(mesh, 2)) + i_cell_start, i_cell_step = index_to_start_step_2d(indices[1], axes(mesh_other, 1)) + j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], axes(mesh_other, 2)) i_cell = i_cell_start j_cell = j_cell_start @@ -492,16 +550,26 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ for cell in cells i_node = i_node_start j_node = j_node_start - - for i in eachnode(solver) - for v in 1:size(u, 1) - boundary_condition.u_boundary[v, i, cell] = u[v, i_node, j_node, - linear_indices[i_cell, - j_cell]] + element_id = linear_indices[i_cell, j_cell] + + for element_id in eachnode(solver_other) + x_other = get_node_coords(node_coordinates_other, equations_other, + solver_other, + i_node, j_node, linear_indices[i_cell, j_cell]) + u_node_other = get_node_vars(u_other, equations_other, solver_other, i_node, + j_node, linear_indices[i_cell, j_cell]) + u_node_converted = coupling_converter(x_other, u_node_other, + equations_other, + equations_own) + + for i in eachindex(u_node_converted) + u_boundary[i, element_id, cell] = u_node_converted[i] end + i_node += i_node_step j_node += j_node_step end + i_cell += i_cell_step j_cell += j_cell_step end @@ -511,7 +579,8 @@ end ### DGSEM/structured ################################################################################ -@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, +@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, + orientation, boundary_condition::BoundaryConditionCoupled, mesh::StructuredMesh, equations, surface_integral, dg::DG, cache, @@ -531,7 +600,8 @@ end sign_jacobian = sign(inverse_jacobian[node_indices..., element]) # Contravariant vector Ja^i is the normal vector - normal = sign_jacobian * get_contravariant_vector(orientation, contravariant_vectors, + normal = sign_jacobian * + get_contravariant_vector(orientation, contravariant_vectors, node_indices..., element) # If the mapping is orientation-reversing, the normal vector will be reversed (see above). @@ -608,3 +678,4 @@ function analyze_convergence(errors_coupled, iterations, return eoc_mean_values end +end # @muladd diff --git a/src/solvers/dgsem_p4est/containers_parallel.jl b/src/solvers/dgsem_p4est/containers_parallel.jl index 7c7bd868457..fd2749155bb 100644 --- a/src/solvers/dgsem_p4est/containers_parallel.jl +++ b/src/solvers/dgsem_p4est/containers_parallel.jl @@ -43,7 +43,8 @@ function Base.resize!(mpi_interfaces::P4estMPIInterfaceContainer, capacity) end # Create MPI interface container and initialize interface data -function init_mpi_interfaces(mesh::ParallelP4estMesh, equations, basis, elements) +function init_mpi_interfaces(mesh::Union{ParallelP4estMesh, ParallelT8codeMesh}, + equations, basis, elements) NDIMS = ndims(elements) uEltype = eltype(elements) @@ -133,7 +134,8 @@ function Base.resize!(mpi_mortars::P4estMPIMortarContainer, capacity) end # Create MPI mortar container and initialize MPI mortar data -function init_mpi_mortars(mesh::ParallelP4estMesh, equations, basis, elements) +function init_mpi_mortars(mesh::Union{ParallelP4estMesh, ParallelT8codeMesh}, equations, + basis, elements) NDIMS = ndims(mesh) RealT = real(mesh) uEltype = eltype(elements) diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index a8887351c46..3bf0cd0cab5 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -6,7 +6,8 @@ #! format: noindent function prolong2mpiinterfaces!(cache, u, - mesh::ParallelP4estMesh{2}, + mesh::Union{ParallelP4estMesh{2}, + ParallelT8codeMesh{2}}, equations, surface_integral, dg::DG) @unpack mpi_interfaces = cache index_range = eachnode(dg) @@ -43,7 +44,8 @@ function prolong2mpiinterfaces!(cache, u, end function calc_mpi_interface_flux!(surface_flux_values, - mesh::ParallelP4estMesh{2}, + mesh::Union{ParallelP4estMesh{2}, + ParallelT8codeMesh{2}}, nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @@ -106,7 +108,8 @@ end # Inlined version of the interface flux computation for conservation laws @inline function calc_mpi_interface_flux!(surface_flux_values, - mesh::P4estMesh{2}, + mesh::Union{ParallelP4estMesh{2}, + ParallelT8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -131,7 +134,8 @@ end end function prolong2mpimortars!(cache, u, - mesh::ParallelP4estMesh{2}, equations, + mesh::Union{ParallelP4estMesh{2}, ParallelT8codeMesh{2}}, + equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DGSEM) @unpack node_indices = cache.mpi_mortars @@ -199,7 +203,7 @@ function prolong2mpimortars!(cache, u, end function calc_mpi_mortar_flux!(surface_flux_values, - mesh::ParallelP4estMesh{2}, + mesh::Union{ParallelP4estMesh{2}, ParallelT8codeMesh{2}}, nonconservative_terms, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @@ -253,7 +257,8 @@ end # Inlined version of the mortar flux computation on small elements for conservation laws @inline function calc_mpi_mortar_flux!(fstar, - mesh::ParallelP4estMesh{2}, + mesh::Union{ParallelP4estMesh{2}, + ParallelT8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, mortar_index, position_index, normal_direction, @@ -271,7 +276,9 @@ end end @inline function mpi_mortar_fluxes_to_elements!(surface_flux_values, - mesh::ParallelP4estMesh{2}, equations, + mesh::Union{ParallelP4estMesh{2}, + ParallelT8codeMesh{2}}, + equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM, cache, mortar, fstar, u_buffer) diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel.jl b/src/solvers/dgsem_p4est/dg_3d_parallel.jl index 13bf2a1a2eb..e504e06d2c4 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parallel.jl @@ -6,7 +6,7 @@ #! format: noindent function rhs!(du, u, t, - mesh::ParallelP4estMesh{3}, equations, + mesh::Union{ParallelP4estMesh{3}, ParallelT8codeMesh{3}}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Start to receive MPI data @@ -113,7 +113,8 @@ function rhs!(du, u, t, end function prolong2mpiinterfaces!(cache, u, - mesh::ParallelP4estMesh{3}, + mesh::Union{ParallelP4estMesh{3}, + ParallelT8codeMesh{3}}, equations, surface_integral, dg::DG) @unpack mpi_interfaces = cache index_range = eachnode(dg) @@ -160,7 +161,8 @@ function prolong2mpiinterfaces!(cache, u, end function calc_mpi_interface_flux!(surface_flux_values, - mesh::ParallelP4estMesh{3}, + mesh::Union{ParallelP4estMesh{3}, + ParallelT8codeMesh{3}}, nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack local_neighbor_ids, node_indices, local_sides = cache.mpi_interfaces @@ -237,7 +239,8 @@ end # Inlined version of the interface flux computation for conservation laws @inline function calc_mpi_interface_flux!(surface_flux_values, - mesh::P4estMesh{3}, + mesh::Union{ParallelP4estMesh{3}, + ParallelT8codeMesh{3}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -265,7 +268,8 @@ end end function prolong2mpimortars!(cache, u, - mesh::ParallelP4estMesh{3}, equations, + mesh::Union{ParallelP4estMesh{3}, ParallelT8codeMesh{3}}, + equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DGSEM) @unpack node_indices = cache.mpi_mortars @@ -374,7 +378,7 @@ function prolong2mpimortars!(cache, u, end function calc_mpi_mortar_flux!(surface_flux_values, - mesh::ParallelP4estMesh{3}, + mesh::Union{ParallelP4estMesh{3}, ParallelT8codeMesh{3}}, nonconservative_terms, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @@ -437,7 +441,8 @@ end # Inlined version of the mortar flux computation on small elements for conservation laws @inline function calc_mpi_mortar_flux!(fstar, - mesh::ParallelP4estMesh{3}, + mesh::Union{ParallelP4estMesh{3}, + ParallelT8codeMesh{3}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, mortar_index, position_index, normal_direction, @@ -456,7 +461,9 @@ end end @inline function mpi_mortar_fluxes_to_elements!(surface_flux_values, - mesh::ParallelP4estMesh{3}, equations, + mesh::Union{ParallelP4estMesh{3}, + ParallelT8codeMesh{3}}, + equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM, cache, mortar, fstar, u_buffer, fstar_tmp) diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl index 712ede2bfce..eaa6ab5cee2 100644 --- a/src/solvers/dgsem_p4est/dg_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_parallel.jl @@ -166,7 +166,8 @@ end # at `index_base`+1 in the MPI buffer. `data_size` is the data size associated with each small # position (i.e. position 1 or 2). The data corresponding to the large side (i.e. position 3) has # size `2 * data_size`. -@inline function buffer_mortar_indices(mesh::ParallelP4estMesh{2}, index_base, +@inline function buffer_mortar_indices(mesh::Union{ParallelP4estMesh{2}, + ParallelT8codeMesh{2}}, index_base, data_size) return ( # first, last for local element in position 1 (small element) @@ -185,7 +186,8 @@ end # at `index_base`+1 in the MPI buffer. `data_size` is the data size associated with each small # position (i.e. position 1 to 4). The data corresponding to the large side (i.e. position 5) has # size `4 * data_size`. -@inline function buffer_mortar_indices(mesh::ParallelP4estMesh{3}, index_base, +@inline function buffer_mortar_indices(mesh::Union{ParallelP4estMesh{3}, + ParallelT8codeMesh{3}}, index_base, data_size) return ( # first, last for local element in position 1 (small element) @@ -491,7 +493,8 @@ end # Exchange normal directions of small elements of the MPI mortars. They are needed on all involved # MPI ranks to calculate the mortar fluxes. -function exchange_normal_directions!(mpi_mortars, mpi_cache, mesh::ParallelP4estMesh, +function exchange_normal_directions!(mpi_mortars, mpi_cache, + mesh::Union{ParallelP4estMesh, ParallelT8codeMesh}, n_nodes) RealT = real(mesh) n_dims = ndims(mesh) diff --git a/src/solvers/dgsem_t8code/containers.jl b/src/solvers/dgsem_t8code/containers.jl index 093feb2985a..d7ff79fbf2f 100644 --- a/src/solvers/dgsem_t8code/containers.jl +++ b/src/solvers/dgsem_t8code/containers.jl @@ -18,19 +18,22 @@ function reinitialize_containers!(mesh::T8codeMesh, equations, dg::DGSEM, cache) @unpack boundaries = cache resize!(boundaries, mesh.nboundaries) - trixi_t8_fill_mesh_info(mesh.forest, elements, interfaces, mortars, boundaries, - mesh.boundary_names) + fill_mesh_info!(mesh, interfaces, mortars, boundaries, + mesh.boundary_names) return nothing end function count_required_surfaces!(mesh::T8codeMesh) - counts = trixi_t8_count_interfaces(mesh.forest) + counts = count_interfaces(mesh) mesh.nmortars = counts.mortars mesh.ninterfaces = counts.interfaces mesh.nboundaries = counts.boundaries + mesh.nmpimortars = counts.mpi_mortars + mesh.nmpiinterfaces = counts.mpi_interfaces + return counts end @@ -38,7 +41,9 @@ end function count_required_surfaces(mesh::T8codeMesh) return (interfaces = mesh.ninterfaces, mortars = mesh.nmortars, - boundaries = mesh.nboundaries) + boundaries = mesh.nboundaries, + mpi_interfaces = mesh.nmpiinterfaces, + mpi_mortars = mesh.nmpimortars) end # Compatibility to `dgsem_p4est/containers.jl`. diff --git a/src/solvers/dgsem_t8code/containers_2d.jl b/src/solvers/dgsem_t8code/containers_2d.jl index bf77826a34b..ce525bfdf65 100644 --- a/src/solvers/dgsem_t8code/containers_2d.jl +++ b/src/solvers/dgsem_t8code/containers_2d.jl @@ -26,6 +26,7 @@ function calc_node_coordinates!(node_coordinates, tree_class = t8_forest_get_tree_class(mesh.forest, itree) eclass_scheme = t8_forest_get_eclass_scheme(mesh.forest, tree_class) num_elements_in_tree = t8_forest_get_tree_num_elements(mesh.forest, itree) + global_itree = t8_forest_global_tree_id(mesh.forest, itree) for ielement in 0:(num_elements_in_tree - 1) element = t8_forest_get_element_in_tree(mesh.forest, itree, ielement) @@ -55,7 +56,7 @@ function calc_node_coordinates!(node_coordinates, multiply_dimensionwise!(view(node_coordinates, :, :, :, current_index += 1), matrix1, matrix2, view(mesh.tree_node_coordinates, :, :, :, - itree + 1), + global_itree + 1), tmp1) end end diff --git a/src/solvers/dgsem_t8code/containers_3d.jl b/src/solvers/dgsem_t8code/containers_3d.jl index f2d54ff07da..4d56bc734aa 100644 --- a/src/solvers/dgsem_t8code/containers_3d.jl +++ b/src/solvers/dgsem_t8code/containers_3d.jl @@ -28,6 +28,7 @@ function calc_node_coordinates!(node_coordinates, tree_class = t8_forest_get_tree_class(mesh.forest, itree) eclass_scheme = t8_forest_get_eclass_scheme(mesh.forest, tree_class) num_elements_in_tree = t8_forest_get_tree_num_elements(mesh.forest, itree) + global_itree = t8_forest_global_tree_id(mesh.forest, itree) for ielement in 0:(num_elements_in_tree - 1) element = t8_forest_get_element_in_tree(mesh.forest, itree, ielement) @@ -63,7 +64,7 @@ function calc_node_coordinates!(node_coordinates, current_index += 1), matrix1, matrix2, matrix3, view(mesh.tree_node_coordinates, :, :, :, :, - itree + 1), + global_itree + 1), tmp1) end end diff --git a/src/solvers/dgsem_t8code/containers_parallel.jl b/src/solvers/dgsem_t8code/containers_parallel.jl new file mode 100644 index 00000000000..0cb3f5887a0 --- /dev/null +++ b/src/solvers/dgsem_t8code/containers_parallel.jl @@ -0,0 +1,65 @@ +function reinitialize_containers!(mesh::ParallelT8codeMesh, equations, dg::DGSEM, cache) + @unpack elements, interfaces, boundaries, mortars, mpi_interfaces, mpi_mortars, + mpi_cache = cache + resize!(elements, ncells(mesh)) + init_elements!(elements, mesh, dg.basis) + + count_required_surfaces!(mesh) + required = count_required_surfaces(mesh) + + resize!(interfaces, required.interfaces) + + resize!(boundaries, required.boundaries) + + resize!(mortars, required.mortars) + + resize!(mpi_interfaces, required.mpi_interfaces) + + resize!(mpi_mortars, required.mpi_mortars) + + mpi_mesh_info = (mpi_mortars = mpi_mortars, + mpi_interfaces = mpi_interfaces, + + # Temporary arrays for updating `mpi_cache`. + global_mortar_ids = fill(UInt64(0), nmpimortars(mpi_mortars)), + global_interface_ids = fill(UInt64(0), nmpiinterfaces(mpi_interfaces)), + neighbor_ranks_mortar = Vector{Vector{Int}}(undef, + nmpimortars(mpi_mortars)), + neighbor_ranks_interface = fill(-1, nmpiinterfaces(mpi_interfaces))) + + fill_mesh_info!(mesh, interfaces, mortars, boundaries, + mesh.boundary_names; mpi_mesh_info = mpi_mesh_info) + + init_mpi_cache!(mpi_cache, mesh, mpi_mesh_info, nvariables(equations), nnodes(dg), + eltype(elements)) + + empty!(mpi_mesh_info.global_mortar_ids) + empty!(mpi_mesh_info.global_interface_ids) + empty!(mpi_mesh_info.neighbor_ranks_mortar) + empty!(mpi_mesh_info.neighbor_ranks_interface) + + # Re-initialize and distribute normal directions of MPI mortars; requires + # MPI communication, so the MPI cache must be re-initialized beforehand. + init_normal_directions!(mpi_mortars, dg.basis, elements) + exchange_normal_directions!(mpi_mortars, mpi_cache, mesh, nnodes(dg)) + + return nothing +end + +# Compatibility to `dgsem_p4est/containers.jl`. +function init_mpi_interfaces!(interfaces, mesh::ParallelT8codeMesh) + # Do nothing. + return nothing +end + +# Compatibility to `dgsem_p4est/containers.jl`. +function init_mpi_mortars!(mortars, mesh::ParallelT8codeMesh) + # Do nothing. + return nothing +end + +# Compatibility to `dgsem_p4est/containers_parallel.jl`. +function init_mpi_mortars!(mpi_mortars, mesh::ParallelT8codeMesh, basis, elements) + # Do nothing. + return nothing +end diff --git a/src/solvers/dgsem_t8code/dg.jl b/src/solvers/dgsem_t8code/dg.jl index 6e9660c917d..e01b12e0f80 100644 --- a/src/solvers/dgsem_t8code/dg.jl +++ b/src/solvers/dgsem_t8code/dg.jl @@ -13,8 +13,8 @@ function create_cache(mesh::T8codeMesh, equations::AbstractEquations, dg::DG, :: boundaries = init_boundaries(mesh, equations, dg.basis, elements) mortars = init_mortars(mesh, equations, dg.basis, elements) - trixi_t8_fill_mesh_info(mesh.forest, elements, interfaces, mortars, boundaries, - mesh.boundary_names) + fill_mesh_info!(mesh, interfaces, mortars, boundaries, + mesh.boundary_names) cache = (; elements, interfaces, boundaries, mortars) @@ -29,4 +29,7 @@ end include("containers.jl") include("containers_2d.jl") include("containers_3d.jl") + +include("containers_parallel.jl") +include("dg_parallel.jl") end # @muladd diff --git a/src/solvers/dgsem_t8code/dg_parallel.jl b/src/solvers/dgsem_t8code/dg_parallel.jl new file mode 100644 index 00000000000..ece614b7d75 --- /dev/null +++ b/src/solvers/dgsem_t8code/dg_parallel.jl @@ -0,0 +1,135 @@ +@muladd begin +#! format: noindent + +# This method is called when a `SemidiscretizationHyperbolic` is constructed. +# It constructs the basic `cache` used throughout the simulation to compute +# the RHS etc. +function create_cache(mesh::ParallelT8codeMesh, equations::AbstractEquations, dg::DG, + ::Any, + ::Type{uEltype}) where {uEltype <: Real} + # Make sure to balance and partition the forest before creating any + # containers in case someone has tampered with forest after creating the + # mesh. + balance!(mesh) + partition!(mesh) + + count_required_surfaces!(mesh) + + elements = init_elements(mesh, equations, dg.basis, uEltype) + mortars = init_mortars(mesh, equations, dg.basis, elements) + interfaces = init_interfaces(mesh, equations, dg.basis, elements) + boundaries = init_boundaries(mesh, equations, dg.basis, elements) + + mpi_mortars = init_mpi_mortars(mesh, equations, dg.basis, elements) + mpi_interfaces = init_mpi_interfaces(mesh, equations, dg.basis, elements) + + mpi_mesh_info = (mpi_mortars = mpi_mortars, + mpi_interfaces = mpi_interfaces, + global_mortar_ids = fill(UInt64(0), nmpimortars(mpi_mortars)), + global_interface_ids = fill(UInt64(0), + nmpiinterfaces(mpi_interfaces)), + neighbor_ranks_mortar = Vector{Vector{Int}}(undef, + nmpimortars(mpi_mortars)), + neighbor_ranks_interface = fill(-1, + nmpiinterfaces(mpi_interfaces))) + + fill_mesh_info!(mesh, interfaces, mortars, boundaries, + mesh.boundary_names; mpi_mesh_info = mpi_mesh_info) + + mpi_cache = init_mpi_cache(mesh, mpi_mesh_info, nvariables(equations), nnodes(dg), + uEltype) + + empty!(mpi_mesh_info.global_mortar_ids) + empty!(mpi_mesh_info.global_interface_ids) + empty!(mpi_mesh_info.neighbor_ranks_mortar) + empty!(mpi_mesh_info.neighbor_ranks_interface) + + init_normal_directions!(mpi_mortars, dg.basis, elements) + exchange_normal_directions!(mpi_mortars, mpi_cache, mesh, nnodes(dg)) + + cache = (; elements, interfaces, mpi_interfaces, boundaries, mortars, mpi_mortars, + mpi_cache) + + # Add specialized parts of the cache required to compute the volume integral etc. + cache = (; cache..., + create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...) + cache = (; cache..., create_cache(mesh, equations, dg.mortar, uEltype)...) + + return cache +end + +function init_mpi_cache(mesh::ParallelT8codeMesh, mpi_mesh_info, nvars, nnodes, uEltype) + mpi_cache = P4estMPICache(uEltype) + init_mpi_cache!(mpi_cache, mesh, mpi_mesh_info, nvars, nnodes, uEltype) + return mpi_cache +end + +function init_mpi_cache!(mpi_cache::P4estMPICache, mesh::ParallelT8codeMesh, + mpi_mesh_info, nvars, nnodes, uEltype) + mpi_neighbor_ranks, mpi_neighbor_interfaces, mpi_neighbor_mortars = init_mpi_neighbor_connectivity(mpi_mesh_info, + mesh) + + mpi_send_buffers, mpi_recv_buffers, mpi_send_requests, mpi_recv_requests = init_mpi_data_structures(mpi_neighbor_interfaces, + mpi_neighbor_mortars, + ndims(mesh), + nvars, + nnodes, + uEltype) + + n_elements_global = Int(t8_forest_get_global_num_elements(mesh.forest)) + n_elements_local = Int(t8_forest_get_local_num_elements(mesh.forest)) + + n_elements_by_rank = Vector{Int}(undef, mpi_nranks()) + n_elements_by_rank[mpi_rank() + 1] = n_elements_local + + MPI.Allgather!(MPI.UBuffer(n_elements_by_rank, 1), mpi_comm()) + + n_elements_by_rank = OffsetArray(n_elements_by_rank, 0:(mpi_nranks() - 1)) + + # Account for 1-based indexing in Julia. + first_element_global_id = sum(n_elements_by_rank[0:(mpi_rank() - 1)]) + 1 + + @assert n_elements_global==sum(n_elements_by_rank) "error in total number of elements" + + @pack! mpi_cache = mpi_neighbor_ranks, mpi_neighbor_interfaces, + mpi_neighbor_mortars, + mpi_send_buffers, mpi_recv_buffers, + mpi_send_requests, mpi_recv_requests, + n_elements_by_rank, n_elements_global, + first_element_global_id + + return mpi_cache +end + +function init_mpi_neighbor_connectivity(mpi_mesh_info, mesh::ParallelT8codeMesh) + @unpack mpi_interfaces, mpi_mortars, global_interface_ids, neighbor_ranks_interface, global_mortar_ids, neighbor_ranks_mortar = mpi_mesh_info + + mpi_neighbor_ranks = vcat(neighbor_ranks_interface, neighbor_ranks_mortar...) |> + sort |> unique + + p = sortperm(global_interface_ids) + + neighbor_ranks_interface .= neighbor_ranks_interface[p] + interface_ids = collect(1:nmpiinterfaces(mpi_interfaces))[p] + + p = sortperm(global_mortar_ids) + neighbor_ranks_mortar .= neighbor_ranks_mortar[p] + mortar_ids = collect(1:nmpimortars(mpi_mortars))[p] + + # For each neighbor rank, init connectivity data structures + mpi_neighbor_interfaces = Vector{Vector{Int}}(undef, length(mpi_neighbor_ranks)) + mpi_neighbor_mortars = Vector{Vector{Int}}(undef, length(mpi_neighbor_ranks)) + for (index, d) in enumerate(mpi_neighbor_ranks) + mpi_neighbor_interfaces[index] = interface_ids[findall(==(d), + neighbor_ranks_interface)] + mpi_neighbor_mortars[index] = mortar_ids[findall(x -> (d in x), + neighbor_ranks_mortar)] + end + + # Check that all interfaces were counted exactly once + @assert mapreduce(length, +, mpi_neighbor_interfaces; init = 0) == + nmpiinterfaces(mpi_interfaces) + + return mpi_neighbor_ranks, mpi_neighbor_interfaces, mpi_neighbor_mortars +end +end # @muladd diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index 8095dae123a..157d462aa2f 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl @@ -446,7 +446,8 @@ function init_mpi_neighbor_connectivity(elements, mpi_interfaces, mpi_mortars, end function rhs!(du, u, t, - mesh::Union{ParallelTreeMesh{2}, ParallelP4estMesh{2}}, equations, + mesh::Union{ParallelTreeMesh{2}, ParallelP4estMesh{2}, + ParallelT8codeMesh{2}}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Start to receive MPI data diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 055e7ce24a4..e433c953779 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -16,18 +16,28 @@ end SubcellLimiterIDP(equations::AbstractEquations, basis; local_minmax_variables_cons = String[], positivity_variables_cons = String[], - positivity_correction_factor = 0.1) + positivity_variables_nonlinear = [], + positivity_correction_factor = 0.1, + max_iterations_newton = 10, + newton_tolerances = (1.0e-12, 1.0e-14), + gamma_constant_newton = 2 * ndims(equations)) Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref) including: - Local maximum/minimum Zalesak-type limiting for conservative variables (`local_minmax_variables_cons`) -- Positivity limiting for conservative variables (`positivity_variables_cons`) +- Positivity limiting for conservative variables (`positivity_variables_cons`) and nonlinear variables +(`positivity_variables_nonlinear`) Conservative variables to be limited are passed as a vector of strings, e.g. `local_minmax_variables_cons = ["rho"]` -and `positivity_variables_cons = ["rho"]`. +and `positivity_variables_cons = ["rho"]`. For nonlinear variables the specific functions are +passed in a vector, e.g. `positivity_variables_nonlinear = [pressure]`. The bounds are calculated using the low-order FV solution. The positivity limiter uses `positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`. +The limiting of nonlinear variables uses a Newton-bisection method with a maximum of +`max_iterations_newton` iterations, relative and absolute tolerances of `newton_tolerances` +and a provisional update constant `gamma_constant_newton` (`gamma_constant_newton>=2*d`, +where `d = #dimensions`). See equation (20) of Pazner (2020) and equation (30) of Rueda-Ramírez et al. (2022). !!! note This limiter and the correction callback [`SubcellLimiterIDPCorrection`](@ref) only work together. @@ -45,22 +55,32 @@ The bounds are calculated using the low-order FV solution. The positivity limite !!! warning "Experimental implementation" This is an experimental feature and may change in future releases. """ -struct SubcellLimiterIDP{RealT <: Real, Cache} <: AbstractSubcellLimiter +struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear, Cache} <: + AbstractSubcellLimiter local_minmax::Bool local_minmax_variables_cons::Vector{Int} # Local mininum/maximum principles for conservative variables positivity::Bool positivity_variables_cons::Vector{Int} # Positivity for conservative variables + positivity_variables_nonlinear::LimitingVariablesNonlinear # Positivity for nonlinear variables positivity_correction_factor::RealT cache::Cache + max_iterations_newton::Int + newton_tolerances::Tuple{RealT, RealT} # Relative and absolute tolerances for Newton's method + gamma_constant_newton::RealT # Constant for the subcell limiting of convex (nonlinear) constraints end # this method is used when the limiter is constructed as for shock-capturing volume integrals function SubcellLimiterIDP(equations::AbstractEquations, basis; local_minmax_variables_cons = String[], positivity_variables_cons = String[], - positivity_correction_factor = 0.1) + positivity_variables_nonlinear = [], + positivity_correction_factor = 0.1, + max_iterations_newton = 10, + newton_tolerances = (1.0e-12, 1.0e-14), + gamma_constant_newton = 2 * ndims(equations)) local_minmax = (length(local_minmax_variables_cons) > 0) - positivity = (length(positivity_variables_cons) > 0) + positivity = (length(positivity_variables_cons) + + length(positivity_variables_nonlinear) > 0) local_minmax_variables_cons_ = get_variable_index.(local_minmax_variables_cons, equations) @@ -80,13 +100,20 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis; bound_keys = (bound_keys..., Symbol(string(v), "_min")) end end + for variable in positivity_variables_nonlinear + bound_keys = (bound_keys..., Symbol(string(variable), "_min")) + end cache = create_cache(SubcellLimiterIDP, equations, basis, bound_keys) SubcellLimiterIDP{typeof(positivity_correction_factor), + typeof(positivity_variables_nonlinear), typeof(cache)}(local_minmax, local_minmax_variables_cons_, positivity, positivity_variables_cons_, - positivity_correction_factor, cache) + positivity_variables_nonlinear, + positivity_correction_factor, cache, + max_iterations_newton, newton_tolerances, + gamma_constant_newton) end function Base.show(io::IO, limiter::SubcellLimiterIDP) @@ -97,10 +124,15 @@ function Base.show(io::IO, limiter::SubcellLimiterIDP) if !(local_minmax || positivity) print(io, "No limiter selected => pure DG method") else - print(io, "limiter=(") - local_minmax && print(io, "min/max limiting, ") - positivity && print(io, "positivity") - print(io, "), ") + features = String[] + if local_minmax + push!(features, "local min/max") + end + if positivity + push!(features, "positivity") + end + join(io, features, ", ") + print(io, "Limiter=($features), ") end print(io, "Local bounds with FV solution") print(io, ")") @@ -120,15 +152,15 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) if local_minmax setup = [ setup..., - "" => "local maximum/minimum bounds for conservative variables $(limiter.local_minmax_variables_cons)", + "" => "Local maximum/minimum limiting for conservative variables $(limiter.local_minmax_variables_cons)", ] end if positivity - string = "positivity for conservative variables $(limiter.positivity_variables_cons)" + string = "Positivity limiting for conservative variables $(limiter.positivity_variables_cons) and $(limiter.positivity_variables_nonlinear)" setup = [setup..., "" => string] setup = [ setup..., - "" => " positivity correction factor = $(limiter.positivity_correction_factor)", + "" => "- with positivity correction factor = $(limiter.positivity_correction_factor)", ] end setup = [ diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 3d272359fe4..3f7954c8958 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -5,6 +5,10 @@ @muladd begin #! format: noindent +############################################################################### +# IDP Limiting +############################################################################### + # this method is used when the limiter is constructed as for shock-capturing volume integrals function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquations{2}, basis::LobattoLegendreBasis, bound_keys) @@ -66,6 +70,9 @@ function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSE return nothing end +############################################################################### +# Calculation of local bounds using low-order FV solution + @inline function calc_bounds_twosided!(var_min, var_max, variable, u, t, semi) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) # Calc bounds inside elements @@ -164,6 +171,9 @@ end return nothing end +############################################################################### +# Local minimum/maximum limiting + @inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi) for variable in limiter.local_minmax_variables_cons idp_local_minmax!(alpha, limiter, u, t, dt, semi, variable) @@ -233,16 +243,36 @@ end return nothing end +############################################################################### +# Global positivity limiting + @inline function idp_positivity!(alpha, limiter, u, dt, semi) # Conservative variables for variable in limiter.positivity_variables_cons - idp_positivity!(alpha, limiter, u, dt, semi, variable) + @trixi_timeit timer() "conservative variables" idp_positivity_conservative!(alpha, + limiter, + u, + dt, + semi, + variable) + end + + # Nonlinear variables + for variable in limiter.positivity_variables_nonlinear + @trixi_timeit timer() "nonlinear variables" idp_positivity_nonlinear!(alpha, + limiter, + u, dt, + semi, + variable) end return nothing end -@inline function idp_positivity!(alpha, limiter, u, dt, semi, variable) +############################################################################### +# Global positivity limiting of conservative variables + +@inline function idp_positivity_conservative!(alpha, limiter, u, dt, semi, variable) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes (; inverse_weights) = dg.basis @@ -256,7 +286,7 @@ end for j in eachnode(dg), i in eachnode(dg) var = u[variable, i, j, element] if var < 0 - error("Safe $variable is not safe. element=$element, node: $i $j, value=$var") + error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.") end # Compute bound @@ -302,4 +332,183 @@ end return nothing end + +@inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; positivity_correction_factor) = limiter + + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + var_min = variable_bounds[Symbol(string(variable), "_min")] + + @threaded for element in eachelement(dg, semi.cache) + inverse_jacobian = cache.elements.inverse_jacobian[element] + for j in eachnode(dg), i in eachnode(dg) + # Compute bound + u_local = get_node_vars(u, equations, dg, i, j, element) + var = variable(u_local, equations) + if var < 0 + error("Safe low-order method produces negative value for variable $variable. Try a smaller time step.") + end + var_min[i, j, element] = positivity_correction_factor * var + + # Perform Newton's bisection method to find new alpha + newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element, + variable, initial_check_nonnegative_newton_idp, + final_check_nonnegative_newton_idp, inverse_jacobian, + dt, equations, dg, cache, limiter) + end + end + + return nothing +end + +@inline function newton_loops_alpha!(alpha, bound, u, i, j, element, variable, + initial_check, final_check, inverse_jacobian, dt, + equations, dg, cache, limiter) + (; inverse_weights) = dg.basis + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes + + (; gamma_constant_newton) = limiter + + # negative xi direction + antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[i] * + get_node_vars(antidiffusive_flux1_R, equations, dg, i, j, + element) + newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, + equations, dt, limiter, antidiffusive_flux) + + # positive xi direction + antidiffusive_flux = -gamma_constant_newton * inverse_jacobian * + inverse_weights[i] * + get_node_vars(antidiffusive_flux1_L, equations, dg, i + 1, j, + element) + newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, + equations, dt, limiter, antidiffusive_flux) + + # negative eta direction + antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[j] * + get_node_vars(antidiffusive_flux2_R, equations, dg, i, j, + element) + newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, + equations, dt, limiter, antidiffusive_flux) + + # positive eta direction + antidiffusive_flux = -gamma_constant_newton * inverse_jacobian * + inverse_weights[j] * + get_node_vars(antidiffusive_flux2_L, equations, dg, i, j + 1, + element) + newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, + equations, dt, limiter, antidiffusive_flux) + + return nothing +end + +@inline function newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, + final_check, equations, dt, limiter, antidiffusive_flux) + newton_reltol, newton_abstol = limiter.newton_tolerances + + beta = 1 - alpha[i, j, element] + + beta_L = 0 # alpha = 1 + beta_R = beta # No higher beta (lower alpha) than the current one + + u_curr = u + beta * dt * antidiffusive_flux + + # If state is valid, perform initial check and return if correction is not needed + if isvalid(u_curr, equations) + goal = goal_function_newton_idp(variable, bound, u_curr, equations) + + initial_check(bound, goal, newton_abstol) && return nothing + end + + # Newton iterations + for iter in 1:(limiter.max_iterations_newton) + beta_old = beta + + # If the state is valid, evaluate d(goal)/d(beta) + if isvalid(u_curr, equations) + dgoal_dbeta = dgoal_function_newton_idp(variable, u_curr, dt, + antidiffusive_flux, equations) + else # Otherwise, perform a bisection step + dgoal_dbeta = 0 + end + + if dgoal_dbeta != 0 + # Update beta with Newton's method + beta = beta - goal / dgoal_dbeta + end + + # Check bounds + if (beta < beta_L) || (beta > beta_R) || (dgoal_dbeta == 0) || isnan(beta) + # Out of bounds, do a bisection step + beta = 0.5 * (beta_L + beta_R) + # Get new u + u_curr = u + beta * dt * antidiffusive_flux + + # If the state is invalid, finish bisection step without checking tolerance and iterate further + if !isvalid(u_curr, equations) + beta_R = beta + continue + end + + # Check new beta for condition and update bounds + goal = goal_function_newton_idp(variable, bound, u_curr, equations) + if initial_check(bound, goal, newton_abstol) + # New beta fulfills condition + beta_L = beta + else + # New beta does not fulfill condition + beta_R = beta + end + else + # Get new u + u_curr = u + beta * dt * antidiffusive_flux + + # If the state is invalid, redefine right bound without checking tolerance and iterate further + if !isvalid(u_curr, equations) + beta_R = beta + continue + end + + # Evaluate goal function + goal = goal_function_newton_idp(variable, bound, u_curr, equations) + end + + # Check relative tolerance + if abs(beta_old - beta) <= newton_reltol + break + end + + # Check absolute tolerance + if final_check(bound, goal, newton_abstol) + break + end + end + + new_alpha = 1 - beta + if alpha[i, j, element] > new_alpha + newton_abstol + error("Alpha is getting smaller. old: $(alpha[i, j, element]), new: $new_alpha") + else + alpha[i, j, element] = new_alpha + end + + return nothing +end + +### Auxiliary routines for Newton's bisection method ### +# Initial checks +@inline initial_check_nonnegative_newton_idp(bound, goal, newton_abstol) = goal <= 0 + +# Goal and d(Goal)d(u) function +@inline goal_function_newton_idp(variable, bound, u, equations) = bound - + variable(u, equations) +@inline function dgoal_function_newton_idp(variable, u, dt, antidiffusive_flux, + equations) + -dot(gradient_conservative(variable, u, equations), dt * antidiffusive_flux) +end + +# Final checks +@inline function final_check_nonnegative_newton_idp(bound, goal, newton_abstol) + (goal <= eps()) && (goal > -max(newton_abstol, abs(bound) * newton_abstol)) +end end # @muladd diff --git a/test/test_mpi.jl b/test/test_mpi.jl index ad1ba4e835d..001d9bff86e 100644 --- a/test/test_mpi.jl +++ b/test/test_mpi.jl @@ -8,6 +8,7 @@ include("test_trixi.jl") # Start with a clean environment: remove Trixi.jl output directory if it exists outdir = "out" Trixi.mpi_isroot() && isdir(outdir) && rm(outdir, recursive = true) +Trixi.MPI.Barrier(Trixi.mpi_comm()) # CI with MPI and some tests fails often on Windows. Thus, we check whether this # is the case here. We use GitHub Actions, so we can check whether we run CI @@ -19,10 +20,12 @@ CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows() # TreeMesh tests include("test_mpi_tree.jl") - # P4estMesh tests + # P4estMesh and T8codeMesh tests include("test_mpi_p4est_2d.jl") + include("test_mpi_t8code_2d.jl") if !CI_ON_WINDOWS # see comment on `CI_ON_WINDOWS` above include("test_mpi_p4est_3d.jl") + include("test_mpi_t8code_3d.jl") end end # MPI @@ -43,5 +46,6 @@ end # MPI supporting functionality # Clean up afterwards: delete Trixi.jl output directory Trixi.mpi_isroot() && @test_nowarn rm(outdir, recursive = true) +Trixi.MPI.Barrier(Trixi.mpi_comm()) end # module diff --git a/test/test_mpi_p4est_2d.jl b/test/test_mpi_p4est_2d.jl index da90537fcfd..6d66bc68a26 100644 --- a/test/test_mpi_p4est_2d.jl +++ b/test/test_mpi_p4est_2d.jl @@ -33,6 +33,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem") @test errors.linf≈[0.00011787417954578494] rtol=1.0e-4 end end + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_advection_nonconforming_flag.jl" begin @@ -40,6 +49,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem") "elixir_advection_nonconforming_flag.jl"), l2=[3.198940059144588e-5], linf=[0.00030636069494005547]) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_advection_unstructured_flag.jl" begin @@ -47,6 +65,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem") "elixir_advection_unstructured_flag.jl"), l2=[0.0005379687442422346], linf=[0.007438525029884735]) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_advection_amr_solution_independent.jl" begin @@ -56,6 +83,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem") l2=[4.949660644033807e-5], linf=[0.0004867846262313763], coverage_override=(maxiters = 6,)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin @@ -64,6 +100,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem") l2=[0.0012766060609964525], linf=[0.01750280631586159], coverage_override=(maxiters = 6,)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_advection_restart.jl" begin @@ -73,6 +118,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem") # With the default `maxiters = 1` in coverage tests, # there would be no time steps after the restart. coverage_override=(maxiters = 100_000,)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin @@ -90,6 +144,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem") 0.03759938693042297, 0.08039824959535657, ]) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end end end # P4estMesh MPI diff --git a/test/test_mpi_p4est_3d.jl b/test/test_mpi_p4est_3d.jl index 75f43650082..cca9093ec51 100644 --- a/test/test_mpi_p4est_3d.jl +++ b/test/test_mpi_p4est_3d.jl @@ -33,6 +33,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem") @test errors.linf≈[0.0014548839020096516] rtol=1.0e-4 end end + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_advection_amr.jl" begin @@ -46,6 +55,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem") initial_refinement_level = 2, base_level = 2, med_level = 3, max_level = 4)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin @@ -58,6 +76,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem") initial_refinement_level = 0, base_level = 0, med_level = 1, max_level = 2)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_advection_restart.jl" begin @@ -67,12 +94,30 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem") # With the default `maxiters = 1` in coverage tests, # there would be no time steps after the restart. coverage_override=(maxiters = 100_000,)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_advection_cubed_sphere.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_cubed_sphere.jl"), l2=[0.002006918015656413], linf=[0.027655117058380085]) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end # Compressible Euler @@ -94,6 +139,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem") 0.008526972236273522, ], tspan=(0.0, 0.01)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin @@ -114,6 +168,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem") 0.01562861968368434, ], tspan=(0.0, 1.0)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_euler_ec.jl" begin @@ -134,6 +197,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem") ], tspan=(0.0, 0.2), coverage_override=(polydeg = 3,)) # Prevent long compile time in CI + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_euler_source_terms_nonperiodic_hohqmesh.jl" begin @@ -153,6 +225,15 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem") 0.048396544302230504, 0.1154589758186293, ]) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end end end # P4estMesh MPI diff --git a/test/test_mpi_t8code_2d.jl b/test/test_mpi_t8code_2d.jl new file mode 100644 index 00000000000..7c7fc03898c --- /dev/null +++ b/test/test_mpi_t8code_2d.jl @@ -0,0 +1,142 @@ +module TestExamplesMPIT8codeMesh2D + +using Test +using Trixi + +include("test_trixi.jl") + +const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_dgsem") + +@testset "T8codeMesh MPI 2D" begin +#! format: noindent + +# Run basic tests +@testset "Examples 2D" begin + # Linear scalar advection + @trixi_testset "elixir_advection_basic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[8.311947673061856e-6], + linf=[6.627000273229378e-5]) + + @testset "error-based step size control" begin + Trixi.mpi_isroot() && println("-"^100) + Trixi.mpi_isroot() && + println("elixir_advection_basic.jl with error-based step size control") + + sol = solve(ode, RDPK3SpFSAL35(); abstol = 1.0e-4, reltol = 1.0e-4, + ode_default_options()..., callback = callbacks) + summary_callback() + errors = analysis_callback(sol) + if Trixi.mpi_isroot() + @test errors.l2≈[3.3022040342579066e-5] rtol=1.0e-4 + @test errors.linf≈[0.00011787417954578494] rtol=1.0e-4 + end + end + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_advection_nonconforming_flag.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_nonconforming_flag.jl"), + l2=[3.198940059144588e-5], + linf=[0.00030636069494005547]) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_advection_unstructured_flag.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_unstructured_flag.jl"), + l2=[0.0005379687442422346], + linf=[0.007438525029884735]) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_advection_amr_solution_independent.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_amr_solution_independent.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[4.933027431215839e-5], + linf=[0.00048678461161243136], + coverage_override=(maxiters = 6,)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_amr_unstructured_flag.jl"), + l2=[0.001980652042312077], + linf=[0.0328882442132265], + coverage_override=(maxiters = 6,)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"), + l2=[ + 0.0034516244508588046, + 0.0023420334036925493, + 0.0024261923964557187, + 0.004731710454271893, + ], + linf=[ + 0.04155789011775046, + 0.024772109862748914, + 0.03759938693042297, + 0.08039824959535657, + ]) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end +end # T8codeMesh MPI + +end # module diff --git a/test/test_mpi_t8code_3d.jl b/test/test_mpi_t8code_3d.jl new file mode 100644 index 00000000000..a15690a7629 --- /dev/null +++ b/test/test_mpi_t8code_3d.jl @@ -0,0 +1,180 @@ +module TestExamplesMPIT8codeMesh3D + +using Test +using Trixi + +include("test_trixi.jl") + +const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_dgsem") + +@testset "T8codeMesh MPI 3D" begin +#! format: noindent + +# Run basic tests +@testset "Examples 3D" begin + # Linear scalar advection + @trixi_testset "elixir_advection_basic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[0.00016263963870641478], + linf=[0.0014537194925779984]) + + @testset "error-based step size control" begin + Trixi.mpi_isroot() && println("-"^100) + Trixi.mpi_isroot() && + println("elixir_advection_basic.jl with error-based step size control") + + sol = solve(ode, RDPK3SpFSAL35(); abstol = 1.0e-4, reltol = 1.0e-4, + ode_default_options()..., callback = callbacks) + summary_callback() + errors = analysis_callback(sol) + if Trixi.mpi_isroot() + @test errors.l2≈[0.00016800412839949264] rtol=1.0e-4 + @test errors.linf≈[0.0014548839020096516] rtol=1.0e-4 + end + end + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_advection_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[1.1302812803902801e-5], + linf=[0.0007889950196294793], + # override values are different from the serial tests to ensure each process holds at least + # one element, otherwise OrdinaryDiffEq fails during initialization + coverage_override=(maxiters = 6, + initial_refinement_level = 2, + base_level = 2, med_level = 3, + max_level = 4)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_amr_unstructured_curved.jl"), + l2=[2.0556575425846923e-5], + linf=[0.00105682693484822], + tspan=(0.0, 1.0), + coverage_override=(maxiters = 6, + initial_refinement_level = 0, + base_level = 0, med_level = 1, + max_level = 2)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + # Compressible Euler + @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_curved.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_nonconforming_unstructured_curved.jl"), + l2=[ + 4.070355207909268e-5, + 4.4993257426833716e-5, + 5.10588457841744e-5, + 5.102840924036687e-5, + 0.00019986264001630542, + ], + linf=[ + 0.0016987332417202072, + 0.003622956808262634, + 0.002029576258317789, + 0.0024206977281964193, + 0.008526972236273522, + ], + tspan=(0.0, 0.01)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_nonperiodic.jl"), + l2=[ + 0.0015106060984283647, + 0.0014733349038567685, + 0.00147333490385685, + 0.001473334903856929, + 0.0028149479453087093, + ], + linf=[ + 0.008070806335238156, + 0.009007245083113125, + 0.009007245083121784, + 0.009007245083102688, + 0.01562861968368434, + ], + tspan=(0.0, 1.0)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_euler_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), + l2=[ + 0.010380390326164493, + 0.006192950051354618, + 0.005970674274073704, + 0.005965831290564327, + 0.02628875593094754, + ], + linf=[ + 0.3326911600075694, + 0.2824952141320467, + 0.41401037398065543, + 0.45574161423218573, + 0.8099577682187109, + ], + tspan=(0.0, 0.2), + coverage_override=(polydeg = 3,)) # Prevent long compile time in CI + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end +end # T8codeMesh MPI + +end # module diff --git a/test/test_mpi_tree.jl b/test/test_mpi_tree.jl index 0831f6a1313..6351a405b5d 100644 --- a/test/test_mpi_tree.jl +++ b/test/test_mpi_tree.jl @@ -76,7 +76,8 @@ CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows() # Here, we also test that SaveSolutionCallback prints multiple mesh files with AMR # Start with a clean environment: remove Trixi.jl output directory if it exists outdir = "out" - isdir(outdir) && rm(outdir, recursive = true) + Trixi.mpi_isroot() && isdir(outdir) && rm(outdir, recursive = true) + Trixi.MPI.Barrier(Trixi.mpi_comm()) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_refine_twice.jl"), l2=[0.00020547512522578292], diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 563bfda5f7a..f5fb169033a 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -31,14 +31,34 @@ end @trixi_testset "elixir_advection_coupled.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled.jl"), - l2=[7.816742843181738e-6, 7.816742843196112e-6], - linf=[6.314906965543265e-5, 6.314906965410039e-5], + l2=[ + 7.816742843336293e-6, + 7.816742843340186e-6, + 7.816742843025513e-6, + 7.816742843061526e-6, + ], + linf=[ + 6.314906965276812e-5, + 6.314906965187994e-5, + 6.31490696496595e-5, + 6.314906965032563e-5, + ], coverage_override=(maxiters = 10^5,)) @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin errors = analysis_callback(sol) - @test errors.l2≈[7.816742843181738e-6, 7.816742843196112e-6] rtol=1.0e-4 - @test errors.linf≈[6.314906965543265e-5, 6.314906965410039e-5] rtol=1.0e-4 + @test errors.l2≈[ + 7.816742843336293e-6, + 7.816742843340186e-6, + 7.816742843025513e-6, + 7.816742843061526e-6, + ] rtol=1.0e-4 + @test errors.linf≈[ + 6.314906965276812e-5, + 6.314906965187994e-5, + 6.31490696496595e-5, + 6.314906965032563e-5, + ] rtol=1.0e-4 # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl index ab95e068d02..d536a6dd73a 100644 --- a/test/test_t8code_2d.jl +++ b/test/test_t8code_2d.jl @@ -33,10 +33,9 @@ end @trixi_testset "test check_for_negative_volumes" begin @test_warn "Discovered negative volumes" begin # Unstructured mesh with six cells which have left-handed node ordering. - mesh_file = joinpath(EXAMPLES_DIR, "rectangle_with_negative_volumes.msh") - isfile(mesh_file) || - download("https://gist.githubusercontent.com/jmark/bfe0d45f8e369298d6cc637733819013/raw/cecf86edecc736e8b3e06e354c494b2052d41f7a/rectangle_with_negative_volumes.msh", - mesh_file) + mesh_file = Trixi.download("https://gist.githubusercontent.com/jmark/bfe0d45f8e369298d6cc637733819013/raw/cecf86edecc736e8b3e06e354c494b2052d41f7a/rectangle_with_negative_volumes.msh", + joinpath(EXAMPLES_DIR, + "rectangle_with_negative_volumes.msh")) # This call should throw a warning about negative volumes detected. mesh = T8codeMesh(mesh_file, 2) diff --git a/test/test_threaded.jl b/test/test_threaded.jl index dbbcbf4c7ce..a8a1b1b425a 100644 --- a/test/test_threaded.jl +++ b/test/test_threaded.jl @@ -8,6 +8,7 @@ include("test_trixi.jl") # Start with a clean environment: remove Trixi.jl output directory if it exists outdir = "out" Trixi.mpi_isroot() && isdir(outdir) && rm(outdir, recursive = true) +Trixi.MPI.Barrier(Trixi.mpi_comm()) @testset "Threaded tests" begin #! format: noindent @@ -471,5 +472,6 @@ end # Clean up afterwards: delete Trixi.jl output directory Trixi.mpi_isroot() && isdir(outdir) && @test_nowarn rm(outdir, recursive = true) +Trixi.MPI.Barrier(Trixi.mpi_comm()) end # module diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 61b5c54b5e9..b937abe92c0 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -581,6 +581,32 @@ end end end +@trixi_testset "elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl"), + l2=[ + 0.42185634563805724, + 0.1686471269704017, + 0.18240674916968103, + 0.17858250604280654, + ], + linf=[ + 1.7012978064377158, + 0.7149714986746726, + 0.5822547982757897, + 0.7300051017382696, + ], + tspan=(0.0, 2.0)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 + end +end + @trixi_testset "elixir_euler_colliding_flow.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_colliding_flow.jl"), l2=[ diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl index 953c077c0a3..1f8458075aa 100644 --- a/test/test_tree_2d_mhd.jl +++ b/test/test_tree_2d_mhd.jl @@ -332,24 +332,28 @@ end @trixi_testset "elixir_mhd_shockcapturing_subcell.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_subcell.jl"), - l2=[2.9974425783503109e-02, - 7.2849646345685956e-02, - 7.2488477174662239e-02, + l2=[ + 3.2064026219236076e-02, + 7.2461094392606618e-02, + 7.2380202888062711e-02, 0.0000000000000000e+00, - 1.2507971380965512e+00, - 1.8929505145499678e-02, - 1.2218606317164420e-02, + 8.6293936673145932e-01, + 8.4091669534557805e-03, + 5.2156364913231732e-03, 0.0000000000000000e+00, - 3.0154796910479838e-03], - linf=[3.2147382412340830e-01, - 1.3709471664007811e+00, - 1.3465154685288383e+00, + 2.0786952301129021e-04, + ], + linf=[ + 3.8778760255775635e-01, + 9.4666683953698927e-01, + 9.4618924645661928e-01, 0.0000000000000000e+00, - 1.6051257523415284e+01, - 3.0564266749926644e-01, - 2.3908016329805595e-01, + 1.0980297261521951e+01, + 1.0264404591009069e-01, + 1.0655686942176350e-01, 0.0000000000000000e+00, - 1.3711262178549158e-01], + 6.1013422157115546e-03, + ], tspan=(0.0, 0.003)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_unit.jl b/test/test_unit.jl index 6c8f5047d97..e5592979cca 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -416,7 +416,8 @@ end indicator_hg = IndicatorHennemannGassner(1.0, 0.0, true, "variable", "cache") @test_nowarn show(stdout, indicator_hg) - limiter_idp = SubcellLimiterIDP(true, [1], true, [1], 0.1, "cache") + limiter_idp = SubcellLimiterIDP(true, [1], true, [1], ["variable"], 0.1, "cache", 1, + (1.0, 1.0), 1.0) @test_nowarn show(stdout, limiter_idp) indicator_loehner = IndicatorLöhner(1.0, "variable", (; cache = nothing)) @@ -1215,6 +1216,26 @@ end end end +@testset "Consistency check for `gradient_conservative` routine" begin + # Set up conservative variables, equations + u = [ + 0.5011914484393387, + 0.8829127712445113, + 0.43024132987932817, + 0.7560616633050348, + ] + + equations = CompressibleEulerEquations2D(1.4) + + # Define wrapper function for pressure in order to call default implementation + function pressure_test(u, equations) + return pressure(u, equations) + end + + @test Trixi.gradient_conservative(pressure_test, u, equations) ≈ + Trixi.gradient_conservative(pressure, u, equations) +end + @testset "Equivalent Fluxes" begin # Set up equations and dummy conservative variables state # Burgers' Equation