Apply suggestions from code review

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

examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl

@@ -1,3 +1,5 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -20,6 +22,7 @@ boundary_conditions = Dict(:outerCircle => boundary_condition_free_stream,
 ###############################################################################
 # Get the Upwind FDSBP approximation space
 
+# TODO: FDSBP
 # Note, one must set `xmin=-1` and `xmax=1` due to the reuse
 # of interpolation routines from `calc_node_coordinates!` to create
 # the physical coordinates in the mappings.
@@ -37,8 +40,9 @@ solver = FDSBP(D_upw,
 ###############################################################################
 # Get the curved quad mesh from a file (downloads the file if not available locally)
 
-# Mesh with second order boundary polynomials requires an upwind SBP operator
-# with (at least) 4th order boundary closure to guarantee the approximation is free-stream
+# Mesh with second-order boundary polynomials requires an upwind SBP operator
+# with (at least) 4th order boundary closure to guarantee the approximation is
+# free-stream preserving
 mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/ec9a345f09199ebe471d35d5c1e4e08f/raw/15975943d8642e42f8292235314b6f1b30aa860d/mesh_inner_outer_boundaries.mesh",
                            joinpath(@__DIR__, "mesh_inner_outer_boundaries.mesh"))
 

examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl

@@ -1,3 +1,5 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
 
 using OrdinaryDiffEq
 using Trixi
@@ -21,6 +23,7 @@ boundary_conditions = Dict(:Top => boundary_condition_eoc,
 ###############################################################################
 # Get the Upwind FDSBP approximation space
 
+# TODO: FDSBP
 # Note, one must set `xmin=-1` and `xmax=1` due to the reuse
 # of interpolation routines from `calc_node_coordinates!` to create
 # the physical coordinates in the mappings.
@@ -38,15 +41,16 @@ solver = FDSBP(D_upw,
 ###############################################################################
 # Get the curved quad mesh from a file (downloads the file if not available locally)
 
-# Mesh with first order boundary polynomials requires an upwind SBP operator
-# with (at least) 2nd order boundary closure to guarantee the approximation is free-stream
+# Mesh with first-order boundary polynomials requires an upwind SBP operator
+# with (at least) 2nd order boundary closure to guarantee the approximation is
+# free-stream preserving
 mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/a4f4743008bf3233957a9ea6ac7a62e0/raw/8b36cc6649153fe0a5723b200368a210a1d74eaf/mesh_refined_box.mesh",
                            joinpath(@__DIR__, "mesh_refined_box.mesh"))
 
 mesh = UnstructuredMesh2D(mesh_file)
 
 ###############################################################################
-# create the semi discretization object
+# create the semidiscretization object
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     source_terms = source_term,

src/solvers/fdsbp_unstructured/fdsbp_2d.jl

@@ -125,16 +125,19 @@ function calc_volume_integral!(du, u,
     # the fluxes line-by-line and add them to `du` for each element.
     @threaded for element in eachelement(dg, cache)
         # f_minus_plus_element wraps the storage provided by f_minus_element and
-        # f_plus_element such that we can use a single plain broadcasting below.
-        # f_minus_element and f_plus_element are updated in broadcasting calls
-        # of the form `@. f_minus_plus_element = ...`.
+        # f_plus_element such that we can use a single assignment below.
+        # f_minus_element and f_plus_element are updated whenever we update
+        # `f_minus_plus_element[i, j] = ...` below.
         f_minus_plus_element = f_minus_plus_threaded[Threads.threadid()]
         f_minus_element = f_minus_threaded[Threads.threadid()]
         f_plus_element = f_plus_threaded[Threads.threadid()]
         u_element = view(u_vectors, :, :, element)
 
         # x direction
-        # @. f_minus_plus_element = splitting(u_element, 1, equations)
+        # We use flux vector splittings in the directions of the contravariant
+        # basis vectors. Thus, we do not use a broadcasting operation like
+        #   @. f_minus_plus_element = splitting(u_element, 1, equations)
+        # in the Cartesian case but loop over all nodes.
         for j in eachnode(dg), i in eachnode(dg)
             # contravariant vectors computed with central D matrix
             Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element)