Merge branch 'main' into ThreadPara_src_EG

NEWS.md

@@ -6,13 +6,27 @@ for human readability.
 
 ## Changes when updating to v0.9 from v0.8.x
 
+#### Added
+
+- Boundary conditions are now supported on nonconservative terms ([#2062]).
+
 #### Changed
 
 - We removed the first argument `semi` corresponding to a `Semidiscretization` from the 
   `AnalysisSurfaceIntegral` constructor, as it is no longer needed (see [#1959]).
-  The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`. ([#2069])
+  The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`.
+  ([#2069])
 - In functions `rhs!`, `rhs_parabolic!`  we removed the unused argument `initial_condition`. ([#2037])
   Users should not be affected by this.
+- Nonconservative terms depend only on `normal_direction_average` instead of both 
+  `normal_direction_average` and `normal_direction_ll`, such that the function signature is now 
+  identical with conservative fluxes. This required a change of the `normal_direction` in
+  `flux_nonconservative_powell` ([#2062]).
+
+#### Deprecated
+
+#### Removed
+
 
 ## Changes in the v0.8 lifecycle
 

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.8.11-DEV"
+version = "0.9.2-DEV"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

benchmark/Project.toml

@@ -8,4 +8,4 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 BenchmarkTools = "0.5, 0.7, 1.0"
 OrdinaryDiffEq = "5.65, 6"
 PkgBenchmark = "0.2.10"
-Trixi = "0.4, 0.5, 0.6, 0.7, 0.8"
+Trixi = "0.4, 0.5, 0.6, 0.7, 0.8, 0.9"

examples/structured_2d_dgsem/elixir_mhd_coupled.jl

@@ -31,14 +31,6 @@ equations = IdealGlmMhdEquations2D(gamma)
 
 cells_per_dimension = (32, 64)
 
-# Extend the definition of the non-conservative Powell flux functions.
-import Trixi.flux_nonconservative_powell
-function flux_nonconservative_powell(u_ll, u_rr,
-                                     normal_direction_ll::AbstractVector,
-                                     equations::IdealGlmMhdEquations2D)
-    flux_nonconservative_powell(u_ll, u_rr, normal_direction_ll, normal_direction_ll,
-                                equations)
-end
 volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
 solver = DGSEM(polydeg = 3,
                surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),

src/basic_types.jl

@@ -71,14 +71,14 @@ struct BoundaryConditionDoNothing end
                                                 orientation_or_normal_direction,
                                                 direction::Integer, x, t, surface_flux,
                                                 equations)
-    return flux(u_inner, orientation_or_normal_direction, equations)
+    return surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations)
 end
 
 # This version can be called by hyperbolic solvers on unstructured, curved meshes
 @inline function (::BoundaryConditionDoNothing)(u_inner,
                                                 outward_direction::AbstractVector,
                                                 x, t, surface_flux, equations)
-    return flux(u_inner, outward_direction, equations)
+    return surface_flux(u_inner, u_inner, outward_direction, equations)
 end
 
 # This version can be called by parabolic solvers

src/callbacks_step/analysis.jl

@@ -314,8 +314,8 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi)
         mpi_println(" #elements:      " *
                     @sprintf("% 14d", nelementsglobal(mesh, solver, cache)))
 
-        # Level information (only show for AMR)
-        print_amr_information(integrator.opts.callback, mesh, solver, cache)
+        # Level information (only for AMR and/or non-uniform `TreeMesh`es)
+        print_level_information(integrator.opts.callback, mesh, solver, cache)
         mpi_println()
 
         # Open file for appending and store time step and time information
@@ -489,21 +489,7 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
     return nothing
 end
 
-# Print level information only if AMR is enabled
-function print_amr_information(callbacks, mesh, solver, cache)
-
-    # Return early if there is nothing to print
-    uses_amr(callbacks) || return nothing
-
-    # Get global minimum and maximum level from the AMRController
-    min_level = max_level = 0
-    for cb in callbacks.discrete_callbacks
-        if cb.affect! isa AMRCallback
-            min_level = cb.affect!.controller.base_level
-            max_level = cb.affect!.controller.max_level
-        end
-    end
-
+function print_level_information(mesh, solver, cache, min_level, max_level)
     # Get local element count per level
     elements_per_level = get_elements_per_level(min_level, max_level, mesh, solver,
                                                 cache)
@@ -522,6 +508,45 @@ function print_amr_information(callbacks, mesh, solver, cache)
     return nothing
 end
 
+function print_level_information(callbacks, mesh::TreeMesh, solver, cache)
+    if uses_amr(callbacks)
+        # Get global minimum and maximum level from the AMRController
+        min_level = max_level = 0
+        for cb in callbacks.discrete_callbacks
+            if cb.affect! isa AMRCallback
+                min_level = cb.affect!.controller.base_level
+                max_level = cb.affect!.controller.max_level
+            end
+        end
+        print_level_information(mesh, solver, cache, min_level, max_level)
+        # Special check for `TreeMesh`es without AMR.
+        # These meshes may still be non-uniform due to `refinement_patches`, see 
+        # `refine_box!`, `coarsen_box!`, and `refine_sphere!`.
+    elseif minimum_level(mesh.tree) != maximum_level(mesh.tree)
+        min_level = minimum_level(mesh.tree)
+        max_level = maximum_level(mesh.tree)
+        print_level_information(mesh, solver, cache, min_level, max_level)
+    else # Uniform mesh
+        return nothing
+    end
+end
+
+function print_level_information(callbacks, mesh, solver, cache)
+    if uses_amr(callbacks)
+        # Get global minimum and maximum level from the AMRController
+        min_level = max_level = 0
+        for cb in callbacks.discrete_callbacks
+            if cb.affect! isa AMRCallback
+                min_level = cb.affect!.controller.base_level
+                max_level = cb.affect!.controller.max_level
+            end
+        end
+        print_level_information(mesh, solver, cache, min_level, max_level)
+    else # Uniform mesh
+        return nothing
+    end
+end
+
 function get_elements_per_level(min_level, max_level, mesh::P4estMesh, solver, cache)
     elements_per_level = zeros(P4EST_MAXLEVEL + 1)
 

src/equations/compressible_euler_1d.jl

@@ -412,10 +412,10 @@ See also
     return SVector(f1, f2, f3)
 end
 
-# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that 
-# the normal component is incorporated into the numerical flux. 
-# 
-# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a 
+# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that
+# the normal component is incorporated into the numerical flux.
+#
+# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a
 # similar implementation.
 @inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector,
                               equations::CompressibleEulerEquations1D)
@@ -943,6 +943,12 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::CompressibleEulerEquations1D)
+    rho = u[1]
+    v1 = u[2] / rho
+    return v1
+end
+
 @inline function pressure(u, equations::CompressibleEulerEquations1D)
     rho, rho_v1, rho_e = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2) / rho)

src/equations/compressible_euler_2d.jl

@@ -1957,6 +1957,19 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::CompressibleEulerEquations2D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    return SVector(v1, v2)
+end
+
+@inline function velocity(u, orientation::Int, equations::CompressibleEulerEquations2D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
 @inline function pressure(u, equations::CompressibleEulerEquations2D)
     rho, rho_v1, rho_v2, rho_e = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho)

src/equations/compressible_euler_3d.jl

@@ -1692,6 +1692,20 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::CompressibleEulerEquations3D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    return SVector(v1, v2, v3)
+end
+
+@inline function velocity(u, orientation::Int, equations::CompressibleEulerEquations3D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
 @inline function pressure(u, equations::CompressibleEulerEquations3D)
     rho, rho_v1, rho_v2, rho_v3, rho_e = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho)

src/equations/compressible_euler_multicomponent_1d.jl

@@ -204,13 +204,12 @@ function initial_condition_weak_blast_wave(x, t,
 
     prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5f0 ?
                                                                 2^(i - 1) * (1 - 2) /
+                                                                (RealT(1) -
+                                                                 2^ncomponents(equations)) :
+                                                                2^(i - 1) * (1 - 2) *
+                                                                RealT(1.1691) /
                                                                 (1 -
-                                                                 2^ncomponents(equations)) *
-                                                                one(RealT) :
-                                                                2^(i - 1) * (1 - 2) /
-                                                                (1 -
-                                                                 2^ncomponents(equations)) *
-                                                                convert(RealT, 1.1691)
+                                                                 2^ncomponents(equations))
                                                                 for i in eachcomponent(equations))
 
     v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi
@@ -618,4 +617,10 @@ end
     return SVector{ncomponents(equations), real(equations)}(u[i + 2] * v
                                                             for i in eachcomponent(equations))
 end
+
+@inline function velocity(u, equations::CompressibleEulerMulticomponentEquations1D)
+    rho = density(u, equations)
+    v1 = u[1] / rho
+    return v1
+end
 end # @muladd

src/equations/compressible_euler_multicomponent_2d.jl

@@ -222,13 +222,12 @@ function initial_condition_weak_blast_wave(x, t,
 
     prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5f0 ?
                                                                 2^(i - 1) * (1 - 2) /
+                                                                (RealT(1) -
+                                                                 2^ncomponents(equations)) :
+                                                                2^(i - 1) * (1 - 2) *
+                                                                RealT(1.1691) /
                                                                 (1 -
-                                                                 2^ncomponents(equations)) *
-                                                                one(RealT) :
-                                                                2^(i - 1) * (1 - 2) /
-                                                                (1 -
-                                                                 2^ncomponents(equations)) *
-                                                                convert(RealT, 1.1691)
+                                                                 2^ncomponents(equations))
                                                                 for i in eachcomponent(equations))
 
     v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi
@@ -825,4 +824,18 @@ end
     return SVector{ncomponents(equations), real(equations)}(u[i + 3] * v
                                                             for i in eachcomponent(equations))
 end
+
+@inline function velocity(u, equations::CompressibleEulerMulticomponentEquations2D)
+    rho = density(u, equations)
+    v1 = u[1] / rho
+    v2 = u[2] / rho
+    return SVector(v1, v2)
+end
+
+@inline function velocity(u, orientation::Int,
+                          equations::CompressibleEulerMulticomponentEquations2D)
+    rho = density(u, equations)
+    v = u[orientation] / rho
+    return v
+end
 end # @muladd

src/equations/equations.jl

@@ -307,6 +307,37 @@ The inverse conversion is performed by [`cons2prim`](@ref).
 """
 function prim2cons end
 
+"""
+    velocity(u, equations)
+
+Return the velocity vector corresponding to the equations, e.g., fluid velocity for
+Euler's equations. The velocity in certain orientation or normal direction (scalar) can be computed
+with `velocity(u, orientation, equations)` or `velocity(u, normal_direction, equations)`
+respectively. The `velocity(u, normal_direction, equations)` function calls
+`velocity(u, equations)` to compute the velocity vector and then normal vector, thus allowing
+a general function to be written for the AbstractEquations type. However, the
+`velocity(u, orientation, equations)` is written for each equation separately to ensure
+only the velocity in the desired direction (orientation) is computed.
+`u` is a vector of the conserved variables at a single node, i.e., a vector
+of the correct length `nvariables(equations)`.
+"""
+function velocity end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::AbstractEquations{2})
+    vel = velocity(u, equations)
+    v = vel[1] * normal_direction[1] + vel[2] * normal_direction[2]
+    return v
+end
+
+@inline function velocity(u, normal_direction::AbstractVector,
+                          equations::AbstractEquations{3})
+    vel = velocity(u, equations)
+    v = vel[1] * normal_direction[1] + vel[2] * normal_direction[2] +
+        vel[3] * normal_direction[3]
+    return v
+end
+
 """
     entropy(u, equations)
 

src/equations/ideal_glm_mhd_1d.jl

@@ -326,7 +326,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
         Sdiff = SsR - SsL
 
         # Compute HLL values for vStar, BStar
-        # These correspond to eq. (28) and (30) from the referenced paper 
+        # 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
 
@@ -394,7 +394,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer,
                          B1_rr * (B1_rr * v1_rr + B2_rr * v2_rr + B3_rr * v3_rr))) /
                        SdiffStar # (23)
 
-            # Classic HLLC update (32)           
+            # 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])
@@ -555,6 +555,20 @@ end
     return rho
 end
 
+@inline function velocity(u, equations::IdealGlmMhdEquations1D)
+    rho = u[1]
+    v1 = u[2] / rho
+    v2 = u[3] / rho
+    v3 = u[4] / rho
+    return SVector(v1, v2, v3)
+end
+
+@inline function velocity(u, orientation::Int, equations::IdealGlmMhdEquations1D)
+    rho = u[1]
+    v = u[orientation + 1] / rho
+    return v
+end
+
 @inline function pressure(u, equations::IdealGlmMhdEquations1D)
     rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u
     p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho