specialize calc_boundary_flux! for nonconservative terms for DGMulti

examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl

@@ -0,0 +1,99 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# 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 = SBP(),
+                 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)
+
+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, .075)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+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,
+                        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

@@ -77,19 +77,22 @@ struct BoundaryConditionDoNothing end
 @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)
+  # TODO: should this be switched to `surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations)` for consistency?
 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)
+
+  # this should reduce to `flux(u_inner, outward_direction, equations)` for a consistent symmetric flux.
+  return surface_flux(u_inner, 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/ideal_glm_mhd_2d.jl

@@ -238,6 +238,11 @@ terms.
   return f
 end
 
+# TODO: is there a way to do this without manually dispatching for every non-conservative flux? We use this for the evaluation of boundary conditions.
+# if only one normal is specified, use that for `normal_direction_average`
+@inline flux_nonconservative_powell(u_ll, u_rr, normal_direction::AbstractVector, equations::IdealGlmMhdEquations2D) =
+  flux_nonconservative_powell(u_ll, u_rr, normal_direction, normal_direction, equations)
+
 @inline function flux_nonconservative_powell(u_ll, u_rr,
                                              normal_direction_ll::AbstractVector,
                                              normal_direction_average::AbstractVector,

src/solvers/dgmulti/dg.jl

@@ -428,23 +428,27 @@ end
 
 # do nothing for periodic (default) boundary conditions
 calc_boundary_flux!(cache, t, boundary_conditions::BoundaryConditionPeriodic,
+                    have_nonconservative_terms,
                     mesh, 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, have_nonconservative_terms,
+                             mesh, equations, dg::DGMulti)
   # peel off first boundary condition
   calc_single_boundary_flux!(cache, t, first(boundary_conditions), first(keys(boundary_conditions)),
-                 mesh, equations, dg)
+                             have_nonconservative_terms, mesh, 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),
+                      have_nonconservative_terms, mesh, equations, dg)
 end
 
 # terminate recursion
 calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple{(),Tuple{}},
-                    mesh, equations, dg::DGMulti) = nothing
+                    have_nonconservative_terms, mesh, equations, dg::DGMulti) = nothing
 
 function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
+                                    have_nonconservative_terms::False,
                                     mesh, equations, dg::DGMulti{NDIMS}) where {NDIMS}
 
   rd = dg.basis
@@ -485,6 +489,55 @@ 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,
+                                    have_nonconservative_terms::True,
+                                    mesh, equations, dg::DGMulti{NDIMS}) where {NDIMS}
+
+  rd = dg.basis
+  md = mesh.md
+  (; u_face_values, flux_face_values) = cache
+  (; xyzf, nxyzJ, Jf) = 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 ÷ rd.Nfaces
+  num_faces_total = rd.Nfaces * 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(u_face_values)
+  flux_face_values = reshape_by_face(flux_face_values)
+  Jf = reshape_by_face(Jf)
+  nxyzJ, xyzf = reshape_by_face.(nxyzJ), reshape_by_face.(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 parts separately
+      cons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t,
+                                                  surface_flux, equations)
+
+      noncons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t,
+                                                     nonconservative_flux, 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 +621,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,
+    have_nonconservative_terms(equations), mesh, 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

@@ -595,6 +595,7 @@ function rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions:
                                                               equations, dg)
 
   @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions,
+                                                            have_nonconservative_terms(equations),
                                                             mesh, equations, dg)
 
   @trixi_timeit timer() "surface integral" calc_surface_integral!(du, u, dg.surface_integral,
@@ -630,7 +631,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, have_nonconservative_terms(equations),
+    mesh, 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

@@ -508,6 +508,7 @@ function rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions:
                                                               equations, dg)
 
   @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions,
+                                                            have_nonconservative_terms(equations),
                                                             mesh, equations, dg)
 
   # `du` is stored at Gauss nodes here

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_BCs.jl (Quad)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_reflective_BCs.jl"),
+      cells_per_dimension = 4, element_type = Quad(), approximation_type = GaussSBP(),
+      l2 = [0.0036019562526885497, 0.0017340971255526598, 0.008375221676923516, 0.0, 0.028596802653952254, 0.0018573697892179616, 0.002080779894047767, 0.0, 5.3012597624420295e-5],
+      linf = [0.01692598382370969, 0.00936965952976002, 0.04145170727841342, 0.0, 0.11569901084179435, 0.009849648257862952, 0.011417088537134745, 0.0, 0.0002992621756978016]
+    )
+  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(),