From 77cd96e681ff68a12c6d259fc806f0503d197c32 Mon Sep 17 00:00:00 2001
From: Andrew Winters <andrew.ross.winters@liu.se>
Date: Mon, 12 Dec 2022 19:50:50 +0100
Subject: [PATCH 01/15] add FDSBP for curvilinear meshes

---
 src/solvers/dg.jl                             |   1 +
 .../dgsem_unstructured/containers_2d.jl       |   4 +-
 src/solvers/fdsbp_tree/fdsbp_2d.jl            |   4 +-
 .../fdsbp_unstructured/containers_2d.jl       | 124 +++++++++++
 src/solvers/fdsbp_unstructured/fdsbp.jl       |  16 ++
 src/solvers/fdsbp_unstructured/fdsbp_2d.jl    | 199 ++++++++++++++++++
 6 files changed, 344 insertions(+), 4 deletions(-)
 create mode 100644 src/solvers/fdsbp_unstructured/containers_2d.jl
 create mode 100644 src/solvers/fdsbp_unstructured/fdsbp.jl
 create mode 100644 src/solvers/fdsbp_unstructured/fdsbp_2d.jl

diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
index d5fa15d3e1c..487e7f1a476 100644
--- a/src/solvers/dg.jl
+++ b/src/solvers/dg.jl
@@ -584,6 +584,7 @@ include("dgsem_p4est/dg.jl")
 # and boundary conditions weakly. Thus, these methods can re-use a lot of
 # functionality implemented for DGSEM.
 include("fdsbp_tree/fdsbp.jl")
+include("fdsbp_unstructured/fdsbp.jl")
 
 
 end # @muladd
diff --git a/src/solvers/dgsem_unstructured/containers_2d.jl b/src/solvers/dgsem_unstructured/containers_2d.jl
index 986ba342e1d..3417adc30d2 100644
--- a/src/solvers/dgsem_unstructured/containers_2d.jl
+++ b/src/solvers/dgsem_unstructured/containers_2d.jl
@@ -62,7 +62,7 @@ function init_elements(mesh::UnstructuredMesh2D, equations, basis, RealT, uEltyp
 end
 
 
-function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis)
+function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis::LobattoLegendreBasis)
   four_corners = zeros(eltype(mesh.corners), 4, 2)
 
   # loop through elements and call the correct constructor based on whether the element is curved
@@ -80,7 +80,7 @@ function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis)
 end
 
 
-# initialize all the values in the container of a general element (either straight sided or curved)
+# initialize all the values in the container of a general DG element (either straight sided or curved)
 function init_element!(elements, element, nodes, corners_or_surface_curves)
 
   calc_node_coordinates!(elements.node_coordinates, element, nodes, corners_or_surface_curves)
diff --git a/src/solvers/fdsbp_tree/fdsbp_2d.jl b/src/solvers/fdsbp_tree/fdsbp_2d.jl
index 83d75ff3e81..e4c4186b884 100644
--- a/src/solvers/fdsbp_tree/fdsbp_2d.jl
+++ b/src/solvers/fdsbp_tree/fdsbp_2d.jl
@@ -9,7 +9,7 @@
 
 
 # 2D caches
-function create_cache(mesh::TreeMesh{2}, equations,
+function create_cache(mesh::Union{TreeMesh{2}, UnstructuredMesh2D}, equations,
                       volume_integral::VolumeIntegralStrongForm, dg, uEltype)
 
   prototype = Array{SVector{nvariables(equations), uEltype}, ndims(mesh)}(
@@ -19,7 +19,7 @@ function create_cache(mesh::TreeMesh{2}, equations,
   return (; f_threaded,)
 end
 
-function create_cache(mesh::TreeMesh{2}, equations,
+function create_cache(mesh::Union{TreeMesh{2}, UnstructuredMesh2D}, equations,
                       volume_integral::VolumeIntegralUpwind, dg, uEltype)
 
   u_node = SVector{nvariables(equations), uEltype}(ntuple(_ -> zero(uEltype), Val{nvariables(equations)}()))
diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl
new file mode 100644
index 00000000000..4e0945dbd6d
--- /dev/null
+++ b/src/solvers/fdsbp_unstructured/containers_2d.jl
@@ -0,0 +1,124 @@
+# !!! warning "Experimental implementation (curvilinear FDSBP)"
+#     This is an experimental feature and may change in future releases.
+
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+
+# Alternative constructor to reuse the element container from `UnstructuredMesh2D`
+# OBS! New element constructors require the derivative operator to compute the metric terms.
+function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis::AbstractDerivativeOperator)
+  four_corners = zeros(eltype(mesh.corners), 4, 2)
+
+  # loop through elements and call the correct constructor based on whether the element is curved
+  for element in eachelement(elements)
+    if mesh.element_is_curved[element]
+      init_element!(elements, element, basis, view(mesh.surface_curves, :, element))
+    else # straight sided element
+      for i in 1:4, j in 1:2
+        # pull the (x,y) values of these corners out of the global corners array
+        four_corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]]
+      end
+      init_element!(elements, element, basis, four_corners)
+    end
+  end
+end
+
+
+# initialize all the values in the container of a general FD block (either straight sided or curved)
+# OBS! Requires the SBP derivative matrix in order to compute metric terms that are free-stream preserving
+function init_element!(elements, element, basis::AbstractDerivativeOperator, corners_or_surface_curves)
+
+  calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis), corners_or_surface_curves)
+
+  calc_metric_terms!(elements.jacobian_matrix, element, basis, elements.node_coordinates)
+
+  calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix)
+
+  calc_contravariant_vectors!(elements.contravariant_vectors, element, elements.jacobian_matrix)
+
+  calc_normal_directions!(elements.normal_directions, element, elements.jacobian_matrix)
+
+  return elements
+end
+
+
+# construct the metric terms for a FDSBP element "block". Directly use the derivative matrix
+# applied to the node coordinates.
+# TODO: FD performance; this is slow and allocates. Rewrite to use `mul!`
+# TODO: FD; How to make this work for the upwind solver because basis has three available derivative matrices
+function calc_metric_terms!(jacobian_matrix, element, basis::AbstractDerivativeOperator, node_coordinates)
+
+  # storage format:
+  #   jacobian_matrix[1,1,:,:,:] <- X_xi
+  #   jacobian_matrix[1,2,:,:,:] <- X_eta
+  #   jacobian_matrix[2,1,:,:,:] <- Y_xi
+  #   jacobian_matrix[2,2,:,:,:] <- Y_eta
+
+  # Compute the xi derivatives by applying D on the left
+  jacobian_matrix[1, 1, :, :, element] = Matrix(basis) * node_coordinates[1, :, :, element]
+  jacobian_matrix[2, 1, :, :, element] = Matrix(basis) * node_coordinates[2, :, :, element]
+
+  # Compute the eta derivatives by applying tranpose of D on the right
+  jacobian_matrix[1, 2, :, :, element] = node_coordinates[1, :, :, element] * Matrix(basis)'
+  jacobian_matrix[2, 2, :, :, element] = node_coordinates[2, :, :, element] * Matrix(basis)'
+
+  return jacobian_matrix
+end
+
+# construct the normal direction vectors (but not actually normalized) for a curved sided FDSBP element "block"
+# normalization occurs on the fly during the surface flux computation
+# OBS! This assumes that the boundary points are included.
+function calc_normal_directions!(normal_directions, element, jacobian_matrix)
+
+  # normal directions on the boundary for the left (local side 4) and right (local side 2)
+  N_ = size(jacobian_matrix, 3)
+  for j in 1:N_
+    # +x side or side 2 in the local indexing
+    X_xi  = jacobian_matrix[1, 1, N_, j, element]
+    X_eta = jacobian_matrix[1, 2, N_, j, element]
+    Y_xi  = jacobian_matrix[2, 1, N_, j, element]
+    Y_eta = jacobian_matrix[2, 2, N_, j, element]
+    Jtemp = X_xi * Y_eta - X_eta * Y_xi
+    normal_directions[1, j, 2, element] = sign(Jtemp) * ( Y_eta)
+    normal_directions[2, j, 2, element] = sign(Jtemp) * (-X_eta)
+
+    # -x side or side 4 in the local indexing
+    X_xi  = jacobian_matrix[1, 1, 1, j, element]
+    X_eta = jacobian_matrix[1, 2, 1, j, element]
+    Y_xi  = jacobian_matrix[2, 1, 1, j, element]
+    Y_eta = jacobian_matrix[2, 2, 1, j, element]
+    Jtemp = X_xi * Y_eta - X_eta * Y_xi
+    normal_directions[1, j, 4, element] = -sign(Jtemp) * ( Y_eta)
+    normal_directions[2, j, 4, element] = -sign(Jtemp) * (-X_eta)
+  end
+
+  # normal directions on the boundary for the top (local side 3) and bottom (local side 1)
+  N_ = size(jacobian_matrix, 4)
+  for i in 1:N_
+    # -y side or side 1 in the local indexing
+    X_xi  = jacobian_matrix[1, 1, i, 1, element]
+    X_eta = jacobian_matrix[1, 2, i, 1, element]
+    Y_xi  = jacobian_matrix[2, 1, i, 1, element]
+    Y_eta = jacobian_matrix[2, 2, i, 1, element]
+    Jtemp = X_xi * Y_eta - X_eta * Y_xi
+    normal_directions[1, i, 1, element] = -sign(Jtemp) * (-Y_xi)
+    normal_directions[2, i, 1, element] = -sign(Jtemp) * ( X_xi)
+
+    # +y side or side 3 in the local indexing
+    X_xi  = jacobian_matrix[1, 1, i, N_, element]
+    X_eta = jacobian_matrix[1, 2, i, N_, element]
+    Y_xi  = jacobian_matrix[2, 1, i, N_, element]
+    Y_eta = jacobian_matrix[2, 2, i, N_, element]
+    Jtemp = X_xi * Y_eta - X_eta * Y_xi
+    normal_directions[1, i, 3, element] = sign(Jtemp) * (-Y_xi)
+    normal_directions[2, i, 3, element] = sign(Jtemp) * ( X_xi)
+  end
+
+  return normal_directions
+end
+
+
+end # @muladd
\ No newline at end of file
diff --git a/src/solvers/fdsbp_unstructured/fdsbp.jl b/src/solvers/fdsbp_unstructured/fdsbp.jl
new file mode 100644
index 00000000000..8645d08de82
--- /dev/null
+++ b/src/solvers/fdsbp_unstructured/fdsbp.jl
@@ -0,0 +1,16 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
+
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+
+
+# dimension specific curvilinear implementations and data structures
+include("containers_2d.jl")
+include("fdsbp_2d.jl")
+
+
+end # @muladd
\ No newline at end of file
diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
new file mode 100644
index 00000000000..fb0d14b1568
--- /dev/null
+++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
@@ -0,0 +1,199 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
+
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+
+
+# 2D unstructured cache
+function create_cache(mesh::UnstructuredMesh2D, equations, dg::FDSBP, RealT, uEltype)
+
+  elements = init_elements(mesh, equations, dg.basis, RealT, uEltype)
+
+  interfaces = init_interfaces(mesh, elements)
+
+  boundaries = init_boundaries(mesh, elements)
+
+  cache = (; elements, interfaces, boundaries)
+
+  # Add specialized parts of the cache required to for efficient flux computations
+  cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...)
+
+  return cache
+end
+
+
+# TODO: FD; Upwind versions of surface / volume integral
+
+# 2D volume integral contributions for `VolumeIntegralStrongForm`
+@inline function calc_volume_integral!(du, u,
+                                       mesh::UnstructuredMesh2D,
+                                       nonconservative_terms::Val{false}, equations,
+                                       volume_integral::VolumeIntegralStrongForm,
+                                       dg::FDSBP, cache)
+  # Pull the derivative matrix
+  # TODO: FD, improve performance to use `mul!`. Current version is slow and allocates
+  D = Matrix(dg.basis) # SBP derivative operator
+  @unpack contravariant_vectors = cache.elements
+
+  @threaded for element in eachelement(dg, cache)
+    for j in eachnode(dg), i in eachnode(dg)
+      u_node = get_node_vars(u, equations, dg, i, j, element)
+
+      flux1 = flux(u_node, 1, equations)
+      flux2 = flux(u_node, 2, equations)
+
+      # Compute the contravariant flux by taking the scalar product of the
+      # first contravariant vector Ja^1 and the flux vector
+      Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
+      contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2
+      for ii in eachnode(dg)
+        multiply_add_to_node_vars!(du, D[ii, i], contravariant_flux1, equations, dg, ii, j, element)
+      end
+
+      # Compute the contravariant flux by taking the scalar product of the
+      # second contravariant vector Ja^2 and the flux vector
+      Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
+      contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2
+      for jj in eachnode(dg)
+        multiply_add_to_node_vars!(du, D[jj, j], contravariant_flux2, equations, dg, i, jj, element)
+      end
+    end
+  end
+
+  return nothing
+end
+
+
+# Note! The local side numbering for the unstructured quadrilateral element implementation differs
+#       from the structured TreeMesh or StructuredMesh local side numbering:
+#
+#      TreeMesh/StructuredMesh sides   versus   UnstructuredMesh sides
+#                  4                                  3
+#          -----------------                  -----------------
+#          |               |                  |               |
+#          | ^ eta         |                  | ^ eta         |
+#        1 | |             | 2              4 | |             | 2
+#          | |             |                  | |             |
+#          | ---> xi       |                  | ---> xi       |
+#          -----------------                  -----------------
+#                  3                                  1
+# Therefore, we require a different surface integral routine here despite their similar structure.
+# Also, the normal directions are already outward pointing for `UnstructuredMesh2D` so all the
+# surface contributions are added.
+function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D,
+                                equations, surface_integral::SurfaceIntegralStrongForm,
+                                dg::DG, cache)
+  inv_weight_left  = inv(left_boundary_weight(dg.basis))
+  inv_weight_right = inv(right_boundary_weight(dg.basis))
+  @unpack normal_directions, surface_flux_values = cache.elements
+
+  @threaded for element in eachelement(dg, cache)
+    for l in eachnode(dg)
+      # surface at -x
+      u_node = get_node_vars(u, equations, dg, 1, l, element)
+      # compute internal flux in normal direction on side 4
+      outward_direction = normal_directions[:, l, 4, element]
+      f_node = flux(u_node, outward_direction, equations)
+      f_num  = get_node_vars(surface_flux_values, equations, dg, l, 4, element)
+      multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node),
+                                 equations, dg, 1, l, element)
+
+      # surface at +x
+      u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element)
+      # compute internal flux in normal direction on side 2
+      outward_direction = normal_directions[:, l, 2, element]
+      f_node = flux(u_node, outward_direction, equations)
+      f_num  = get_node_vars(surface_flux_values, equations, dg, l, 2, element)
+      multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node),
+                                 equations, dg, nnodes(dg), l, element)
+
+      # surface at -y
+      u_node = get_node_vars(u, equations, dg, l, 1, element)
+      # compute internal flux in normal direction on side 1
+      outward_direction = normal_directions[:, l, 1, element]
+      f_node = flux(u_node, outward_direction, equations)
+      f_num  = get_node_vars(surface_flux_values, equations, dg, l, 1, element)
+      multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node),
+                                 equations, dg, l, 1, element)
+
+      # surface at +y
+      u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element)
+      # compute internal flux in normal direction on side 3
+      outward_direction = normal_directions[:, l, 3, element]
+      f_node = flux(u_node, outward_direction, equations)
+      f_num  = get_node_vars(surface_flux_values, equations, dg, l, 3, element)
+      multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node),
+                                 equations, dg, l, nnodes(dg), element)
+    end
+  end
+
+  return nothing
+end
+
+
+# AnalysisCallback
+function integrate_via_indices(func::Func, u,
+                               mesh::UnstructuredMesh2D, equations,
+                               dg::FDSBP, cache, args...; normalize=true) where {Func}
+  # TODO: FD. This is rather inefficient right now and allocates...
+  weights = diag(SummationByPartsOperators.mass_matrix(dg.basis))
+
+  # Initialize integral with zeros of the right shape
+  integral = zero(func(u, 1, 1, 1, equations, dg, args...))
+  total_volume = zero(real(mesh))
+
+  # Use quadrature to numerically integrate over entire domain
+  for element in eachelement(dg, cache)
+    for j in eachnode(dg), i in eachnode(dg)
+      volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, element]))
+      integral += volume_jacobian * weights[i] * weights[j] * func(u, i, j, element, equations, dg, args...)
+      total_volume += volume_jacobian * weights[i] * weights[j]
+    end
+  end
+
+  # Normalize with total volume
+  if normalize
+    integral = integral / total_volume
+  end
+
+  return integral
+end
+
+
+function calc_error_norms(func, u, t, analyzer,
+                          mesh::UnstructuredMesh2D, equations, initial_condition,
+                          dg::FDSBP, cache, cache_analysis)
+  # TODO: FD. This is rather inefficient right now and allocates...
+  weights = diag(SummationByPartsOperators.mass_matrix(dg.basis))
+  @unpack node_coordinates, inverse_jacobian = cache.elements
+
+  # Set up data structures
+  l2_error   = zero(func(get_node_vars(u, equations, dg, 1, 1, 1), equations))
+  linf_error = copy(l2_error)
+  total_volume = zero(real(mesh))
+
+  # Iterate over all elements for error calculations
+  for element in eachelement(dg, cache)
+    for j in eachnode(analyzer), i in eachnode(analyzer)
+      volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, element]))
+      u_exact = initial_condition(
+        get_node_coords(node_coordinates, equations, dg, i, j, element), t, equations)
+      diff = func(u_exact, equations) - func(
+        get_node_vars(u, equations, dg, i, j, element), equations)
+      l2_error += diff.^2 * (weights[i] * weights[j] * volume_jacobian)
+      linf_error = @. max(linf_error, abs(diff))
+      total_volume += weights[i] * weights[j] * volume_jacobian
+    end
+  end
+
+  # For L2 error, divide by total volume
+  l2_error = @. sqrt(l2_error / total_volume)
+
+  return l2_error, linf_error
+end
+
+end # @muladd

From 00b0effd2b039b221d8930b5b026723db2f15158 Mon Sep 17 00:00:00 2001
From: Andrew Winters <andrew.ross.winters@liu.se>
Date: Mon, 12 Dec 2022 19:51:33 +0100
Subject: [PATCH 02/15] new elixirs for free-stream and convergence test

---
 .../elixir_euler_convergence.jl               | 64 +++++++++++++++
 .../elixir_euler_free_stream.jl               | 78 +++++++++++++++++++
 2 files changed, 142 insertions(+)
 create mode 100644 examples/unstructured_2d_fdsbp/elixir_euler_convergence.jl
 create mode 100644 examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl

diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_convergence.jl b/examples/unstructured_2d_fdsbp/elixir_euler_convergence.jl
new file mode 100644
index 00000000000..23c07b3300b
--- /dev/null
+++ b/examples/unstructured_2d_fdsbp/elixir_euler_convergence.jl
@@ -0,0 +1,64 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the FDSBP approximation operator
+
+D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
+                            derivative_order=1, accuracy_order=4,
+                            xmin=-1.0, xmax=1.0, N=12)
+solver = FDSBP(D_SBP,
+               surface_integral=SurfaceIntegralStrongForm(flux_lax_friedrichs),
+               volume_integral=VolumeIntegralStrongForm())
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+
+default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")
+isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh",
+                                       default_mesh_file)
+mesh_file = default_mesh_file
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity=true)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms=source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(), abstol=1.0e-9, reltol=1.0e-9,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary
diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl
new file mode 100644
index 00000000000..e2846e24eb5
--- /dev/null
+++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl
@@ -0,0 +1,78 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Free-stream initial condition
+initial_condition = initial_condition_constant
+
+# Boundary conditions for free-stream testing
+boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict( :Body    => boundary_condition_free_stream,
+                            :Button1 => boundary_condition_free_stream,
+                            :Button2 => boundary_condition_free_stream,
+                            :Eye1    => boundary_condition_free_stream,
+                            :Eye2    => boundary_condition_free_stream,
+                            :Smile   => boundary_condition_free_stream,
+                            :Bowtie  => boundary_condition_free_stream )
+
+###############################################################################
+# Get the FDSBP approximation space
+
+D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
+                            derivative_order=1, accuracy_order=4,
+                            xmin=-1.0, xmax=1.0, N=9)
+solver = FDSBP(D_SBP,
+               surface_integral=SurfaceIntegralStrongForm(flux_hll),
+               volume_integral=VolumeIntegralStrongForm())
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+
+default_mesh_file = joinpath(@__DIR__, "mesh_gingerbread_man.mesh")
+isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/2c6440b5f8a57db131061ad7aa78ee2b/raw/1f89fdf2c874ff678c78afb6fe8dc784bdfd421f/mesh_gingerbread_man.mesh",
+                                       default_mesh_file)
+mesh_file = default_mesh_file
+
+mesh = UnstructuredMesh2D(mesh_file)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# set small tolerances for the free-stream preservation test
+sol = solve(ode, SSPRK43(), abstol=1.0e-12, reltol=1.0e-12,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary

From 6229c8572766edcc16c09058ba3e1d7b44b57249 Mon Sep 17 00:00:00 2001
From: Andrew Winters <andrew.ross.winters@liu.se>
Date: Mon, 12 Dec 2022 19:51:57 +0100
Subject: [PATCH 03/15] add tests

---
 test/test_unstructured_2d.jl | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl
index 882309d70da..6cfe2de6115 100644
--- a/test/test_unstructured_2d.jl
+++ b/test/test_unstructured_2d.jl
@@ -155,6 +155,21 @@ isdir(outdir) && rm(outdir, recursive=true)
       linf = [2.7639232436274392, 3.3985508653311767, 3.3330308209196224, 2.052861364219655],
       tspan = (0.0, 0.25))
   end
+
+  # TODO: FD; for now put the unstructured tests for the 2D FDSBP here.
+  @trixi_testset "FDSBP (central): elixir_euler_free_stream.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "unstructured_2d_fdsbp", "elixir_euler_free_stream.jl"),
+      l2   = [1.0782928209546407e-14, 7.961856579424368e-14, 5.54300886580838e-14, 1.160015193708739e-13],
+      linf = [3.410882687404637e-11, 1.5696624555694427e-11, 9.794248745365053e-12, 2.504663143554353e-11],
+      tspan=(0.0, 0.2))
+  end
+
+  @trixi_testset "FDSBP (central): elixir_euler_convergence.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "unstructured_2d_fdsbp", "elixir_euler_convergence.jl"),
+      l2   = [3.372206837157435e-5, 6.001247046050437e-5, 6.0012470460493534e-5, 0.00018518205373473714],
+      linf = [0.00015638063264633573, 0.0003329723828695563, 0.00033297238286467135, 0.0009179128246792345],
+      tspan=(0.0, 0.025))
+  end
 end
 
 # Clean up afterwards: delete Trixi output directory

From ef47e92e07addd085205a53b10e7aa7e0fc490e9 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Tue, 13 Dec 2022 09:49:58 +0100
Subject: [PATCH 04/15] avoid allocations in init_element!

---
 .../fdsbp_unstructured/containers_2d.jl       | 29 +++++++++++++++----
 1 file changed, 23 insertions(+), 6 deletions(-)

diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl
index 4e0945dbd6d..c18a5b24745 100644
--- a/src/solvers/fdsbp_unstructured/containers_2d.jl
+++ b/src/solvers/fdsbp_unstructured/containers_2d.jl
@@ -47,9 +47,8 @@ end
 
 # construct the metric terms for a FDSBP element "block". Directly use the derivative matrix
 # applied to the node coordinates.
-# TODO: FD performance; this is slow and allocates. Rewrite to use `mul!`
 # TODO: FD; How to make this work for the upwind solver because basis has three available derivative matrices
-function calc_metric_terms!(jacobian_matrix, element, basis::AbstractDerivativeOperator, node_coordinates)
+function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeOperator, node_coordinates)
 
   # storage format:
   #   jacobian_matrix[1,1,:,:,:] <- X_xi
@@ -58,12 +57,30 @@ function calc_metric_terms!(jacobian_matrix, element, basis::AbstractDerivativeO
   #   jacobian_matrix[2,2,:,:,:] <- Y_eta
 
   # Compute the xi derivatives by applying D on the left
-  jacobian_matrix[1, 1, :, :, element] = Matrix(basis) * node_coordinates[1, :, :, element]
-  jacobian_matrix[2, 1, :, :, element] = Matrix(basis) * node_coordinates[2, :, :, element]
+  # This is basically the same as
+  # jacobian_matrix[1, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[1, :, :, element]
+  # but uses only matrix-vector products instead of a matrix-matrix product.
+  for j in eachnode(D_SBP)
+    mul!(view(jacobian_matrix, 1, 1, :, j, element), D_SBP,
+         view(node_coordinates, 1, :, j, element))
+  end
+  # jacobian_matrix[2, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[2, :, :, element]
+  for j in eachnode(D_SBP)
+    mul!(view(jacobian_matrix, 2, 1, :, j, element), D_SBP,
+         view(node_coordinates, 2, :, j, element))
+  end
 
   # Compute the eta derivatives by applying tranpose of D on the right
-  jacobian_matrix[1, 2, :, :, element] = node_coordinates[1, :, :, element] * Matrix(basis)'
-  jacobian_matrix[2, 2, :, :, element] = node_coordinates[2, :, :, element] * Matrix(basis)'
+  # jacobian_matrix[1, 2, :, :, element] = node_coordinates[1, :, :, element] * Matrix(D_SBP)'
+  for i in eachnode(D_SBP)
+    mul!(view(jacobian_matrix, 1, 2, i, :, element), D_SBP,
+         view(node_coordinates, 1, i, :, element))
+  end
+  # jacobian_matrix[2, 2, :, :, element] = node_coordinates[2, :, :, element] * Matrix(D_SBP)'
+  for i in eachnode(D_SBP)
+    mul!(view(jacobian_matrix, 2, 2, i, :, element), D_SBP,
+         view(node_coordinates, 2, i, :, element))
+  end
 
   return jacobian_matrix
 end

From f512d0358f0da8de69a2733fe0475234d45c6b2c Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Tue, 13 Dec 2022 10:04:30 +0100
Subject: [PATCH 05/15] avoid allocations in calc_surface_integral!

---
 src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
index fb0d14b1568..60ac7cdb35e 100644
--- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
+++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
@@ -96,7 +96,7 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D,
       # surface at -x
       u_node = get_node_vars(u, equations, dg, 1, l, element)
       # compute internal flux in normal direction on side 4
-      outward_direction = normal_directions[:, l, 4, element]
+      outward_direction = get_node_coords(normal_directions, equations, dg, l, 4, element)
       f_node = flux(u_node, outward_direction, equations)
       f_num  = get_node_vars(surface_flux_values, equations, dg, l, 4, element)
       multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node),
@@ -105,7 +105,7 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D,
       # surface at +x
       u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element)
       # compute internal flux in normal direction on side 2
-      outward_direction = normal_directions[:, l, 2, element]
+      outward_direction = get_node_coords(normal_directions, equations, dg, l, 2, element)
       f_node = flux(u_node, outward_direction, equations)
       f_num  = get_node_vars(surface_flux_values, equations, dg, l, 2, element)
       multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node),
@@ -114,7 +114,7 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D,
       # surface at -y
       u_node = get_node_vars(u, equations, dg, l, 1, element)
       # compute internal flux in normal direction on side 1
-      outward_direction = normal_directions[:, l, 1, element]
+      outward_direction = get_node_coords(normal_directions, equations, dg, l, 1, element)
       f_node = flux(u_node, outward_direction, equations)
       f_num  = get_node_vars(surface_flux_values, equations, dg, l, 1, element)
       multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node),
@@ -123,7 +123,7 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D,
       # surface at +y
       u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element)
       # compute internal flux in normal direction on side 3
-      outward_direction = normal_directions[:, l, 3, element]
+      outward_direction = get_node_coords(normal_directions, equations, dg, l, 3, element)
       f_node = flux(u_node, outward_direction, equations)
       f_num  = get_node_vars(surface_flux_values, equations, dg, l, 3, element)
       multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node),

From 5f968f75ee89e9fd99d9b64bc201a2f7b005638f Mon Sep 17 00:00:00 2001
From: Andrew Winters <andrew.ross.winters@liu.se>
Date: Tue, 13 Dec 2022 12:01:42 +0100
Subject: [PATCH 06/15] Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl | 2 --
 src/solvers/fdsbp_unstructured/containers_2d.jl            | 2 +-
 2 files changed, 1 insertion(+), 3 deletions(-)

diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl
index e2846e24eb5..e9aa345e80a 100644
--- a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl
+++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl
@@ -64,8 +64,6 @@ save_solution = SaveSolutionCallback(interval=100,
                                      save_initial_solution=true,
                                      save_final_solution=true)
 
-stepsize_callback = StepsizeCallback(cfl=1.0)
-
 callbacks = CallbackSet(summary_callback, analysis_callback,
                         alive_callback, save_solution)
 
diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl
index c18a5b24745..6426cda293f 100644
--- a/src/solvers/fdsbp_unstructured/containers_2d.jl
+++ b/src/solvers/fdsbp_unstructured/containers_2d.jl
@@ -70,7 +70,7 @@ function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeO
          view(node_coordinates, 2, :, j, element))
   end
 
-  # Compute the eta derivatives by applying tranpose of D on the right
+  # Compute the eta derivatives by applying transpose of D on the right
   # jacobian_matrix[1, 2, :, :, element] = node_coordinates[1, :, :, element] * Matrix(D_SBP)'
   for i in eachnode(D_SBP)
     mul!(view(jacobian_matrix, 1, 2, i, :, element), D_SBP,

From 3cf1886cebaa1301770de44d2281ba658c90449f Mon Sep 17 00:00:00 2001
From: Andrew Winters <andrew.ross.winters@liu.se>
Date: Tue, 13 Dec 2022 14:22:56 +0100
Subject: [PATCH 07/15] add helper function to get the derivative matrix from
 basis

---
 src/solvers/fdsbp_tree/fdsbp.jl                 | 2 ++
 src/solvers/fdsbp_unstructured/containers_2d.jl | 3 ++-
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/solvers/fdsbp_tree/fdsbp.jl b/src/solvers/fdsbp_tree/fdsbp.jl
index b89d59c1156..c31e2c12748 100644
--- a/src/solvers/fdsbp_tree/fdsbp.jl
+++ b/src/solvers/fdsbp_tree/fdsbp.jl
@@ -36,6 +36,8 @@ end
 nnodes(D::AbstractDerivativeOperator) = size(D, 1)
 eachnode(D::AbstractDerivativeOperator) = Base.OneTo(nnodes(D))
 get_nodes(D::AbstractDerivativeOperator) = grid(D)
+get_derivative_matrix(D::DerivativeOperator) = D
+get_derivative_matrix(D::SummationByPartsOperators.UpwindOperators) = D.central
 
 # TODO: This is hack to enable the FDSBP solver to use the
 #       `SaveSolutionCallback`.
diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl
index 6426cda293f..ba0ef86b099 100644
--- a/src/solvers/fdsbp_unstructured/containers_2d.jl
+++ b/src/solvers/fdsbp_unstructured/containers_2d.jl
@@ -48,7 +48,7 @@ end
 # construct the metric terms for a FDSBP element "block". Directly use the derivative matrix
 # applied to the node coordinates.
 # TODO: FD; How to make this work for the upwind solver because basis has three available derivative matrices
-function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeOperator, node_coordinates)
+function calc_metric_terms!(jacobian_matrix, element, basis::AbstractDerivativeOperator, node_coordinates)
 
   # storage format:
   #   jacobian_matrix[1,1,:,:,:] <- X_xi
@@ -56,6 +56,7 @@ function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeO
   #   jacobian_matrix[2,1,:,:,:] <- Y_xi
   #   jacobian_matrix[2,2,:,:,:] <- Y_eta
 
+  D_SBP = get_derivative_matrix(basis)
   # Compute the xi derivatives by applying D on the left
   # This is basically the same as
   # jacobian_matrix[1, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[1, :, :, element]

From 14b1d96bdc4669251ff29d3d21e7dd7d3c76becb Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Tue, 13 Dec 2022 19:37:02 +0100
Subject: [PATCH 08/15] increased tolerance for free stream test

---
 test/test_unstructured_2d.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl
index 6cfe2de6115..ab15032a02e 100644
--- a/test/test_unstructured_2d.jl
+++ b/test/test_unstructured_2d.jl
@@ -161,14 +161,14 @@ isdir(outdir) && rm(outdir, recursive=true)
     @test_trixi_include(joinpath(examples_dir(), "unstructured_2d_fdsbp", "elixir_euler_free_stream.jl"),
       l2   = [1.0782928209546407e-14, 7.961856579424368e-14, 5.54300886580838e-14, 1.160015193708739e-13],
       linf = [3.410882687404637e-11, 1.5696624555694427e-11, 9.794248745365053e-12, 2.504663143554353e-11],
-      tspan=(0.0, 0.2))
+      tspan = (0.0, 0.2), atol = 1.0e-11)
   end
 
   @trixi_testset "FDSBP (central): elixir_euler_convergence.jl" begin
     @test_trixi_include(joinpath(examples_dir(), "unstructured_2d_fdsbp", "elixir_euler_convergence.jl"),
       l2   = [3.372206837157435e-5, 6.001247046050437e-5, 6.0012470460493534e-5, 0.00018518205373473714],
       linf = [0.00015638063264633573, 0.0003329723828695563, 0.00033297238286467135, 0.0009179128246792345],
-      tspan=(0.0, 0.025))
+      tspan = (0.0, 0.025))
   end
 end
 

From 46a986649fea9501227ad8f3c27b0d2b4f08f60d Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Tue, 13 Dec 2022 19:58:41 +0100
Subject: [PATCH 09/15] WIP: improve performance of volume terms

---
 src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 41 ++++++++++++++++++++++
 1 file changed, 41 insertions(+)

diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
index 60ac7cdb35e..18bc114b580 100644
--- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
+++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
@@ -34,6 +34,47 @@ end
                                        nonconservative_terms::Val{false}, equations,
                                        volume_integral::VolumeIntegralStrongForm,
                                        dg::FDSBP, cache)
+  # D = dg.basis # SBP derivative operator
+  # @unpack f_threaded = cache
+
+  # # SBP operators from SummationByPartsOperators.jl implement the basic interface
+  # # of matrix-vector multiplication. Thus, we pass an "array of structures",
+  # # packing all variables per node in an `SVector`.
+  # if nvariables(equations) == 1
+  #   # `reinterpret(reshape, ...)` removes the leading dimension only if more
+  #   # than one variable is used.
+  #   u_vectors  = reshape(reinterpret(SVector{nvariables(equations), eltype(u)}, u),
+  #                        nnodes(dg), nnodes(dg), nelements(dg, cache))
+  #   du_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(du)}, du),
+  #                        nnodes(dg), nnodes(dg), nelements(dg, cache))
+  # else
+  #   u_vectors  = reinterpret(reshape, SVector{nvariables(equations), eltype(u)}, u)
+  #   du_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(du)}, du)
+  # end
+
+  # # Use the tensor product structure to compute the discrete derivatives of
+  # # the contravariant fluxes line-by-line and add them to `du` for each element.
+  # @threaded for element in eachelement(dg, cache)
+  #   f_element = f_threaded[Threads.threadid()]
+  #   u_element = view(u_vectors,  :, :, element)
+
+  #   # x direction
+  #   @. f_element = flux(u_element, 1, equations)
+  #   for j in eachnode(dg)
+  #     mul!(view(du_vectors, :, j, element), D, view(f_element, :, j),
+  #          one(eltype(du)), one(eltype(du)))
+  #   end
+
+  #   # y direction
+  #   @. f_element = flux(u_element, 2, equations)
+  #   for i in eachnode(dg)
+  #     mul!(view(du_vectors, i, :, element), D, view(f_element, i, :),
+  #          one(eltype(du)), one(eltype(du)))
+  #   end
+  # end
+
+
+
   # Pull the derivative matrix
   # TODO: FD, improve performance to use `mul!`. Current version is slow and allocates
   D = Matrix(dg.basis) # SBP derivative operator

From e3b27721bf1ca148abec1f60da8ca86194a3630f Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Tue, 13 Dec 2022 20:25:29 +0100
Subject: [PATCH 10/15] remove allocations in volume terms

---
 src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 103 ++++++++-------------
 1 file changed, 39 insertions(+), 64 deletions(-)

diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
index 18bc114b580..73c6e1524a1 100644
--- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
+++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
@@ -34,74 +34,49 @@ end
                                        nonconservative_terms::Val{false}, equations,
                                        volume_integral::VolumeIntegralStrongForm,
                                        dg::FDSBP, cache)
-  # D = dg.basis # SBP derivative operator
-  # @unpack f_threaded = cache
-
-  # # SBP operators from SummationByPartsOperators.jl implement the basic interface
-  # # of matrix-vector multiplication. Thus, we pass an "array of structures",
-  # # packing all variables per node in an `SVector`.
-  # if nvariables(equations) == 1
-  #   # `reinterpret(reshape, ...)` removes the leading dimension only if more
-  #   # than one variable is used.
-  #   u_vectors  = reshape(reinterpret(SVector{nvariables(equations), eltype(u)}, u),
-  #                        nnodes(dg), nnodes(dg), nelements(dg, cache))
-  #   du_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(du)}, du),
-  #                        nnodes(dg), nnodes(dg), nelements(dg, cache))
-  # else
-  #   u_vectors  = reinterpret(reshape, SVector{nvariables(equations), eltype(u)}, u)
-  #   du_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(du)}, du)
-  # end
-
-  # # Use the tensor product structure to compute the discrete derivatives of
-  # # the contravariant fluxes line-by-line and add them to `du` for each element.
-  # @threaded for element in eachelement(dg, cache)
-  #   f_element = f_threaded[Threads.threadid()]
-  #   u_element = view(u_vectors,  :, :, element)
-
-  #   # x direction
-  #   @. f_element = flux(u_element, 1, equations)
-  #   for j in eachnode(dg)
-  #     mul!(view(du_vectors, :, j, element), D, view(f_element, :, j),
-  #          one(eltype(du)), one(eltype(du)))
-  #   end
-
-  #   # y direction
-  #   @. f_element = flux(u_element, 2, equations)
-  #   for i in eachnode(dg)
-  #     mul!(view(du_vectors, i, :, element), D, view(f_element, i, :),
-  #          one(eltype(du)), one(eltype(du)))
-  #   end
-  # end
-
-
-
-  # Pull the derivative matrix
-  # TODO: FD, improve performance to use `mul!`. Current version is slow and allocates
-  D = Matrix(dg.basis) # SBP derivative operator
+  D = dg.basis # SBP derivative operator
+  @unpack f_threaded = cache
   @unpack contravariant_vectors = cache.elements
 
+  # SBP operators from SummationByPartsOperators.jl implement the basic interface
+  # of matrix-vector multiplication. Thus, we pass an "array of structures",
+  # packing all variables per node in an `SVector`.
+  if nvariables(equations) == 1
+    # `reinterpret(reshape, ...)` removes the leading dimension only if more
+    # than one variable is used.
+    u_vectors  = reshape(reinterpret(SVector{nvariables(equations), eltype(u)}, u),
+                         nnodes(dg), nnodes(dg), nelements(dg, cache))
+    du_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(du)}, du),
+                         nnodes(dg), nnodes(dg), nelements(dg, cache))
+  else
+    u_vectors  = reinterpret(reshape, SVector{nvariables(equations), eltype(u)}, u)
+    du_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(du)}, du)
+  end
+
+  # Use the tensor product structure to compute the discrete derivatives of
+  # the contravariant fluxes line-by-line and add them to `du` for each element.
   @threaded for element in eachelement(dg, cache)
+    f_element = f_threaded[Threads.threadid()]
+    u_element = view(u_vectors,  :, :, element)
+
+    # x direction
+    for j in eachnode(dg), i in eachnode(dg)
+      Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
+      f_element[i, j] = flux(u_element[i, j], Ja1, equations)
+    end
+    for j in eachnode(dg)
+      mul!(view(du_vectors, :, j, element), D, view(f_element, :, j),
+           one(eltype(du)), one(eltype(du)))
+    end
+
+    # y direction
     for j in eachnode(dg), i in eachnode(dg)
-      u_node = get_node_vars(u, equations, dg, i, j, element)
-
-      flux1 = flux(u_node, 1, equations)
-      flux2 = flux(u_node, 2, equations)
-
-      # Compute the contravariant flux by taking the scalar product of the
-      # first contravariant vector Ja^1 and the flux vector
-      Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
-      contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2
-      for ii in eachnode(dg)
-        multiply_add_to_node_vars!(du, D[ii, i], contravariant_flux1, equations, dg, ii, j, element)
-      end
-
-      # Compute the contravariant flux by taking the scalar product of the
-      # second contravariant vector Ja^2 and the flux vector
-      Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
-      contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2
-      for jj in eachnode(dg)
-        multiply_add_to_node_vars!(du, D[jj, j], contravariant_flux2, equations, dg, i, jj, element)
-      end
+      Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
+      f_element[i, j] = flux(u_element[i, j], Ja2, equations)
+    end
+    for i in eachnode(dg)
+      mul!(view(du_vectors, i, :, element), D, view(f_element, i, :),
+           one(eltype(du)), one(eltype(du)))
     end
   end
 

From 6e6aedc4bdf83c3b21312b9839da9f1b8ede7f11 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Tue, 13 Dec 2022 20:27:13 +0100
Subject: [PATCH 11/15] improve performance of volume terms

---
 src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
index 73c6e1524a1..2bc0fe66225 100644
--- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
+++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
@@ -60,21 +60,21 @@ end
     u_element = view(u_vectors,  :, :, element)
 
     # x direction
-    for j in eachnode(dg), i in eachnode(dg)
-      Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
-      f_element[i, j] = flux(u_element[i, j], Ja1, equations)
-    end
     for j in eachnode(dg)
+      for i in eachnode(dg)
+        Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
+        f_element[i, j] = flux(u_element[i, j], Ja1, equations)
+      end
       mul!(view(du_vectors, :, j, element), D, view(f_element, :, j),
            one(eltype(du)), one(eltype(du)))
     end
 
     # y direction
-    for j in eachnode(dg), i in eachnode(dg)
-      Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
-      f_element[i, j] = flux(u_element[i, j], Ja2, equations)
-    end
     for i in eachnode(dg)
+      for j in eachnode(dg)
+        Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
+        f_element[i, j] = flux(u_element[i, j], Ja2, equations)
+      end
       mul!(view(du_vectors, i, :, element), D, view(f_element, i, :),
            one(eltype(du)), one(eltype(du)))
     end

From c2d187722cff556d8852dc4e0bbbe1864e09d4c4 Mon Sep 17 00:00:00 2001
From: Andrew Winters <andrew.ross.winters@liu.se>
Date: Thu, 15 Dec 2022 12:37:59 +0100
Subject: [PATCH 12/15] fix summary output when discretization uses
 VolumeIntegralStrongForm

---
 src/solvers/dg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
index 487e7f1a476..a262889872b 100644
--- a/src/solvers/dg.jl
+++ b/src/solvers/dg.jl
@@ -349,7 +349,7 @@ function Base.show(io::IO, mime::MIME"text/plain", dg::DG)
     summary_line(io, "surface integral", dg.surface_integral |> typeof |> nameof)
     show(increment_indent(io), mime, dg.surface_integral)
     summary_line(io, "volume integral", dg.volume_integral |> typeof |> nameof)
-    if !(dg.volume_integral isa VolumeIntegralWeakForm)
+    if !(dg.volume_integral isa VolumeIntegralWeakForm) && !(dg.volume_integral isa VolumeIntegralStrongForm)
       show(increment_indent(io), mime, dg.volume_integral)
     end
     summary_footer(io)

From 101aa5757b61734d573ac442b25470a1c90b1280 Mon Sep 17 00:00:00 2001
From: Andrew Winters <andrew.ross.winters@liu.se>
Date: Thu, 29 Dec 2022 19:36:38 +0100
Subject: [PATCH 13/15] unify UnstructuredElement2D constructors for DGSEM and
 FDSBP. Removes duplicate code

---
 .../dgsem_unstructured/containers_2d.jl       | 14 ++++++-------
 .../fdsbp_unstructured/containers_2d.jl       | 20 -------------------
 2 files changed, 7 insertions(+), 27 deletions(-)

diff --git a/src/solvers/dgsem_unstructured/containers_2d.jl b/src/solvers/dgsem_unstructured/containers_2d.jl
index 3417adc30d2..3ac674ab9b8 100644
--- a/src/solvers/dgsem_unstructured/containers_2d.jl
+++ b/src/solvers/dgsem_unstructured/containers_2d.jl
@@ -62,36 +62,36 @@ function init_elements(mesh::UnstructuredMesh2D, equations, basis, RealT, uEltyp
 end
 
 
-function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis::LobattoLegendreBasis)
+function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis)
   four_corners = zeros(eltype(mesh.corners), 4, 2)
 
   # loop through elements and call the correct constructor based on whether the element is curved
   for element in eachelement(elements)
     if mesh.element_is_curved[element]
-      init_element!(elements, element, basis.nodes, view(mesh.surface_curves, :, element))
+      init_element!(elements, element, basis, view(mesh.surface_curves, :, element))
     else # straight sided element
       for i in 1:4, j in 1:2
         # pull the (x,y) values of these corners out of the global corners array
         four_corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]]
       end
-      init_element!(elements, element, basis.nodes, four_corners)
+      init_element!(elements, element, basis, four_corners)
     end
   end
 end
 
 
 # initialize all the values in the container of a general DG element (either straight sided or curved)
-function init_element!(elements, element, nodes, corners_or_surface_curves)
+function init_element!(elements, element, basis::LobattoLegendreBasis, corners_or_surface_curves)
 
-  calc_node_coordinates!(elements.node_coordinates, element, nodes, corners_or_surface_curves)
+  calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis), corners_or_surface_curves)
 
-  calc_metric_terms!(elements.jacobian_matrix, element, nodes, corners_or_surface_curves)
+  calc_metric_terms!(elements.jacobian_matrix, element, get_nodes(basis), corners_or_surface_curves)
 
   calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix)
 
   calc_contravariant_vectors!(elements.contravariant_vectors, element, elements.jacobian_matrix)
 
-  calc_normal_directions!(elements.normal_directions, element, nodes, corners_or_surface_curves)
+  calc_normal_directions!(elements.normal_directions, element, get_nodes(basis), corners_or_surface_curves)
 
   return elements
 end
diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl
index ba0ef86b099..20a056ad64d 100644
--- a/src/solvers/fdsbp_unstructured/containers_2d.jl
+++ b/src/solvers/fdsbp_unstructured/containers_2d.jl
@@ -7,26 +7,6 @@
 # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
 @muladd begin
 
-# Alternative constructor to reuse the element container from `UnstructuredMesh2D`
-# OBS! New element constructors require the derivative operator to compute the metric terms.
-function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis::AbstractDerivativeOperator)
-  four_corners = zeros(eltype(mesh.corners), 4, 2)
-
-  # loop through elements and call the correct constructor based on whether the element is curved
-  for element in eachelement(elements)
-    if mesh.element_is_curved[element]
-      init_element!(elements, element, basis, view(mesh.surface_curves, :, element))
-    else # straight sided element
-      for i in 1:4, j in 1:2
-        # pull the (x,y) values of these corners out of the global corners array
-        four_corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]]
-      end
-      init_element!(elements, element, basis, four_corners)
-    end
-  end
-end
-
-
 # initialize all the values in the container of a general FD block (either straight sided or curved)
 # OBS! Requires the SBP derivative matrix in order to compute metric terms that are free-stream preserving
 function init_element!(elements, element, basis::AbstractDerivativeOperator, corners_or_surface_curves)

From 3c2b1e64ec3ac41233ecbe4571811397d2869c0d Mon Sep 17 00:00:00 2001
From: Andrew Winters <andrew.ross.winters@liu.se>
Date: Mon, 16 Jan 2023 21:09:04 +0100
Subject: [PATCH 14/15] new container type for biased metric terms. Upwind
 volume integral with biased metric terms

---
 src/equations/compressible_euler_2d.jl        | 42 ++++++++-
 src/equations/numerical_fluxes.jl             | 14 ++-
 .../dgsem_unstructured/containers_2d.jl       | 92 ++++++++++++++++++-
 src/solvers/fdsbp_tree/fdsbp.jl               |  2 -
 .../fdsbp_unstructured/containers_2d.jl       | 71 +++++++++++++-
 src/solvers/fdsbp_unstructured/fdsbp_2d.jl    | 85 ++++++++++++++++-
 6 files changed, 292 insertions(+), 14 deletions(-)

diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
index 65d5310d5b1..be9fa4eb0b2 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -780,7 +780,7 @@ end
     rho_2gamma = 0.5 * rho / equations.gamma
     f1m = rho_2gamma * alpha_m
     f2m = rho_2gamma * alpha_m * v1
-    f3m = rho_2gamma * (alpha_m * v2 + a * (lambda2_m-lambda3_m))
+    f3m = rho_2gamma * (alpha_m * v2 + a * (lambda2_m - lambda3_m))
     f4m = rho_2gamma * (alpha_m * 0.5 * (v1^2 + v2^2) + a * v2 * (lambda2_m - lambda3_m)
                         + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one)
   end
@@ -970,6 +970,46 @@ end
   return SVector(f1m, f2m, f3m, f4m)
 end
 
+@inline function splitting_lax_friedrichs(u, ::Val{:plus}, normal_direction::AbstractVector,
+                                          equations::CompressibleEulerEquations2D)
+  rho, rho_v1, rho_v2, rho_e = u
+  v1 = rho_v1 / rho
+  v2 = rho_v2 / rho
+  v_dot_n = normal_direction[1] * v1 + normal_direction[2] * v2
+  p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2))
+
+  a = sqrt(equations.gamma * p / rho)
+  H = (rho_e + p) / rho
+  lambda = 0.5 * (v_dot_n + a * norm(normal_direction))
+
+  f1p = 0.5 * rho * v_dot_n + lambda * u[1]
+  f2p = 0.5 * rho * v_dot_n * v1 + 0.5 * p * normal_direction[1] + lambda * u[2]
+  f3p = 0.5 * rho * v_dot_n * v2 + 0.5 * p * normal_direction[2] + lambda * u[3]
+  f4p = 0.5 * rho * v_dot_n * H + lambda * u[4]
+
+  return SVector(f1p, f2p, f3p, f4p)
+end
+
+@inline function splitting_lax_friedrichs(u, ::Val{:minus}, normal_direction::AbstractVector,
+                                          equations::CompressibleEulerEquations2D)
+  rho, rho_v1, rho_v2, rho_e = u
+  v1 = rho_v1 / rho
+  v2 = rho_v2 / rho
+  v_dot_n = normal_direction[1] * v1 + normal_direction[2] * v2
+  p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2))
+
+  a = sqrt(equations.gamma * p / rho)
+  H = (rho_e + p) / rho
+  lambda = 0.5 * (v_dot_n + a * norm(normal_direction))
+
+  f1m = 0.5 * rho * v_dot_n - lambda * u[1]
+  f2m = 0.5 * rho * v_dot_n * v1 + 0.5 * p * normal_direction[1] - lambda * u[2]
+  f3m = 0.5 * rho * v_dot_n * v2 + 0.5 * p * normal_direction[2] - lambda * u[3]
+  f4m = 0.5 * rho * v_dot_n * H - lambda * u[4]
+
+  return SVector(f1m, f2m, f3m, f4m)
+end
+
 
 # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the
 # maximum velocity magnitude plus the maximum speed of sound
diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl
index ff9596848bb..31722bfc4d9 100644
--- a/src/equations/numerical_fluxes.jl
+++ b/src/equations/numerical_fluxes.jl
@@ -333,13 +333,25 @@ struct FluxUpwind{Splitting}
   splitting::Splitting
 end
 
-@inline function (numflux::FluxUpwind)(u_ll, u_rr, orientation::Int, equations)
+@inline function (numflux::FluxUpwind)(u_ll, u_rr, orientation, equations)
   @unpack splitting = numflux
   fm = splitting(u_rr, Val{:minus}(), orientation, equations)
   fp = splitting(u_ll, Val{:plus}(),  orientation, equations)
   return fm + fp
 end
 
+# TODO: FD. Works properly for `splitting_lax_friedrichs` but needs further
+#       testing for other splittings.
+@inline function (numflux::FluxUpwind)(u_ll, u_rr, normal_direction::AbstractVector,
+                                       equations::AbstractEquations{2})
+  @unpack splitting = numflux
+  # Compute splitting in generic normal direction with specialized
+  # eigenvalues estimates calculated inside the `splitting` function
+  f_tilde_m = splitting(u_rr, Val{:minus}(), normal_direction, equations)
+  f_tilde_p = splitting(u_ll, Val{:plus}() , normal_direction, equations)
+  return f_tilde_m + f_tilde_p
+end
+
 Base.show(io::IO, f::FluxUpwind) = print(io, "FluxUpwind(",  f.splitting, ")")
 
 
diff --git a/src/solvers/dgsem_unstructured/containers_2d.jl b/src/solvers/dgsem_unstructured/containers_2d.jl
index 340838ae5d8..31cb2228b1a 100644
--- a/src/solvers/dgsem_unstructured/containers_2d.jl
+++ b/src/solvers/dgsem_unstructured/containers_2d.jl
@@ -4,6 +4,92 @@
 # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
 @muladd begin
 
+######################################################################################################
+# TODO: FD; where should this live? Want to put it into `solvers/fdsbp_unstructured/containers_2d.jl`
+#       but then dispatching is not possible to reuse functionality for interfaces/boundaries here
+#       unless the FDSBP solver files are included before the DG files, which seems strange.
+# Container data structure (structure-of-arrays style) for FDSBP upwind solver elements on curved unstructured mesh
+struct UpwindElementContainer2D{RealT<:Real, uEltype<:Real}
+  node_coordinates           ::Array{RealT, 4}   # [ndims, nnodes, nnodes, nelement]
+  jacobian_matrix            ::Array{RealT, 5}   # [ndims, ndims, nnodes, nnodes, nelement]
+  inverse_jacobian           ::Array{RealT, 3}   # [nnodes, nnodes, nelement]
+  contravariant_vectors      ::Array{RealT, 5}   # [ndims, ndims, nnodes, nnodes, nelement]
+  contravariant_vectors_plus ::Array{RealT, 5}   # [ndims, ndims, nnodes, nnodes, nelement]
+  contravariant_vectors_minus::Array{RealT, 5}   # [ndims, ndims, nnodes, nnodes, nelement]
+  normal_directions          ::Array{RealT, 4}   # [ndims, nnodes, local sides, nelement]
+  rotations                  ::Vector{Int}       # [nelement]
+  surface_flux_values        ::Array{uEltype, 4} # [variables, nnodes, local sides, elements]
+end
+
+
+# construct an empty curved element container for the `UpwindOperators` type solver
+# to be filled later with geometries in the unstructured mesh constructor
+# OBS! Extended version of the `UnstructuredElementContainer2D` with additional arrays to hold
+#      the contravariant vectors created with the biased derivative operators.
+function UpwindElementContainer2D{RealT, uEltype}(capacity::Integer, n_variables, n_nodes) where {RealT<:Real, uEltype<:Real}
+  nan_RealT = convert(RealT, NaN)
+  nan_uEltype = convert(uEltype, NaN)
+
+  node_coordinates            = fill(nan_RealT, (2, n_nodes, n_nodes, capacity))
+  jacobian_matrix             = fill(nan_RealT, (2, 2, n_nodes, n_nodes, capacity))
+  inverse_jacobian            = fill(nan_RealT, (n_nodes, n_nodes, capacity))
+  contravariant_vectors       = fill(nan_RealT, (2, 2, n_nodes, n_nodes, capacity))
+  contravariant_vectors_plus  = fill(nan_RealT, (2, 2, n_nodes, n_nodes, capacity))
+  contravariant_vectors_minus = fill(nan_RealT, (2, 2, n_nodes, n_nodes, capacity))
+  normal_directions           = fill(nan_RealT, (2, n_nodes, 4, capacity))
+  rotations                   = fill(typemin(Int), capacity) # Fill with "nonsense" rotation values
+  surface_flux_values         = fill(nan_uEltype, (n_variables, n_nodes, 4, capacity))
+
+  return UpwindElementContainer2D{RealT, uEltype}(node_coordinates,
+                                                  jacobian_matrix,
+                                                  inverse_jacobian,
+                                                  contravariant_vectors,
+                                                  contravariant_vectors_plus,
+                                                  contravariant_vectors_minus,
+                                                  normal_directions,
+                                                  rotations,
+                                                  surface_flux_values)
+end
+
+
+@inline nelements(elements::UpwindElementContainer2D) = size(elements.surface_flux_values, 4)
+@inline eachelement(elements::UpwindElementContainer2D) = Base.OneTo(nelements(elements))
+
+@inline nvariables(elements::UpwindElementContainer2D) = size(elements.surface_flux_values, 1)
+@inline nnodes(elements::UpwindElementContainer2D) = size(elements.surface_flux_values, 2)
+
+Base.real(elements::UpwindElementContainer2D) = eltype(elements.node_coordinates)
+Base.eltype(elements::UpwindElementContainer2D) = eltype(elements.surface_flux_values)
+
+function init_elements(mesh::UnstructuredMesh2D, equations, basis::SummationByPartsOperators.UpwindOperators, RealT, uEltype)
+  elements = UpwindElementContainer2D{RealT, uEltype}(
+    mesh.n_elements, nvariables(equations), nnodes(basis))
+  init_elements!(elements, mesh, basis)
+  return elements
+end
+
+function init_elements!(elements::UpwindElementContainer2D, mesh, basis)
+  four_corners = zeros(eltype(mesh.corners), 4, 2)
+
+  # loop through elements and call the correct constructor based on whether the element is curved
+  for element in eachelement(elements)
+    if mesh.element_is_curved[element]
+      init_element!(elements, element, basis, view(mesh.surface_curves, :, element))
+    else # straight sided element
+      for i in 1:4, j in 1:2
+        # pull the (x,y) values of these corners out of the global corners array
+        four_corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]]
+      end
+      init_element!(elements, element, basis, four_corners)
+    end
+  end
+
+  # Use the mesh element information to determine the rotations, if any,
+  # of the local coordinate system in each element
+  calc_element_rotations!(elements, mesh)
+end
+
+######################################################################################################
 
 # Container data structure (structure-of-arrays style) for DG elements on curved unstructured mesh
 struct UnstructuredElementContainer2D{RealT<:Real, uEltype<:Real}
@@ -43,7 +129,7 @@ end
     eachelement(elements::UnstructuredElementContainer2D)
 
 Return an iterator over the indices that specify the location in relevant data structures
-for the elements in `elements`. 
+for the elements in `elements`.
 In particular, not the elements themselves are returned.
 """
 @inline eachelement(elements::UnstructuredElementContainer2D) = Base.OneTo(nelements(elements))
@@ -135,7 +221,7 @@ end
 @inline nnodes(interfaces::UnstructuredInterfaceContainer2D) = size(interfaces.u, 3)
 
 
-function init_interfaces(mesh::UnstructuredMesh2D, elements::UnstructuredElementContainer2D)
+function init_interfaces(mesh::UnstructuredMesh2D, elements::Union{UnstructuredElementContainer2D, UpwindElementContainer2D})
 
   interfaces = UnstructuredInterfaceContainer2D{eltype(elements)}(
     mesh.n_interfaces, nvariables(elements), nnodes(elements))
@@ -274,7 +360,7 @@ end
 @inline nboundaries(boundaries::UnstructuredBoundaryContainer2D) = length(boundaries.name)
 
 
-function init_boundaries(mesh::UnstructuredMesh2D, elements::UnstructuredElementContainer2D)
+function init_boundaries(mesh::UnstructuredMesh2D, elements::Union{UnstructuredElementContainer2D, UpwindElementContainer2D})
 
   boundaries = UnstructuredBoundaryContainer2D{real(elements), eltype(elements)}(
     mesh.n_boundaries, nvariables(elements), nnodes(elements))
diff --git a/src/solvers/fdsbp_tree/fdsbp.jl b/src/solvers/fdsbp_tree/fdsbp.jl
index c31e2c12748..b89d59c1156 100644
--- a/src/solvers/fdsbp_tree/fdsbp.jl
+++ b/src/solvers/fdsbp_tree/fdsbp.jl
@@ -36,8 +36,6 @@ end
 nnodes(D::AbstractDerivativeOperator) = size(D, 1)
 eachnode(D::AbstractDerivativeOperator) = Base.OneTo(nnodes(D))
 get_nodes(D::AbstractDerivativeOperator) = grid(D)
-get_derivative_matrix(D::DerivativeOperator) = D
-get_derivative_matrix(D::SummationByPartsOperators.UpwindOperators) = D.central
 
 # TODO: This is hack to enable the FDSBP solver to use the
 #       `SaveSolutionCallback`.
diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl
index 20a056ad64d..b3b3f07ec89 100644
--- a/src/solvers/fdsbp_unstructured/containers_2d.jl
+++ b/src/solvers/fdsbp_unstructured/containers_2d.jl
@@ -9,7 +9,8 @@
 
 # initialize all the values in the container of a general FD block (either straight sided or curved)
 # OBS! Requires the SBP derivative matrix in order to compute metric terms that are free-stream preserving
-function init_element!(elements, element, basis::AbstractDerivativeOperator, corners_or_surface_curves)
+# function init_element!(elements, element, basis::AbstractDerivativeOperator, corners_or_surface_curves)
+function init_element!(elements, element, basis::DerivativeOperator, corners_or_surface_curves)
 
   calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis), corners_or_surface_curves)
 
@@ -25,10 +26,36 @@ function init_element!(elements, element, basis::AbstractDerivativeOperator, cor
 end
 
 
+# initialize all the values in the container of a general FD block (either straight sided or curved)
+# for a set of upwind finite difference operators
+# OBS! Requires the biased derivative matrices in order to compute proper metric terms that are free-stream preserving
+function init_element!(elements, element, basis::SummationByPartsOperators.UpwindOperators, corners_or_surface_curves)
+
+  calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis), corners_or_surface_curves)
+
+  # Create the metric terms and contravariant vectors with the D⁺ operator
+  calc_metric_terms!(elements.jacobian_matrix, element, basis.plus, elements.node_coordinates)
+  calc_contravariant_vectors!(elements.contravariant_vectors_plus, element, elements.jacobian_matrix)
+
+  # Create the metric terms and contravariant vectors with the D⁻ operator
+  calc_metric_terms!(elements.jacobian_matrix, element, basis.minus, elements.node_coordinates)
+  calc_contravariant_vectors!(elements.contravariant_vectors_minus, element, elements.jacobian_matrix)
+
+  # Create the metric terms, Jacobian information, and normals for analysis
+  # and the SATs with the central D operator.
+  calc_metric_terms!(elements.jacobian_matrix, element, basis.central, elements.node_coordinates)
+  calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix)
+  calc_contravariant_vectors!(elements.contravariant_vectors, element, elements.jacobian_matrix)
+  calc_normal_directions!(elements.normal_directions, element, elements.jacobian_matrix)
+
+  return elements
+end
+
+
 # construct the metric terms for a FDSBP element "block". Directly use the derivative matrix
 # applied to the node coordinates.
 # TODO: FD; How to make this work for the upwind solver because basis has three available derivative matrices
-function calc_metric_terms!(jacobian_matrix, element, basis::AbstractDerivativeOperator, node_coordinates)
+function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeOperator, node_coordinates)
 
   # storage format:
   #   jacobian_matrix[1,1,:,:,:] <- X_xi
@@ -36,7 +63,6 @@ function calc_metric_terms!(jacobian_matrix, element, basis::AbstractDerivativeO
   #   jacobian_matrix[2,1,:,:,:] <- Y_xi
   #   jacobian_matrix[2,2,:,:,:] <- Y_eta
 
-  D_SBP = get_derivative_matrix(basis)
   # Compute the xi derivatives by applying D on the left
   # This is basically the same as
   # jacobian_matrix[1, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[1, :, :, element]
@@ -66,6 +92,7 @@ function calc_metric_terms!(jacobian_matrix, element, basis::AbstractDerivativeO
   return jacobian_matrix
 end
 
+
 # construct the normal direction vectors (but not actually normalized) for a curved sided FDSBP element "block"
 # normalization occurs on the fly during the surface flux computation
 # OBS! This assumes that the boundary points are included.
@@ -119,4 +146,42 @@ function calc_normal_directions!(normal_directions, element, jacobian_matrix)
 end
 
 
+# Compute the rotation, if any, for the local coordinate system on each `element`
+# with respect to the standard x-axis. This is necessary because if the local
+# element axis is rotated it influences the computation of the +/- directions
+# for the flux vector splitting of the upwind scheme.
+# Local principle x-axis is computed from the four corners of a particular `element`,
+# see page 35 of the "Verdict Library" for details.
+# From the local axis the `atan` function is used to determine the value of the local rotation.
+#
+# - C. J. Stimpson, C. D. Ernst, P. Knupp, P. P. Pébay, and D. Thompson (2007)
+#   The Verdict Geometric Quality Library
+#   Technical Report. Sandia National Laboraties
+#   [SAND2007-1751](https://coreform.com/papers/verdict_quality_library.pdf)
+# - `atan(y, x)` function (https://docs.julialang.org/en/v1/base/math/#Base.atan-Tuple%7BNumber%7D)
+function calc_element_rotations!(elements, mesh::UnstructuredMesh2D)
+
+  for element in 1:size(elements.node_coordinates, 4)
+    # Pull the corners for the current element
+    corner_points = mesh.corners[:, mesh.element_node_ids[:, element]]
+
+    # Compute the principle x-axis of a right-handed quadrilateral element
+    local_x_axis = (  (corner_points[:, 2] - corner_points[:, 1])
+                    + (corner_points[:, 3] - corner_points[:, 4]) )
+
+    # Two argument `atan` function retuns the angle in radians in [−pi, pi]
+    # between the positive x-axis and the point (x, y)
+    local_angle = atan(local_x_axis[2], local_x_axis[1])
+
+    # Possible reference rotations of the local axis for a given quadrilateral element in the mesh
+    #                   0°, 90°, 180°, -90° (or 270°), -180°
+    reference_angles = [0.0, 0.5*pi, pi, -0.5*pi, -pi]
+
+    # Find the nearest value and save the amount of rotation
+    elements.rotations[element] = argmin( abs.(reference_angles .- local_angle) ) - 1
+  end
+
+  return nothing
+end
+
 end # @muladd
\ No newline at end of file
diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
index 2bc0fe66225..0d0bece473b 100644
--- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
+++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl
@@ -26,12 +26,10 @@ function create_cache(mesh::UnstructuredMesh2D, equations, dg::FDSBP, RealT, uEl
 end
 
 
-# TODO: FD; Upwind versions of surface / volume integral
-
 # 2D volume integral contributions for `VolumeIntegralStrongForm`
 @inline function calc_volume_integral!(du, u,
                                        mesh::UnstructuredMesh2D,
-                                       nonconservative_terms::Val{false}, equations,
+                                       nonconservative_terms::False, equations,
                                        volume_integral::VolumeIntegralStrongForm,
                                        dg::FDSBP, cache)
   D = dg.basis # SBP derivative operator
@@ -84,6 +82,85 @@ end
 end
 
 
+# 2D volume integral contributions for `VolumeIntegralUpwind`
+# Note that the plus / minus notation of the operators does not refer to the
+# upwind / downwind directions of the fluxes.
+# Instead, the plus / minus refers to the direction of the biasing within
+# the finite difference stencils. Thus, the D^- operator acts on the positive
+# part of the flux splitting f^+ and the D^+ operator acts on the negative part
+# of the flux splitting f^-.
+# Further note that appropriate plus / minus contravariant vector terms are applied
+# to the corresponding plus / minus derivative operator.
+@inline function calc_volume_integral!(du, u,
+                                       mesh::UnstructuredMesh2D,
+                                       nonconservative_terms::False, equations,
+                                       volume_integral::VolumeIntegralUpwind,
+                                       dg::FDSBP, cache)
+  # Pull the derivative matrix assuming that
+  # dg.basis isa SummationByPartsOperators.UpwindOperators
+  # TODO: FD, improve performance to use `mul!`. Current version allocates.
+  #       Also, would be good if the logic below for local rotations could be avoided
+  D_minus = Matrix(dg.basis.minus) # Upwind SBP D^- derivative operator
+  D_plus  = Matrix(dg.basis.plus)  # Upwind SBP D^+ derivative operator
+  @unpack contravariant_vectors_plus, contravariant_vectors_minus, rotations = cache.elements
+  @unpack splitting = volume_integral
+
+  @threaded for element in eachelement(dg, cache)
+    for j in eachnode(dg), i in eachnode(dg)
+      u_node = get_node_vars(u, equations, dg, i, j, element)
+
+      # Compute the flux vector splitting. If the local coordinate
+      # system in an element is rotated, then the role of +/- swaps
+      if rotations[element] == 0
+        # standard flux splittings
+        flux1_minus, flux1_plus = splitting(u_node, 1, equations)
+        flux2_minus, flux2_plus = splitting(u_node, 2, equations)
+      elseif rotations[element] == 1 # corresponds to 90° rotation
+        # `i` direction flipped so swap x-direction
+        flux1_plus, flux1_minus = splitting(u_node, 1, equations)
+        flux2_minus, flux2_plus = splitting(u_node, 2, equations)
+      elseif rotations[element] == 2 || rotations[element] == 4 # corresponds to 180° rotation
+        # `i` and `j` directions flipped so swap x- and y-directions
+        flux1_plus, flux1_minus = splitting(u_node, 1, equations)
+        flux2_plus, flux2_minus = splitting(u_node, 2, equations)
+      else # rotations[element] == 3 # corresponds to 270° (or -90°) rotation
+        # `j` direction flipped so swap y-direction
+        flux1_minus, flux1_plus = splitting(u_node, 1, equations)
+        flux2_plus, flux2_minus = splitting(u_node, 2, equations)
+      end
+
+      # Compute the contravariant fluxes by taking the scalar product of the
+      # appropriate contravariant vector Ja^1 and the flux vector splitting
+      Ja11_plus, Ja12_plus = get_contravariant_vector(1, contravariant_vectors_plus, i, j, element)
+      contravariant_flux1_minus = Ja11_plus * flux1_minus + Ja12_plus * flux2_minus
+
+      Ja11_minus, Ja12_minus = get_contravariant_vector(1, contravariant_vectors_minus, i, j, element)
+      contravariant_flux1_plus = Ja11_minus * flux1_plus + Ja12_minus * flux2_plus
+
+      for ii in eachnode(dg)
+        multiply_add_to_node_vars!(du, D_minus[ii, i], contravariant_flux1_plus , equations, dg, ii, j, element)
+        multiply_add_to_node_vars!(du, D_plus[ii, i] , contravariant_flux1_minus, equations, dg, ii, j, element)
+      end
+
+      # Compute the contravariant fluxes by taking the scalar product of the
+      # appropriate contravariant vector Ja^2 and the flux vector splitting
+      Ja21_plus, Ja22_plus = get_contravariant_vector(2, contravariant_vectors_plus, i, j, element)
+      contravariant_flux2_minus = Ja21_plus * flux1_minus + Ja22_plus * flux2_minus
+
+      Ja21_minus, Ja22_minus = get_contravariant_vector(2, contravariant_vectors_minus, i, j, element)
+      contravariant_flux2_plus = Ja21_minus * flux1_plus + Ja22_minus * flux2_plus
+
+      for jj in eachnode(dg)
+        multiply_add_to_node_vars!(du, D_minus[jj, j], contravariant_flux2_plus , equations, dg, i, jj, element)
+        multiply_add_to_node_vars!(du, D_plus[jj, j] , contravariant_flux2_minus, equations, dg, i, jj, element)
+      end
+    end
+  end
+
+  return nothing
+end
+
+
 # Note! The local side numbering for the unstructured quadrilateral element implementation differs
 #       from the structured TreeMesh or StructuredMesh local side numbering:
 #
@@ -102,7 +179,7 @@ end
 # surface contributions are added.
 function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D,
                                 equations, surface_integral::SurfaceIntegralStrongForm,
-                                dg::DG, cache)
+                                dg::FDSBP, cache)
   inv_weight_left  = inv(left_boundary_weight(dg.basis))
   inv_weight_right = inv(right_boundary_weight(dg.basis))
   @unpack normal_directions, surface_flux_values = cache.elements

From 588453e9e0694951b66367509773c091d24ee443 Mon Sep 17 00:00:00 2001
From: Andrew Winters <andrew.ross.winters@liu.se>
Date: Mon, 16 Jan 2023 21:10:08 +0100
Subject: [PATCH 15/15] add two tests for the upwind curvilinear solver

---
 .../elixir_euler_convergence_upwind.jl        | 111 ++++++++++++++++++
 .../elixir_euler_free_stream_upwind.jl        |  78 ++++++++++++
 test/test_trixi.jl                            |   3 +
 test/test_unstructured_2d.jl                  |  16 +++
 4 files changed, 208 insertions(+)
 create mode 100644 examples/unstructured_2d_fdsbp/elixir_euler_convergence_upwind.jl
 create mode 100644 examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl

diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_convergence_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_convergence_upwind.jl
new file mode 100644
index 00000000000..261483a8edd
--- /dev/null
+++ b/examples/unstructured_2d_fdsbp/elixir_euler_convergence_upwind.jl
@@ -0,0 +1,111 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Modify the manufactured solution test to use `L = sqrt(2)` in the initial and source terms
+# such that testing works on the flipped mesh
+function initial_condition_convergence_upwind(x, t, equations::CompressibleEulerEquations2D)
+  c = 2
+  A = 0.1
+  L = sqrt(2)
+  f = 1/L
+  ω = 2 * pi * f
+  ini = c + A * sin(ω * (x[1] + x[2] - t))
+
+  rho = ini
+  rho_v1 = ini
+  rho_v2 = ini
+  rho_e = ini^2
+
+  return SVector(rho, rho_v1, rho_v2, rho_e)
+end
+
+@inline function source_terms_convergence_upwind(u, x, t, equations::CompressibleEulerEquations2D)
+  # Same settings as in `initial_condition`
+  c = 2
+  A = 0.1
+  L = sqrt(2)
+  f = 1/L
+  ω = 2 * pi * f
+  γ = equations.gamma
+
+  x1, x2 = x
+  si, co = sincos(ω * (x1 + x2 - t))
+  rho = c + A * si
+  rho_x = ω * A * co
+  # Note that d/dt rho = -d/dx rho = -d/dy rho.
+
+  tmp = (2 * rho - 1) * (γ - 1)
+
+  du1 = rho_x
+  du2 = rho_x * (1 + tmp)
+  du3 = du2
+  du4 = 2 * rho_x * (rho + tmp)
+
+  return SVector(du1, du2, du3, du4)
+end
+
+initial_condition = initial_condition_convergence_upwind
+
+###############################################################################
+# Get the FDSBP approximation operator
+
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order=1,
+                         accuracy_order=4,
+                         xmin=-1.0, xmax=1.0,
+                         N=10)
+
+flux_splitting = splitting_lax_friedrichs
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+default_mesh_file = joinpath(@__DIR__, "mesh_multiple_flips.mesh")
+isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/b434e724e3972a9c4ee48d58c80cdcdb/raw/9a967f066bc5bf081e77ef2519b3918e40a964ed/mesh_multiple_flips.mesh",
+                                       default_mesh_file)
+
+mesh_file = default_mesh_file
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity=true)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms=source_terms_convergence_upwind)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(), abstol=1.0e-9, reltol=1.0e-9,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary
diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl
new file mode 100644
index 00000000000..59e44376cd4
--- /dev/null
+++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl
@@ -0,0 +1,78 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Free-stream initial condition
+initial_condition = initial_condition_constant
+
+# Boundary conditions for free-stream testing
+boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition)
+
+boundary_conditions = Dict( :Top    => boundary_condition_free_stream,
+                            :Bottom => boundary_condition_free_stream,
+                            :Right  => boundary_condition_free_stream,
+                            :Left   => boundary_condition_free_stream )
+
+###############################################################################
+# Get the Upwind FDSBP approximation space
+
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order=1,
+                         accuracy_order=4,
+                         xmin=-1.0, xmax=1.0,
+                         N=9)
+
+flux_splitting = splitting_lax_friedrichs
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+default_mesh_file = joinpath(@__DIR__, "mesh_multiple_flips.mesh")
+isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/b434e724e3972a9c4ee48d58c80cdcdb/raw/9a967f066bc5bf081e77ef2519b3918e40a964ed/mesh_multiple_flips.mesh",
+                                       default_mesh_file)
+
+mesh_file = default_mesh_file
+
+mesh = UnstructuredMesh2D(mesh_file)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=1000,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# set small tolerances for the free-stream preservation test
+sol = solve(ode, SSPRK43(), abstol=1.0e-12, reltol=1.0e-12,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary
diff --git a/test/test_trixi.jl b/test/test_trixi.jl
index 044e421107d..c1a5d991e4c 100644
--- a/test/test_trixi.jl
+++ b/test/test_trixi.jl
@@ -153,6 +153,9 @@ macro test_nowarn_mod(expr, additional_ignore_content=String[])
           "WARNING: importing deprecated binding Colors.RGB4 into PlotUtils.\n",
           r"┌ Warning: Keyword argument letter not supported with Plots.+\n└ @ Plots.+\n",
           r"┌ Warning: `parse\(::Type, ::Coloarant\)` is deprecated.+\n│.+\n│.+\n└ @ Plots.+\n",
+          # Ignore the message about ODE Symbols nonsense. DO NOT COMMIT! This is local to my machine
+          "┌ Warning: Backwards compatability support of the new return codes to Symbols will be deprecated with the Julia v1.9 release. Please see https://docs.sciml.ai/SciMLBase/stable/interfaces/Solutions/#retcodes for more information
+└ @ SciMLBase ~/.julia/packages/SciMLBase/0IYdl/src/retcodes.jl:360\n"
         ]
         append!(ignore_content, $additional_ignore_content)
         for pattern in ignore_content
diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl
index ab15032a02e..4f68e1424d4 100644
--- a/test/test_unstructured_2d.jl
+++ b/test/test_unstructured_2d.jl
@@ -170,6 +170,22 @@ isdir(outdir) && rm(outdir, recursive=true)
       linf = [0.00015638063264633573, 0.0003329723828695563, 0.00033297238286467135, 0.0009179128246792345],
       tspan = (0.0, 0.025))
   end
+
+  # TODO: this free-stream test might change once proper boundary conditions
+  #       are implemented for the upwind solver
+  @trixi_testset "FDSBP (upwind): elixir_euler_free_stream_upwind.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "unstructured_2d_fdsbp", "elixir_euler_free_stream_upwind.jl"),
+      l2   = [1.317918039557329e-13, 1.0704743312905878e-13, 8.090840230283804e-14, 1.6278848188744677e-12],
+      linf = [6.595390900088205e-12, 6.5621397205006815e-12, 4.15131817810277e-12, 8.673950446791423e-11],
+      tspan = (0.0, 0.2))
+  end
+
+  @trixi_testset "FDSBP (upwind): elixir_euler_convergence_upwind.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "unstructured_2d_fdsbp", "elixir_euler_convergence_upwind.jl"),
+      l2   = [7.953799252662559e-5, 0.00012546042648804053, 0.00012638745384152826, 0.00040076939261004185],
+      linf = [0.0005051201525521076, 0.0008676201698107899, 0.0008294017394752107, 0.0023917887756228495],
+      tspan = (0.0, 0.025))
+  end
 end
 
 # Clean up afterwards: delete Trixi output directory