Merge branch 'main' into subcell-limiting-p4est-all

bennibolm · Dec 3, 2023 · 0d39ec4 · 0d39ec4
2 parents 000e2da + 17d5070
commit 0d39ec4
Show file tree

Hide file tree

Showing 3 changed files with 49 additions and 19 deletions.
diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl
@@ -272,6 +272,10 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
     rho_ll, v1_ll, v2_ll, v3_ll, p_ll, B1_ll, B2_ll, B3_ll = cons2prim(u_ll, equations)
     rho_rr, v1_rr, v2_rr, v3_rr, p_rr, B1_rr, B2_rr, B3_rr = cons2prim(u_rr, equations)
 
+    # Total pressure, i.e., thermal + magnetic pressures (eq. (12))
+    p_tot_ll = p_ll + 0.5 * (B1_ll^2 + B2_ll^2 + B3_ll^2)
+    p_tot_rr = p_rr + 0.5 * (B1_rr^2 + B2_rr^2 + B3_rr^2)
+
     # Conserved variables
     rho_v1_ll = u_ll[2]
     rho_v2_ll = u_ll[3]
@@ -309,11 +313,11 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
     else
         # Compute the "HLLC-speed", eq. (14) from paper mentioned above
         #=
-        SStar = (rho_rr * v1_rr * sMu_R - rho_ll * v1_ll * sMu_L + p_ll - p_rr - B1_ll^2 + B1_rr^2 ) /
+        SStar = (rho_rr * v1_rr * sMu_R - rho_ll * v1_ll * sMu_L + p_tot_ll - p_tot_rr - B1_ll^2 + B1_rr^2 ) /
                 (rho_rr * sMu_R - rho_ll * sMu_L)
         =#
         # Simplification for 1D: B1 is constant
-        SStar = (rho_rr * v1_rr * sMu_R - rho_ll * v1_ll * sMu_L + p_ll - p_rr) /
+        SStar = (rho_rr * v1_rr * sMu_R - rho_ll * v1_ll * sMu_L + p_tot_ll - p_tot_rr) /
                 (rho_rr * sMu_R - rho_ll * sMu_L)
 
         Sdiff = SsR - SsL
@@ -347,12 +351,12 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
             mom_3_Star = densStar * v3_ll -
                          (B1Star * B3Star - B1_ll * B3_ll) / SdiffStar # (22)
 
-            #pstar = rho_ll * sMu_L * (SStar - v1_ll) + p_ll - B1_ll^2 + B1Star^2 # (17)
+            #p_tot_Star = rho_ll * sMu_L * (SStar - v1_ll) + p_tot_ll - B1_ll^2 + B1Star^2 # (17)
             # 1D B1 = constant => B1_ll = B1_rr = B1Star
-            pstar = rho_ll * sMu_L * (SStar - v1_ll) + p_ll # (17)
+            p_tot_Star = rho_ll * sMu_L * (SStar - v1_ll) + p_tot_ll # (17)
 
             enerStar = u_ll[5] * sMu_L / SdiffStar +
-                       (pstar * SStar - p_ll * v1_ll - (B1Star *
+                       (p_tot_Star * SStar - p_tot_ll * v1_ll - (B1Star *
                          (B1Star * v1Star + B2Star * v2Star + B3Star * v3Star) -
                          B1_ll * (B1_ll * v1_ll + B2_ll * v2_ll + B3_ll * v3_ll))) /
                        SdiffStar # (23)
@@ -377,12 +381,12 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
             mom_3_Star = densStar * v3_rr -
                          (B1Star * B3Star - B1_rr * B3_rr) / SdiffStar # (22)
 
-            #pstar = rho_rr * sMu_R * (SStar - v1_rr) + p_rr - B1_rr^2 + B1Star^2 # (17)
+            #p_tot_Star = rho_rr * sMu_R * (SStar - v1_rr) + p_tot_rr - B1_rr^2 + B1Star^2 # (17)
             # 1D B1 = constant => B1_ll = B1_rr = B1Star
-            pstar = rho_rr * sMu_R * (SStar - v1_rr) + p_rr # (17)
+            p_tot_Star = rho_rr * sMu_R * (SStar - v1_rr) + p_tot_rr # (17)
 
             enerStar = u_rr[5] * sMu_R / SdiffStar +
-                       (pstar * SStar - p_rr * v1_rr - (B1Star *
+                       (p_tot_Star * SStar - p_tot_rr * v1_rr - (B1Star *
                          (B1Star * v1Star + B2Star * v2Star + B3Star * v3Star) -
                          B1_rr * (B1_rr * v1_rr + B2_rr * v2_rr + B3_rr * v3_rr))) /
                        SdiffStar # (23)

diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl
@@ -30,6 +30,7 @@ function save_mesh_file(mesh::TreeMesh, output_directory, timestep,
         attributes(file)["mesh_type"] = get_name(mesh)
         attributes(file)["ndims"] = ndims(mesh)
         attributes(file)["n_cells"] = n_cells
+        attributes(file)["capacity"] = mesh.tree.capacity
         attributes(file)["n_leaf_cells"] = count_leaf_cells(mesh.tree)
         attributes(file)["minimum_level"] = minimum_level(mesh.tree)
         attributes(file)["maximum_level"] = maximum_level(mesh.tree)
@@ -249,10 +250,10 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT)
     end
 
     if mesh_type == "TreeMesh"
-        n_cells = h5open(mesh_file, "r") do file
-            return read(attributes(file)["n_cells"])
+        capacity = h5open(mesh_file, "r") do file
+            return read(attributes(file)["capacity"])
         end
-        mesh = TreeMesh(SerialTree{ndims}, max(n_cells, n_cells_max))
+        mesh = TreeMesh(SerialTree{ndims}, max(n_cells_max, capacity))
         load_mesh!(mesh, mesh_file)
     elseif mesh_type == "StructuredMesh"
         size_, mapping_as_string = h5open(mesh_file, "r") do file

diff --git a/test/test_tree_1d_mhd.jl b/test/test_tree_1d_mhd.jl
@@ -109,6 +109,31 @@ end
     end
 end
 
+@trixi_testset "elixir_mhd_alfven_wave.jl with flux_hllc" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
+                        l2=[
+                            1.036850596986597e-5, 1.965192583650368e-6,
+                            3.5882124656715505e-5, 3.5882124656638764e-5,
+                            5.270975504780837e-6, 1.1963224165731992e-16,
+                            3.595811808912869e-5, 3.5958118089159453e-5,
+                        ],
+                        linf=[
+                            2.887280521446378e-5, 7.310580790352001e-6,
+                            0.00012390046377899755, 0.00012390046377787345,
+                            1.5102711136583125e-5, 2.220446049250313e-16,
+                            0.0001261935452181312, 0.0001261935452182006,
+                        ],
+                        surface_flux=flux_hllc)
+    # 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_mhd_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
                         l2=[
@@ -210,16 +235,16 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_torrilhon_shock_tube.jl"),
                         surface_flux=flux_hllc,
                         l2=[
-                            0.45738965718253993, 0.479402222862685, 0.34069729746967664,
-                            0.44795514335568865, 0.9206813325913135,
-                            1.3216517820475193e-16, 0.2889672868491632,
-                            0.2552794220777942,
+                            0.4573799618744708, 0.4792633358230866, 0.34064852506872795,
+                            0.4479668434955162, 0.9203891782415092,
+                            1.3216517820475193e-16, 0.28887826520860815,
+                            0.255281629265771,
                         ],
                         linf=[
-                            1.2181099854251536, 0.8869319941747589, 0.8763562906332134,
-                            0.9712221036087284, 1.6734231113527818,
-                            2.220446049250313e-16, 0.7035011427822779,
-                            0.6562884129650286,
+                            1.2382842201671505, 0.8929169308132259, 0.871298623806198,
+                            0.9822415614542821, 1.6726170732132717,
+                            2.220446049250313e-16, 0.7016155888023747,
+                            0.6556091522071984,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)