specialize calc_boundary_flux! for nonconservative terms for DGMulti

examples/dgmulti_2d/elixir_mhd_reflective_wall.jl

@@ -0,0 +1,105 @@
+
+using OrdinaryDiffEq
+using Trixi
+using LinearAlgebra: norm, dot # for use in the MHD boundary condition
+
+###############################################################################
+# semidiscretization of the compressible ideal GLM-MHD equations
+equations = IdealGlmMhdEquations2D(1.4)
+
+function initial_condition_perturbation(x, t, equations::IdealGlmMhdEquations2D)
+  # pressure perturbation in a vertically magnetized field on the domain [-1, 1]^2
+
+  r2 = (x[1] + 0.25)^2 + (x[2] + 0.25)^2
+
+  rho = 1.0
+  v1  = 0.0
+  v2  = 0.0
+  v3  = 0.0
+  p   = 1 + 0.5 * exp(-100 * r2)
+
+  # the pressure and magnetic field are chosen to be strongly
+  # magnetized, such that p / ||B||^2 ≈ 0.01.
+  B1  = 0.0
+  B2  = 40.0 / sqrt(4.0 * pi)
+  B3  = 0.0
+
+  psi = 0.0
+  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
+end
+initial_condition = initial_condition_perturbation
+
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell)
+volume_flux  = (flux_hindenlang_gassner, flux_nonconservative_powell)
+
+solver = DGMulti(polydeg=3, element_type = Quad(), approximation_type = GaussSBP(),
+                 surface_integral = SurfaceIntegralWeakForm(surface_flux),
+                 volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+x_neg(x, tol=50*eps()) = abs(x[1] + 1) < tol
+x_pos(x, tol=50*eps()) = abs(x[1] - 1) < tol
+y_neg(x, tol=50*eps()) = abs(x[2] + 1) < tol
+y_pos(x, tol=50*eps()) = abs(x[2] - 1) < tol
+is_on_boundary = Dict(:x_neg => x_neg, :x_pos => x_pos, :y_neg => y_neg, :y_pos => y_pos)
+
+cells_per_dimension = (16, 16)
+mesh = DGMultiMesh(solver, cells_per_dimension; periodicity=(false, false), is_on_boundary)
+
+# Create a "reflective-like" boundary condition by mirroring the velocity but leaving the magnetic field alone.
+# Note that this boundary condition is probably not entropy stable.
+function boundary_condition_velocity_slip_wall(u_inner, normal_direction::AbstractVector, x, t,
+                                               surface_flux_function, equations::IdealGlmMhdEquations2D)
+
+  # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later
+  norm_ = norm(normal_direction)
+  normal = normal_direction / norm_
+
+  # compute the primitive variables
+  rho, v1, v2, v3, p, B1, B2, B3, psi = cons2prim(u_inner, equations)
+
+  v_normal = dot(normal, SVector(v1, v2))
+  u_mirror = prim2cons(SVector(rho, v1 - 2 * v_normal * normal[1],
+                                    v2 - 2 * v_normal * normal[2],
+                                    v3, p, B1, B2, B3, psi), equations)
+
+  return surface_flux_function(u_inner, u_mirror, normal, equations) * norm_
+end
+
+boundary_conditions = (; x_neg=boundary_condition_velocity_slip_wall,
+                         x_pos=boundary_condition_velocity_slip_wall,
+                         y_neg=boundary_condition_do_nothing,
+                         y_pos=BoundaryConditionDirichlet(initial_condition))
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    boundary_conditions=boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.075)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(solver))
+alive_callback = AliveCallback(alive_interval=10)
+
+cfl = 0.5
+stepsize_callback = StepsizeCallback(cfl=cfl)
+glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        stepsize_callback,
+                        glm_speed_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1e-5, # 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/basic_types.jl

@@ -76,20 +76,22 @@ struct BoundaryConditionDoNothing end
 # This version can be called by hyperbolic solvers on logically Cartesian meshes
 @inline function (::BoundaryConditionDoNothing)(
     u_inner, orientation_or_normal_direction, direction::Integer, x, t, surface_flux, equations)
+
   return flux(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)
+@inline function (::BoundaryConditionDoNothing)(u_inner, outward_direction::AbstractVector,
+                                                x, t, surface_flux, equations)
+
   return flux(u_inner, outward_direction, equations)
 end
 
 # This version can be called by parabolic solvers
 @inline function (::BoundaryConditionDoNothing)(inner_flux_or_state, other_args...)
   return inner_flux_or_state
 end
-    
+
 """
     boundary_condition_do_nothing = Trixi.BoundaryConditionDoNothing()
 

src/callbacks_step/glm_speed_dg.jl

@@ -23,7 +23,7 @@ function calc_dt_for_cleaning_speed(cfl::Real, mesh,
 
   # Compute time step for GLM linear advection equation with c_h=1 for a DGMulti discretization.
   # Copies implementation behavior of `calc_dt_for_cleaning_speed` for DGSEM discretizations.
-  max_scaled_speed_for_c_h = (1 / minimum(md.J)) * ndims(equations)
+  max_scaled_speed_for_c_h = inv(minimum(md.J)) * ndims(equations)
 
   # This mimics `max_dt` for `TreeMesh`, except that `nnodes(dg)` is replaced by
   # `polydeg+1`. This is because `nnodes(dg)` returns the total number of

src/equations/equations.jl

@@ -171,7 +171,8 @@ end
 @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
                                                                   normal_direction::AbstractVector,
                                                                   x, t,
-                                                                  surface_flux_function, equations)
+                                                                  surface_flux_function,
+                                                                  equations)
   # get the external value of the solution
   u_boundary = boundary_condition.boundary_value_function(x, t, equations)
 
@@ -328,7 +329,7 @@ include("compressible_euler_multicomponent_2d.jl")
     eachcomponent(equations::AbstractCompressibleEulerMulticomponentEquations)
 
 Return an iterator over the indices that specify the location in relevant data structures
-for the components in `AbstractCompressibleEulerMulticomponentEquations`. 
+for the components in `AbstractCompressibleEulerMulticomponentEquations`.
 In particular, not the components themselves are returned.
 """
 @inline eachcomponent(equations::AbstractCompressibleEulerMulticomponentEquations) = Base.OneTo(ncomponents(equations))
@@ -350,7 +351,7 @@ include("ideal_glm_mhd_multicomponent_2d.jl")
     eachcomponent(equations::AbstractIdealGlmMhdMulticomponentEquations)
 
 Return an iterator over the indices that specify the location in relevant data structures
-for the components in `AbstractIdealGlmMhdMulticomponentEquations`. 
+for the components in `AbstractIdealGlmMhdMulticomponentEquations`.
 In particular, not the components themselves are returned.
 """
 @inline eachcomponent(equations::AbstractIdealGlmMhdMulticomponentEquations) = Base.OneTo(ncomponents(equations))

src/solvers/dgmulti/dg.jl

@@ -428,24 +428,27 @@ end
 
 # do nothing for periodic (default) boundary conditions
 calc_boundary_flux!(cache, t, boundary_conditions::BoundaryConditionPeriodic,
-                    mesh, equations, dg::DGMulti) = nothing
+                    mesh, have_nonconservative_terms, equations, dg::DGMulti) = nothing
 
 # "lispy tuple programming" instead of for loop for type stability
-function calc_boundary_flux!(cache, t, boundary_conditions, mesh, equations, dg::DGMulti)
+function calc_boundary_flux!(cache, t, boundary_conditions, mesh,
+                             have_nonconservative_terms, equations, dg::DGMulti)
+
   # peel off first boundary condition
   calc_single_boundary_flux!(cache, t, first(boundary_conditions), first(keys(boundary_conditions)),
-                 mesh, equations, dg)
+                             mesh, have_nonconservative_terms, equations, dg)
 
   # recurse on the remainder of the boundary conditions
-  calc_boundary_flux!(cache, t, Base.tail(boundary_conditions), mesh, equations, dg)
+  calc_boundary_flux!(cache, t, Base.tail(boundary_conditions),
+                      mesh, have_nonconservative_terms, equations, dg)
 end
 
 # terminate recursion
 calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple{(),Tuple{}},
-                    mesh, equations, dg::DGMulti) = nothing
+                    mesh, have_nonconservative_terms, equations, dg::DGMulti) = nothing
 
-function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
-                                    mesh, equations, dg::DGMulti{NDIMS}) where {NDIMS}
+function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, mesh,
+                                    have_nonconservative_terms::False, equations, dg::DGMulti{NDIMS}) where {NDIMS}
 
   rd = dg.basis
   md = mesh.md
@@ -455,8 +458,9 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
 
   # reshape face/normal arrays to have size = (num_points_on_face, num_faces_total).
   # mesh.boundary_faces indexes into the columns of these face-reshaped arrays.
-  num_pts_per_face = rd.Nfq ÷ rd.Nfaces
-  num_faces_total = rd.Nfaces * md.num_elements
+  num_faces = StartUpDG.num_faces(rd.element_type)
+  num_pts_per_face = rd.Nfq ÷ num_faces
+  num_faces_total = num_faces * md.num_elements
 
   # This function was originally defined as
   # `reshape_by_face(u) = reshape(view(u, :), num_pts_per_face, num_faces_total)`.
@@ -485,6 +489,58 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
   # However, we don't have to re-reshape, since cache.flux_face_values still retains its original shape.
 end
 
+function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, mesh,
+                                    have_nonconservative_terms::True, equations, dg::DGMulti{NDIMS}) where {NDIMS}
+
+  rd = dg.basis
+  md = mesh.md
+  surface_flux, nonconservative_flux = dg.surface_integral.surface_flux
+
+  # reshape face/normal arrays to have size = (num_points_on_face, num_faces_total).
+  # mesh.boundary_faces indexes into the columns of these face-reshaped arrays.
+  num_pts_per_face = rd.Nfq ÷ StartUpDG.num_faces(rd.element_type)
+  num_faces_total = StartUpDG.num_faces(rd.element_type) * md.num_elements
+
+  # This function was originally defined as
+  # `reshape_by_face(u) = reshape(view(u, :), num_pts_per_face, num_faces_total)`.
+  # This results in allocations due to https://github.com/JuliaLang/julia/issues/36313.
+  # To avoid allocations, we use Tim Holy's suggestion:
+  # https://github.com/JuliaLang/julia/issues/36313#issuecomment-782336300.
+  reshape_by_face(u) = Base.ReshapedArray(u, (num_pts_per_face, num_faces_total), ())
+
+  u_face_values = reshape_by_face(cache.u_face_values)
+  flux_face_values = reshape_by_face(cache.flux_face_values)
+  Jf = reshape_by_face(md.Jf)
+  nxyzJ, xyzf = reshape_by_face.(md.nxyzJ), reshape_by_face.(md.xyzf) # broadcast over nxyzJ::NTuple{NDIMS,Matrix}
+
+  # loop through boundary faces, which correspond to columns of reshaped u_face_values, ...
+  for f in mesh.boundary_faces[boundary_key]
+    for i in Base.OneTo(num_pts_per_face)
+      face_normal = SVector{NDIMS}(getindex.(nxyzJ, i, f)) / Jf[i,f]
+      face_coordinates = SVector{NDIMS}(getindex.(xyzf, i, f))
+
+      # Compute conservative and non-conservative fluxes separately.
+      # This imposes boundary conditions on the conservative part of the flux.
+      cons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t,
+                                                  surface_flux, equations)
+
+      # Compute pointwise nonconservative numerical flux at the boundary.
+      # In general, nonconservative fluxes can depend on both the contravariant
+      # vectors (normal direction) at the current node and the averaged ones.
+      # However, there is only one `face_normal` at boundaries, which we pass in twice.
+      # Note: This does not set any type of boundary condition for the nonconservative term
+      noncons_flux_at_face_node = nonconservative_flux(u_face_values[i,f], u_face_values[i,f],
+                                                       face_normal, face_normal, equations)
+
+      flux_face_values[i,f] = (cons_flux_at_face_node + 0.5 * noncons_flux_at_face_node) * Jf[i,f]
+
+    end
+  end
+
+  # Note: modifying the values of the reshaped array modifies the values of cache.flux_face_values.
+  # However, we don't have to re-reshape, since cache.flux_face_values still retains its original shape.
+end
+
 
 # inverts Jacobian and scales by -1.0
 function invert_jacobian!(du, mesh::DGMultiMesh, equations, dg::DGMulti, cache; scaling=-1)
@@ -568,7 +624,8 @@ function rhs!(du, u, t, mesh, equations,
     have_nonconservative_terms(equations), equations, dg)
 
   @trixi_timeit timer() "boundary flux" calc_boundary_flux!(
-    cache, t, boundary_conditions, mesh, equations, dg)
+    cache, t, boundary_conditions, mesh,
+    have_nonconservative_terms(equations), equations, dg)
 
   @trixi_timeit timer() "surface integral" calc_surface_integral!(
     du, u, dg.surface_integral, mesh, equations, dg, cache)

src/solvers/dgmulti/flux_differencing.jl

@@ -594,8 +594,8 @@ function rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions:
                                                               have_nonconservative_terms(equations),
                                                               equations, dg)
 
-  @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions,
-                                                            mesh, equations, dg)
+  @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, mesh,
+                                                            have_nonconservative_terms(equations), equations, dg)
 
   @trixi_timeit timer() "surface integral" calc_surface_integral!(du, u, dg.surface_integral,
                                                                   mesh, equations, dg, cache)
@@ -630,7 +630,8 @@ function rhs!(du, u, t, mesh, equations,
     have_nonconservative_terms(equations), equations, dg)
 
   @trixi_timeit timer() "boundary flux" calc_boundary_flux!(
-    cache, t, boundary_conditions, mesh, equations, dg)
+    cache, t, boundary_conditions, mesh,
+    have_nonconservative_terms(equations), equations, dg)
 
   @trixi_timeit timer() "surface integral" calc_surface_integral!(
     du, u, dg.surface_integral, mesh, equations, dg, cache)

src/solvers/dgmulti/flux_differencing_gauss_sbp.jl

@@ -61,7 +61,8 @@ function TensorProductGaussFaceOperator(operator::AbstractGaussOperator,
   nnodes_1d = length(rq1D)
 
   # Permutation of indices in a tensor product form
-  indices = reshape(1:length(rd.rf), nnodes_1d, rd.Nfaces)
+  num_faces = StartUpDG.num_faces(rd.element_type)
+  indices = reshape(1:length(rd.rf), nnodes_1d, num_faces)
   face_indices_tensor_product = zeros(Int, 2, nnodes_1d, ndims(rd.element_type))
   for i in 1:nnodes_1d # loop over nodes in one face
     face_indices_tensor_product[:, i, 1] .= indices[i, 1:2]
@@ -76,7 +77,7 @@ function TensorProductGaussFaceOperator(operator::AbstractGaussOperator,
   return TensorProductGaussFaceOperator{2, T_op, Tm, Tw, Tf, Ti}(interp_matrix_gauss_to_face_1d,
                                                                  inv.(wq1D), rd.wf,
                                                                  face_indices_tensor_product,
-                                                                 nnodes_1d, rd.Nfaces)
+                                                                 nnodes_1d, num_faces)
 end
 
 # constructor for a 3D operator
@@ -90,7 +91,8 @@ function TensorProductGaussFaceOperator(operator::AbstractGaussOperator,
   nnodes_1d = length(rq1D)
 
   # Permutation of indices in a tensor product form
-  indices = reshape(1:length(rd.rf), nnodes_1d, nnodes_1d, rd.Nfaces)
+  num_faces = StartUpDG.num_faces(rd.element_type)
+  indices = reshape(1:length(rd.rf), nnodes_1d, nnodes_1d, num_faces)
   face_indices_tensor_product = zeros(Int, 2, nnodes_1d, nnodes_1d, ndims(rd.element_type))
   for j in 1:nnodes_1d, i in 1:nnodes_1d # loop over nodes in one face
     face_indices_tensor_product[:, i, j, 1] .= indices[i, j, 1:2]
@@ -106,7 +108,7 @@ function TensorProductGaussFaceOperator(operator::AbstractGaussOperator,
   return TensorProductGaussFaceOperator{3, T_op, Tm, Tw, Tf, Ti}(interp_matrix_gauss_to_face_1d,
                                                                  inv.(wq1D), rd.wf,
                                                                  face_indices_tensor_product,
-                                                                 nnodes_1d, rd.Nfaces)
+                                                                 nnodes_1d, num_faces)
 end
 
 # specialize behavior of `mul_by!(A)` where `A isa TensorProductGaussFaceOperator)`
@@ -507,8 +509,8 @@ function rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions:
                                                               have_nonconservative_terms(equations),
                                                               equations, dg)
 
-  @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions,
-                                                            mesh, equations, dg)
+  @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, mesh,
+                                                            have_nonconservative_terms(equations), equations, dg)
 
   # `du` is stored at Gauss nodes here
   @trixi_timeit timer() "surface integral" calc_surface_integral!(du, u, dg.surface_integral,

src/solvers/dgsem_tree/dg_2d.jl

@@ -702,7 +702,7 @@ end
 
 function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, t,
                                           boundary_condition, nonconservative_terms::True, equations,
-                                          surface_integral ,dg::DG, cache,
+                                          surface_integral, dg::DG, cache,
                                           direction, first_boundary, last_boundary)
   surface_flux, nonconservative_flux = surface_integral.surface_flux
   @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries

test/test_dgmulti_2d.jl

@@ -292,6 +292,14 @@ isdir(outdir) && rm(outdir, recursive=true)
     )
   end
 
+  @trixi_testset "elixir_mhd_reflective_wall.jl (Quad)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_reflective_wall.jl"),
+      cells_per_dimension = 4,
+      l2 = [0.0036019536614619687, 0.001734097206958611, 0.008375221008997178, 0.0, 0.028596796602124414, 0.0018573693138866614, 0.0020807798141551166, 0.0, 5.301188920230166e-5],
+      linf = [0.01692601228199253, 0.009369662298436778, 0.04145169295835428, 0.0, 0.11569908670112738, 0.00984964453299233, 0.01141708032148614, 0.0, 0.0002992631411931389]
+    )
+  end
+
   @trixi_testset "elixir_shallowwater_source_terms.jl (Quad, SBP)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
       cells_per_dimension = 8, element_type = Quad(), approximation_type = SBP(),