Add flexibility to nonconservative BCs (#2200)

* adding flexibility to noncons BCs

* minor fix

* Dirichlet BC for noncons

* formatting

* Noncons Dirichlet BC

* test_type.jl fix

* dg multi fix

* moving 0.5 factor internally

* Update src/solvers/dgsem_p4est/dg_2d.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/shallow_water_2d.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/shallow_water_2d.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/equations.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update examples/dgmulti_2d/elixir_mhd_reflective_wall.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/shallow_water_2d.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/shallow_water_2d.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/shallow_water_2d.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/equations.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/equations.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/equations.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/shallow_water_1d.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/shallow_water_1d.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update src/equations/shallow_water_1d.jl

Co-authored-by: Daniel Doehring <[email protected]>

* minor fixes

* fix test type stability

* fix test type stability

* Update src/equations/shallow_water_1d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/equations/shallow_water_2d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/basic_types.jl

Co-authored-by: Patrick Ersing <[email protected]>

* Update src/basic_types.jl

Co-authored-by: Patrick Ersing <[email protected]>

* Update src/solvers/dgsem_p4est/dg_2d.jl

Co-authored-by: Patrick Ersing <[email protected]>

* Update src/solvers/dgsem_unstructured/dg_2d.jl

Co-authored-by: Patrick Ersing <[email protected]>

* Update src/solvers/dgsem_unstructured/dg_2d.jl

Co-authored-by: Patrick Ersing <[email protected]>

* Update src/basic_types.jl

Co-authored-by: Daniel Doehring <[email protected]>

* formatting

* Update src/basic_types.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Update test/test_type.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update test/test_type.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/equations/equations.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/equations/equations.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* add docs

---------

Co-authored-by: Daniel Doehring <[email protected]>
Co-authored-by: Andrew Winters <[email protected]>
Co-authored-by: Patrick Ersing <[email protected]>
Co-authored-by: Hendrik Ranocha <[email protected]>

Author: 5 people

NEWS.md

@@ -23,6 +23,10 @@ for human readability.
   Previously, only the first 8 digits were printed to file.
   Furthermore, the names of the printed fields are now only separated by a single white space, 
   in contrast to before where this were multiple, depending on the actual name of the printed data.
+- The boundary conditions for non-conservative equations can now be defined separately from the conservative part.
+  The `surface_flux_functions` tuple is now passed directly to the boundary condition call,
+  returning a tuple with boundary condition values for both the conservative and non-conservative parts ([#2200]).
+
 
 #### Deprecated
 

docs/literate/src/files/adding_nonconservative_equation.jl

@@ -177,6 +177,9 @@ error_1 / error_2
 # As expected, the new error is roughly reduced by a factor of 16, corresponding
 # to an experimental order of convergence of 4 (for polynomials of degree 3).
 
+# For non-trivial boundary conditions involving non-conservative terms, 
+# please refer to the section on [Other available example elixirs with non-trivial BC](https://trixi-framework.github.io/Trixi.jl/stable/tutorials/non_periodic_boundaries/#Other-available-example-elixirs-with-non-trivial-BC).
+
 # ## Summary of the code
 
 # Here is the complete code that we used (without the callbacks since these

docs/literate/src/files/non_periodic_boundaries.jl

@@ -158,6 +158,11 @@ end
 # ```
 # Source: [`Video`](https://www.youtube.com/watch?v=w0A9X38cSe4) on Trixi.jl's YouTube channel [`Trixi Framework`](https://www.youtube.com/watch?v=WElqqdMhY4A)
 
+# Furthermore, Trixi.jl also handles equations that include non-conservative terms. 
+# For such equations, the tuple of conservative and non-conservative surfaces fluxes is passed to the boundary condition, 
+# which then returns a tuple containing the boundary condition values for both the conservative and non-conservative terms.
+# For instance, a 2D ideal compressible GLM-MHD setup with reflective walls can be found in the elixir ['elixir_mhd_reflective_wall.jl](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/dgmulti_2d/elixir_mhd_reflective_wall.jl).   
+
 # ## Package versions
 
 # These results were obtained using the following versions.

examples/dgmulti_2d/elixir_mhd_reflective_wall.jl

@@ -48,9 +48,9 @@ mesh = DGMultiMesh(solver, cells_per_dimension; periodicity = (false, false),
 # 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,
+                                               x, t, surface_flux_functions,
                                                equations::IdealGlmMhdEquations2D)
-
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
     # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later
     norm_ = norm(normal_direction)
     normal = normal_direction / norm_
@@ -62,8 +62,10 @@ function boundary_condition_velocity_slip_wall(u_inner, normal_direction::Abstra
     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_
+    flux = surface_flux_function(u_inner, u_mirror, normal, equations) * norm_
+    noncons_flux = nonconservative_flux_function(u_inner, u_mirror, normal, equations) *
+                   norm_
+    return flux, noncons_flux
 end
 
 boundary_conditions = (; x_neg = boundary_condition_velocity_slip_wall,

src/basic_types.jl

@@ -71,14 +71,39 @@ struct BoundaryConditionDoNothing end
                                                 orientation_or_normal_direction,
                                                 direction::Integer, x, t, surface_flux,
                                                 equations)
-    return surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations)
+    return flux(u_inner, orientation_or_normal_direction, equations)
+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_functions::Tuple,
+                                                equations)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+    return surface_flux_function(u_inner, u_inner,
+                                 orientation_or_normal_direction, equations),
+           nonconservative_flux_function(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 surface_flux(u_inner, u_inner, outward_direction, equations)
+    return flux(u_inner, outward_direction, equations)
+end
+
+# Version for equations involving nonconservative fluxes
+@inline function (::BoundaryConditionDoNothing)(u_inner,
+                                                outward_direction::AbstractVector,
+                                                x, t, surface_flux_functions::Tuple,
+                                                equations)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
+    return surface_flux_function(u_inner, u_inner, outward_direction, equations),
+           nonconservative_flux_function(u_inner, u_inner, outward_direction,
+                                         equations)
 end
 
 # This version can be called by parabolic solvers

src/equations/equations.jl

@@ -181,6 +181,36 @@ end
     return flux
 end
 
+# Dirichlet-type boundary condition for use with TreeMesh or StructuredMesh
+# passing a tuple of surface flux functions for nonconservative terms
+@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
+                                                                  orientation_or_normal,
+                                                                  direction,
+                                                                  x, t,
+                                                                  surface_flux_functions::Tuple,
+                                                                  equations)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
+    u_boundary = boundary_condition.boundary_value_function(x, t, equations)
+
+    # Calculate boundary flux
+    if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
+        flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal,
+                                     equations)
+        noncons_flux = nonconservative_flux_function(u_inner, u_boundary,
+                                                     orientation_or_normal,
+                                                     equations)
+    else # u_boundary is "left" of boundary, u_inner is "right" of boundary
+        flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal,
+                                     equations)
+        noncons_flux = nonconservative_flux_function(u_boundary, u_inner,
+                                                     orientation_or_normal,
+                                                     equations)
+    end
+
+    return flux, noncons_flux
+end
+
 # Dirichlet-type boundary condition for use with UnstructuredMesh2D
 # Note: For unstructured we lose the concept of an "absolute direction"
 @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
@@ -197,6 +227,26 @@ end
     return flux
 end
 
+# Dirichlet-type boundary condition for equations with non-conservative terms for use with UnstructuredMesh2D
+# passing a tuple of surface flux functions for nonconservative terms
+# Note: For unstructured we lose the concept of an "absolute direction"
+@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
+                                                                  normal_direction::AbstractVector,
+                                                                  x, t,
+                                                                  surface_flux_functions::Tuple,
+                                                                  equations)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
+    # get the external value of the solution
+    u_boundary = boundary_condition.boundary_value_function(x, t, equations)
+
+    # Calculate boundary flux
+    flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+    noncons_flux = nonconservative_flux_function(u_inner, u_boundary, normal_direction,
+                                                 equations)
+    return flux, noncons_flux
+end
+
 # operator types used for dispatch on parabolic boundary fluxes
 struct Gradient end
 struct Divergence end

src/equations/shallow_water_1d.jl

@@ -160,9 +160,12 @@ For details see Section 9.2.5 of the book:
 """
 @inline function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction,
                                               x, t,
-                                              surface_flux_function,
+                                              surface_flux_functions,
                                               equations::ShallowWaterEquations1D)
 
+    # The boundary conditions for the non-conservative term are identically 0 here.
+    # Bottom topography is assumed to be continuous at the boundary.
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
     # create the "external" boundary solution state
     u_boundary = SVector(u_inner[1],
                          -u_inner[2],
@@ -172,12 +175,18 @@ For details see Section 9.2.5 of the book:
     if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
         flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal,
                                      equations)
+        noncons_flux = nonconservative_flux_function(u_inner, u_boundary,
+                                                     orientation_or_normal,
+                                                     equations)
     else # u_boundary is "left" of boundary, u_inner is "right" of boundary
         flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal,
                                      equations)
+        noncons_flux = nonconservative_flux_function(u_boundary, u_inner,
+                                                     orientation_or_normal,
+                                                     equations)
     end
 
-    return flux
+    return flux, noncons_flux
 end
 
 # Calculate 1D flux for a single point

src/equations/shallow_water_2d.jl

@@ -180,8 +180,10 @@ For details see Section 9.2.5 of the book:
 """
 @inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector,
                                               x, t,
-                                              surface_flux_function,
+                                              surface_flux_functions,
                                               equations::ShallowWaterEquations2D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
     # normalize the outward pointing direction
     normal = normal_direction / norm(normal_direction)
 
@@ -196,8 +198,10 @@ For details see Section 9.2.5 of the book:
 
     # calculate the boundary flux
     flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+    noncons_flux = nonconservative_flux_function(u_inner, u_boundary, normal_direction,
+                                                 equations)
 
-    return flux
+    return flux, noncons_flux
 end
 
 """
@@ -208,8 +212,11 @@ Should be used together with [`TreeMesh`](@ref).
 """
 @inline function boundary_condition_slip_wall(u_inner, orientation,
                                               direction, x, t,
-                                              surface_flux_function,
+                                              surface_flux_functions,
                                               equations::ShallowWaterEquations2D)
+    # The boundary conditions for the non-conservative term are identically 0 here.
+    # Bottom topography is assumed to be continuous at the boundary.
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
     ## get the appropriate normal vector from the orientation
     if orientation == 1
         u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4])
@@ -220,11 +227,15 @@ Should be used together with [`TreeMesh`](@ref).
     # Calculate boundary flux
     if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
         flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
+        noncons_flux = nonconservative_flux_function(u_inner, u_boundary, orientation,
+                                                     equations)
     else # u_boundary is "left" of boundary, u_inner is "right" of boundary
         flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
+        noncons_flux = nonconservative_flux_function(u_boundary, u_inner, orientation,
+                                                     equations)
     end
 
-    return flux
+    return flux, noncons_flux
 end
 
 # Calculate 1D flux for a single point

src/solvers/dgmulti/dg.jl

@@ -523,7 +523,6 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
                                     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.
@@ -550,21 +549,16 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
 
             # 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.
-            noncons_flux_at_face_node = boundary_condition(u_face_values[i, f],
-                                                           face_normal,
-                                                           face_coordinates,
-                                                           t,
-                                                           nonconservative_flux,
-                                                           equations)
+            cons_flux_at_face_node, noncons_flux_at_face_node = boundary_condition(u_face_values[i,
+                                                                                                 f],
+                                                                                   face_normal,
+                                                                                   face_coordinates,
+                                                                                   t,
+                                                                                   dg.surface_integral.surface_flux,
+                                                                                   equations)
 
             flux_face_values[i, f] = (cons_flux_at_face_node +
-                                      0.5 * noncons_flux_at_face_node) * Jf[i, f]
+                                      0.5f0 * noncons_flux_at_face_node) * Jf[i, f]
         end
     end
 

src/solvers/dgsem_p4est/dg_2d.jl

@@ -349,7 +349,6 @@ end
                                      boundary_index)
     @unpack boundaries = cache
     @unpack node_coordinates, contravariant_vectors = cache.elements
-    surface_flux, nonconservative_flux = surface_integral.surface_flux
 
     # Extract solution data from boundary container
     u_inner = get_node_vars(boundaries.u, equations, dg, node_index, boundary_index)
@@ -362,22 +361,19 @@ end
     x = get_node_coords(node_coordinates, equations, dg, i_index, j_index,
                         element_index)
 
-    # Call pointwise numerical flux function for the conservative part
+    # Call pointwise numerical flux functions for the conservative and nonconservative part
     # in the normal direction on the boundary
-    flux_ = boundary_condition(u_inner, normal_direction, x, t, surface_flux, equations)
-
-    # Compute pointwise nonconservative numerical flux at the boundary.
-    noncons_ = boundary_condition(u_inner, normal_direction, x, t, nonconservative_flux,
-                                  equations)
+    flux, noncons_flux = boundary_condition(u_inner, normal_direction, x, t,
+                                            surface_integral.surface_flux, equations)
 
     # Copy flux to element storage in the correct orientation
     for v in eachvariable(equations)
         # Note the factor 0.5 necessary for the nonconservative fluxes based on
         # the interpretation of global SBP operators coupled discontinuously via
         # central fluxes/SATs
-        surface_flux_values[v, node_index, direction_index, element_index] = flux_[v] +
+        surface_flux_values[v, node_index, direction_index, element_index] = flux[v] +
                                                                              0.5f0 *
-                                                                             noncons_[v]
+                                                                             noncons_flux[v]
     end
 end