Merge branch 'main' into subcell-limiting-positivity-nonlinear

Project.toml

@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-version = "0.6.1-pre"
+version = "0.6.2-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"

examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl

@@ -0,0 +1,63 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the polytropic Euler equations
+
+gamma = 1.4
+kappa = 0.5     # Scaling factor for the pressure.
+equations = PolytropicEulerEquations2D(gamma, kappa)
+
+initial_condition = initial_condition_convergence_test
+
+volume_flux = flux_winters_etal
+solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_davis),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 2,
+                periodicity = true,
+                n_cells_max = 30_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:l2_error_primitive,
+                                                              :linf_error_primitive))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.1)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

src/equations/ideal_glm_mhd_1d.jl

@@ -259,6 +259,148 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat
     return SVector(f1, f2, f3, f4, f5, f6, f7, f8)
 end
 
+"""
+    flux_hllc(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations1D)
+
+- Li (2005)
+An HLLC Riemann solver for magneto-hydrodynamics
+[DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020).
+"""
+function flux_hllc(u_ll, u_rr, orientation::Integer,
+                   equations::IdealGlmMhdEquations1D)
+    # Unpack left and right states
+    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)
+
+    # Conserved variables
+    rho_v1_ll = u_ll[2]
+    rho_v2_ll = u_ll[3]
+    rho_v3_ll = u_ll[4]
+
+    rho_v1_rr = u_rr[2]
+    rho_v2_rr = u_rr[3]
+    rho_v3_rr = u_rr[4]
+
+    # Obtain left and right fluxes
+    f_ll = flux(u_ll, orientation, equations)
+    f_rr = flux(u_rr, orientation, equations)
+
+    SsL, SsR = min_max_speed_naive(u_ll, u_rr, orientation, equations)
+    sMu_L = SsL - v1_ll
+    sMu_R = SsR - v1_rr
+    if SsL >= 0
+        f1 = f_ll[1]
+        f2 = f_ll[2]
+        f3 = f_ll[3]
+        f4 = f_ll[4]
+        f5 = f_ll[5]
+        f6 = f_ll[6]
+        f7 = f_ll[7]
+        f8 = f_ll[8]
+    elseif SsR <= 0
+        f1 = f_rr[1]
+        f2 = f_rr[2]
+        f3 = f_rr[3]
+        f4 = f_rr[4]
+        f5 = f_rr[5]
+        f6 = f_rr[6]
+        f7 = f_rr[7]
+        f8 = f_rr[8]
+    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 ) /
+                (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) /
+                (rho_rr * sMu_R - rho_ll * sMu_L)
+
+        Sdiff = SsR - SsL
+
+        # Compute HLL values for vStar, BStar
+        # These correspond to eq. (28) and (30) from the referenced paper 
+        # and the classic HLL intermediate state given by (2)
+        rho_HLL = (SsR * rho_rr - SsL * rho_ll - (f_rr[1] - f_ll[1])) / Sdiff
+
+        v1Star = (SsR * rho_v1_rr - SsL * rho_v1_ll - (f_rr[2] - f_ll[2])) /
+                 (Sdiff * rho_HLL)
+        v2Star = (SsR * rho_v2_rr - SsL * rho_v2_ll - (f_rr[3] - f_ll[3])) /
+                 (Sdiff * rho_HLL)
+        v3Star = (SsR * rho_v3_rr - SsL * rho_v3_ll - (f_rr[4] - f_ll[4])) /
+                 (Sdiff * rho_HLL)
+
+        #B1Star = (SsR * B1_rr - SsL * B1_ll - (f_rr[6] - f_ll[6])) / Sdiff
+        # 1D B1 = constant => B1_ll = B1_rr = B1Star
+        B1Star = B1_ll
+
+        B2Star = (SsR * B2_rr - SsL * B2_ll - (f_rr[7] - f_ll[7])) / Sdiff
+        B3Star = (SsR * B3_rr - SsL * B3_ll - (f_rr[8] - f_ll[8])) / Sdiff
+        if SsL <= SStar
+            SdiffStar = SsL - SStar
+
+            densStar = rho_ll * sMu_L / SdiffStar # (19)
+
+            mom_1_Star = densStar * SStar # (20)
+            mom_2_Star = densStar * v2_ll -
+                         (B1Star * B2Star - B1_ll * B2_ll) / SdiffStar # (21)
+            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)
+            # 1D B1 = constant => B1_ll = B1_rr = B1Star
+            pstar = rho_ll * sMu_L * (SStar - v1_ll) + p_ll # (17)
+
+            enerStar = u_ll[5] * sMu_L / SdiffStar +
+                       (pstar * SStar - p_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)
+
+            # Classic HLLC update (32)
+            f1 = f_ll[1] + SsL * (densStar - u_ll[1])
+            f2 = f_ll[2] + SsL * (mom_1_Star - u_ll[2])
+            f3 = f_ll[3] + SsL * (mom_2_Star - u_ll[3])
+            f4 = f_ll[4] + SsL * (mom_3_Star - u_ll[4])
+            f5 = f_ll[5] + SsL * (enerStar - u_ll[5])
+            f6 = f_ll[6] + SsL * (B1Star - u_ll[6])
+            f7 = f_ll[7] + SsL * (B2Star - u_ll[7])
+            f8 = f_ll[8] + SsL * (B3Star - u_ll[8])
+        else # SStar <= Ssr
+            SdiffStar = SsR - SStar
+
+            densStar = rho_rr * sMu_R / SdiffStar # (19)
+
+            mom_1_Star = densStar * SStar # (20)
+            mom_2_Star = densStar * v2_rr -
+                         (B1Star * B2Star - B1_rr * B2_rr) / SdiffStar # (21)
+            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)
+            # 1D B1 = constant => B1_ll = B1_rr = B1Star
+            pstar = rho_rr * sMu_R * (SStar - v1_rr) + p_rr # (17)
+
+            enerStar = u_rr[5] * sMu_R / SdiffStar +
+                       (pstar * SStar - p_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)
+
+            # Classic HLLC update (32)           
+            f1 = f_rr[1] + SsR * (densStar - u_rr[1])
+            f2 = f_rr[2] + SsR * (mom_1_Star - u_rr[2])
+            f3 = f_rr[3] + SsR * (mom_2_Star - u_rr[3])
+            f4 = f_rr[4] + SsR * (mom_3_Star - u_rr[4])
+            f5 = f_rr[5] + SsR * (enerStar - u_rr[5])
+            f6 = f_rr[6] + SsR * (B1Star - u_rr[6])
+            f7 = f_rr[7] + SsR * (B2Star - u_rr[7])
+            f8 = f_rr[8] + SsR * (B3Star - u_rr[8])
+        end
+    end
+    return SVector(f1, f2, f3, f4, f5, f6, f7, f8)
+end
+
 # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
 @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
                                      equations::IdealGlmMhdEquations1D)

src/equations/polytropic_euler_2d.jl

@@ -120,7 +120,7 @@ function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquat
     return prim2cons(SVector(rho, v1, v2), equations)
 end
 
-# Calculate 1D flux for a single point in the normal direction
+# Calculate 2D flux for a single point in the normal direction
 # Note, this directional vector is not normalized
 @inline function flux(u, normal_direction::AbstractVector,
                       equations::PolytropicEulerEquations2D)
@@ -135,8 +135,28 @@ end
     return SVector(f1, f2, f3)
 end
 
+# Calculate 2D flux for a single point
+@inline function flux(u, orientation::Integer, equations::PolytropicEulerEquations2D)
+    _, v1, v2 = cons2prim(u, equations)
+    p = pressure(u, equations)
+
+    rho_v1 = u[2]
+    rho_v2 = u[3]
+
+    if orientation == 1
+        f1 = rho_v1
+        f2 = rho_v1 * v1 + p
+        f3 = rho_v1 * v2
+    else
+        f1 = rho_v2
+        f2 = rho_v2 * v1
+        f3 = rho_v2 * v2 + p
+    end
+    return SVector(f1, f2, f3)
+end
+
 """
-    flux_winters_etal(u_ll, u_rr, normal_direction,
+    flux_winters_etal(u_ll, u_rr, orientation_or_normal_direction,
                       equations::PolytropicEulerEquations2D)
 
 Entropy conserving two-point flux for isothermal or polytropic gases.
@@ -178,6 +198,37 @@ For details see Section 3.2 of the following reference
     return SVector(f1, f2, f3)
 end
 
+@inline function flux_winters_etal(u_ll, u_rr, orientation::Integer,
+                                   equations::PolytropicEulerEquations2D)
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+    p_ll = equations.kappa * rho_ll^equations.gamma
+    p_rr = equations.kappa * rho_rr^equations.gamma
+
+    # Compute the necessary mean values
+    if equations.gamma == 1.0 # isothermal gas
+        rho_mean = ln_mean(rho_ll, rho_rr)
+    else # equations.gamma > 1 # polytropic gas
+        rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma)
+    end
+    v1_avg = 0.5 * (v1_ll + v1_rr)
+    v2_avg = 0.5 * (v2_ll + v2_rr)
+    p_avg = 0.5 * (p_ll + p_rr)
+
+    if orientation == 1 # x-direction
+        f1 = rho_mean * 0.5 * (v1_ll + v1_rr)
+        f2 = f1 * v1_avg + p_avg
+        f3 = f1 * v2_avg
+    else # y-direction
+        f1 = rho_mean * 0.5 * (v2_ll + v2_rr)
+        f2 = f1 * v1_avg
+        f3 = f1 * v2_avg + p_avg
+    end
+
+    return SVector(f1, f2, f3)
+end
+
 @inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
                                      equations::PolytropicEulerEquations2D)
     rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
@@ -196,6 +247,53 @@ end
     return lambda_min, lambda_max
 end
 
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
+                                     equations::PolytropicEulerEquations2D)
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+    # Pressure for polytropic Euler
+    p_ll = equations.kappa * rho_ll^equations.gamma
+    p_rr = equations.kappa * rho_rr^equations.gamma
+
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    if orientation == 1 # x-direction
+        λ_min = min(v1_ll - c_ll, v1_rr - c_rr)
+        λ_max = max(v1_ll + c_ll, v1_rr + c_rr)
+    else # y-direction
+        λ_min = min(v2_ll - c_ll, v2_rr - c_rr)
+        λ_max = max(v2_ll + c_ll, v2_rr + c_rr)
+    end
+
+    return λ_min, λ_max
+end
+
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::PolytropicEulerEquations2D)
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+    # Pressure for polytropic Euler
+    p_ll = equations.kappa * rho_ll^equations.gamma
+    p_rr = equations.kappa * rho_rr^equations.gamma
+
+    norm_ = norm(normal_direction)
+
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_
+
+    v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
+
+    # The v_normals are already scaled by the norm
+    λ_min = min(v_normal_ll - c_ll, v_normal_rr - c_rr)
+    λ_max = max(v_normal_ll + c_ll, v_normal_rr + c_rr)
+
+    return λ_min, λ_max
+end
+
 @inline function max_abs_speeds(u, equations::PolytropicEulerEquations2D)
     rho, v1, v2 = cons2prim(u, equations)
     c = sqrt(equations.gamma * equations.kappa * rho^(equations.gamma - 1))

test/Project.toml

@@ -12,7 +12,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-Aqua = "0.7"
+Aqua = "0.8"
 CairoMakie = "0.6, 0.7, 0.8, 0.9, 0.10"
 Downloads = "1"
 ForwardDiff = "0.10"

test/test_aqua.jl

@@ -11,8 +11,8 @@ include("test_trixi.jl")
                   ambiguities = false,
                   # exceptions necessary for adding a new method `StartUpDG.estimate_h`
                   # in src/solvers/dgmulti/sbp.jl
-                  piracy = (treat_as_own = [Trixi.StartUpDG.RefElemData,
-                                Trixi.StartUpDG.MeshData],))
+                  piracies = (treat_as_own = [Trixi.StartUpDG.RefElemData,
+                                  Trixi.StartUpDG.MeshData],))
 end
 
 end #module