From 86d14f08bcc81a849eaa2f891f69e0545167de50 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Thu, 25 May 2023 15:43:37 -0500
Subject: [PATCH 01/46] generalize function signatures to P4estMesh

---
 src/solvers/dgsem_tree/dg_2d_parabolic.jl | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
index ca6394172ad..b7cb05856f4 100644
--- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
@@ -12,10 +12,11 @@
 #               2. compute f(u, grad(u))
 #               3. compute div(f(u, grad(u))) (i.e., the "regular" rhs! call)
 # boundary conditions will be applied to both grad(u) and div(f(u, grad(u))).
-function rhs_parabolic!(du, u, t, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic,
+function rhs_parabolic!(du, u, t, mesh::Union{TreeMesh{2}, P4estMesh{2}}, 
+                        equations_parabolic::AbstractEquationsParabolic,
                         initial_condition, boundary_conditions_parabolic, source_terms,
                         dg::DG, parabolic_scheme, cache, cache_parabolic)
-  @unpack u_transformed, gradients, flux_viscous = cache_parabolic
+  (; u_transformed, gradients, flux_viscous) = cache_parabolic
 
   # Convert conservative variables to a form more suitable for viscous flux calculations
   @trixi_timeit timer() "transform variables" transform_variables!(
@@ -85,7 +86,7 @@ end
 # Transform solution variables prior to taking the gradient
 # (e.g., conservative to primitive variables). Defaults to doing nothing.
 # TODO: can we avoid copying data?
-function transform_variables!(u_transformed, u, mesh::TreeMesh{2},
+function transform_variables!(u_transformed, u, mesh::Union{TreeMesh{2}, P4estMesh{2}},
                               equations_parabolic::AbstractEquationsParabolic,
                               dg::DG, parabolic_scheme, cache, cache_parabolic)
   @threaded for element in eachelement(dg, cache)
@@ -245,7 +246,8 @@ function prolong2boundaries!(cache_parabolic, flux_viscous,
 end
 
 
-function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::TreeMesh{2},
+function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, 
+                              mesh::Union{TreeMesh{2}, P4estMesh{2}}, 
                               equations_parabolic::AbstractEquationsParabolic,
                               dg::DG, cache, cache_parabolic)
   gradients_x, gradients_y = gradients

From b5a3ede720e2306095c94bfc0fb782783164dd27 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Thu, 25 May 2023 15:43:59 -0500
Subject: [PATCH 02/46] add specializations for P4estMesh

d
---
 src/solvers/dgsem_p4est/dg.jl              |   3 +
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 242 +++++++++++++++++++++
 2 files changed, 245 insertions(+)
 create mode 100644 src/solvers/dgsem_p4est/dg_2d_parabolic.jl

diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl
index 22f847dbf3e..dabaa896fbf 100644
--- a/src/solvers/dgsem_p4est/dg.jl
+++ b/src/solvers/dgsem_p4est/dg.jl
@@ -46,7 +46,10 @@ end
 
 
 include("containers.jl")
+
 include("dg_2d.jl")
+include("dg_2d_parabolic.jl")
+
 include("dg_3d.jl")
 include("dg_parallel.jl")
 
diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
new file mode 100644
index 00000000000..0f0a0134b25
--- /dev/null
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -0,0 +1,242 @@
+# This method is called when a SemidiscretizationHyperbolic is constructed.
+# It constructs the basic `cache` used throughout the simulation to compute
+# the RHS etc.
+function create_cache_parabolic(mesh::P4estMesh, equations_hyperbolic::AbstractEquations,
+                                equations_parabolic::AbstractEquationsParabolic,
+                                dg::DG, parabolic_scheme, RealT, uEltype)
+
+  balance!(mesh)
+
+  elements   = init_elements(mesh, equations_hyperbolic, dg.basis, uEltype)
+  interfaces = init_interfaces(mesh, equations_hyperbolic, dg.basis, elements)
+  boundaries = init_boundaries(mesh, equations_hyperbolic, dg.basis, elements)
+ 
+  n_vars = nvariables(equations_hyperbolic)
+  n_elements = nelements(elements)
+  n_nodes = nnodes(dg.basis) # nodes in one direction
+  u_transformed = Array{uEltype}(undef, n_vars, n_nodes, n_nodes, n_elements)
+  gradients = ntuple(_ -> similar(u_transformed), ndims(mesh))
+  flux_viscous = ntuple(_ -> similar(u_transformed), ndims(mesh))
+
+  cache = (; elements, interfaces, boundaries, gradients, flux_viscous, u_transformed)
+
+  return cache
+end
+
+function calc_gradient!(gradients, u_transformed, t,
+                        mesh::P4estMesh{2}, equations_parabolic,
+                        boundary_conditions_parabolic, dg::DG, 
+                        cache, cache_parabolic)
+
+  gradients_x, gradients_y = gradients
+
+  # Reset du
+  @trixi_timeit timer() "reset gradients" begin
+    reset_du!(gradients_x, dg, cache)
+    reset_du!(gradients_y, dg, cache)
+  end
+
+  # Calculate volume integral
+  @trixi_timeit timer() "volume integral" begin
+    (; derivative_dhat) = dg.basis
+    (; contravariant_vectors) = cache.elements
+    
+    @threaded for element in eachelement(dg, cache)
+
+      # Calculate gradients with respect to reference coordinates in one element
+      for j in eachnode(dg), i in eachnode(dg)
+        u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, element)
+
+        for ii in eachnode(dg)
+          multiply_add_to_node_vars!(gradients_x, derivative_dhat[ii, i], u_node, equations_parabolic, dg, ii, j, element)
+        end
+
+        for jj in eachnode(dg)
+          multiply_add_to_node_vars!(gradients_y, derivative_dhat[jj, j], u_node, equations_parabolic, dg, i, jj, element)
+        end
+      end
+
+      # now that the reference coordinate gradients are computed, transform them node-by-node to physical gradients
+      # using the contravariant vectors
+      for j in eachnode(dg), i in eachnode(dg)
+        Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
+        Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
+
+        gradients_reference_1 = get_node_vars(gradients_x, equations_parabolic, dg, i, j, element)
+        gradients_reference_2 = get_node_vars(gradients_y, equations_parabolic, dg, i, j, element)
+
+        gradient_x_node = Ja11 * gradients_reference_1 + Ja21 * gradients_reference_2
+        gradient_y_node = Ja12 * gradients_reference_1 + Ja22 * gradients_reference_2
+
+        set_node_vars!(gradients_x, gradient_x_node, equations_parabolic, dg, i, j, element)
+        set_node_vars!(gradients_y, gradient_y_node, equations_parabolic, dg, i, j, element)
+      end
+
+    end
+  end
+
+  # Prolong solution to interfaces
+  @trixi_timeit timer() "prolong2interfaces" prolong2interfaces!(
+    cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg)
+
+  # # Calculate interface fluxes
+  # @trixi_timeit timer() "interface flux" begin
+  #   (; surface_flux_values) = cache_parabolic.elements
+  #   (; neighbor_ids, orientations) = cache_parabolic.interfaces
+
+  #   @threaded for interface in eachinterface(dg, cache_parabolic)
+  #     # Get neighboring elements
+  #     left_id  = neighbor_ids[1, interface]
+  #     right_id = neighbor_ids[2, interface]
+
+  #     # Determine interface direction with respect to elements:
+  #     # orientation = 1: left -> 2, right -> 1
+  #     # orientation = 2: left -> 4, right -> 3
+  #     left_direction  = 2 * orientations[interface]
+  #     right_direction = 2 * orientations[interface] - 1
+
+  #     for i in eachnode(dg)
+  #       # Call pointwise Riemann solver
+  #       u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u,
+  #                                          equations_parabolic, dg, i, interface)
+  #       flux = 0.5 * (u_ll + u_rr)
+
+  #       # Copy flux to left and right element storage
+  #       for v in eachvariable(equations_parabolic)
+  #         surface_flux_values[v, i, left_direction,  left_id]  = flux[v]
+  #         surface_flux_values[v, i, right_direction, right_id] = flux[v]
+  #       end
+  #     end
+  #   end
+  # end
+
+  (; surface_flux_values) = cache_parabolic.elements
+  (; neighbor_ids, node_indices) = cache.interfaces
+  index_range = eachnode(dg)
+  index_end = last(index_range)
+
+  @threaded for interface in eachinterface(dg, cache)
+    # Get element and side index information on the primary element
+    primary_element = neighbor_ids[1, interface]
+    primary_indices = node_indices[1, interface]
+    primary_direction = indices2direction(primary_indices)
+
+    # Create the local i,j indexing on the primary element used to pull normal direction information
+    i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
+    j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)
+
+    i_primary = i_primary_start
+    j_primary = j_primary_start
+
+    # Get element and side index information on the secondary element
+    secondary_element = neighbor_ids[2, interface]
+    secondary_indices = node_indices[2, interface]
+    secondary_direction = indices2direction(secondary_indices)
+
+    # Initiate the secondary index to be used in the surface for loop.
+    # This index on the primary side will always run forward but
+    # the secondary index might need to run backwards for flipped sides.
+    if :i_backward in secondary_indices
+      node_secondary = index_end
+      node_secondary_step = -1
+    else
+      node_secondary = 1
+      node_secondary_step = 1
+    end
+
+    for node in eachnode(dg)
+      u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u,
+                                         equations_parabolic, dg, node, interface)                                              
+      flux = 0.5 * (u_ll + u_rr)
+
+      for v in eachvariable(equations_parabolic)
+        surface_flux_values[v, node, primary_direction, primary_element] = flux[v]
+        surface_flux_values[v, node_secondary, secondary_direction, secondary_element] = flux[v]
+      end
+  
+      # Increment primary element indices to pull the normal direction
+      i_primary += i_primary_step
+      j_primary += j_primary_step
+      # Increment the surface node index along the secondary element
+      node_secondary += node_secondary_step
+    end
+  end
+
+  # # Prolong solution to boundaries
+  # @trixi_timeit timer() "prolong2boundaries" prolong2boundaries!(
+  #   cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg)
+
+  # # Calculate boundary fluxes
+  # @trixi_timeit timer() "boundary flux" calc_boundary_flux_gradients!(
+  #   cache_parabolic, t, boundary_conditions_parabolic, mesh, equations_parabolic,
+  #   dg.surface_integral, dg)
+
+  # TODO: parabolic; mortars
+
+  # Calculate surface integrals
+  @trixi_timeit timer() "surface integral" begin
+    (; boundary_interpolation) = dg.basis
+    (; surface_flux_values) = cache_parabolic.elements
+    (; contravariant_vectors) = cache.elements
+
+    # Note that all fluxes have been computed with outward-pointing normal vectors.
+    # Access the factors only once before beginning the loop to increase performance.
+    # We also use explicit assignments instead of `+=` to let `@muladd` turn these
+    # into FMAs (see comment at the top of the file).
+    factor_1 = boundary_interpolation[1,          1]
+    factor_2 = boundary_interpolation[nnodes(dg), 2]
+    @threaded for element in eachelement(dg, cache)
+      for l in eachnode(dg)
+        for v in eachvariable(equations_parabolic)
+          # surface at -x
+          normal_direction_x, _ = get_normal_direction(1, contravariant_vectors, l, 1, element)
+          gradients_x[v, 1,          l, element] = (
+            gradients_x[v, 1,          l, element] + surface_flux_values[v, l, 1, element] * factor_1 * normal_direction_x)
+
+          # surface at +x
+          normal_direction_x, _ = get_normal_direction(2, contravariant_vectors, l, 1, element)
+          gradients_x[v, nnodes(dg), l, element] = (
+            gradients_x[v, nnodes(dg), l, element] + surface_flux_values[v, l, 2, element] * factor_2 * normal_direction_x)
+
+          # surface at -y
+          _, normal_direction_y = get_normal_direction(3, contravariant_vectors, l, 1, element)
+          gradients_y[v, l, 1,          element] = (
+            gradients_y[v, l, 1,          element] + surface_flux_values[v, l, 3, element] * factor_1 * normal_direction_y)
+
+          # surface at +y
+          _, normal_direction_y = get_normal_direction(4, contravariant_vectors, l, 1, element)
+          gradients_y[v, l, nnodes(dg), element] = (
+            gradients_y[v, l, nnodes(dg), element] + surface_flux_values[v, l, 4, element] * factor_2 * normal_direction_y)
+        end
+      end
+    end
+  end
+
+  # Apply Jacobian from mapping to reference element
+  @trixi_timeit timer() "Jacobian" begin
+    apply_jacobian!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic)
+  end
+
+  return nothing
+end
+
+# Needed to *not* flip the sign of the inverse Jacobian.
+# This is because the parabolic fluxes are assumed to be of the form
+#   `du/dt + df/dx = dg/dx + source(x,t)`,
+# where f(u) is the inviscid flux and g(u) is the viscous flux.
+function apply_jacobian!(du, mesh::TreeMesh{2},
+                         equations::AbstractEquationsParabolic, dg::DG, cache)
+
+  @threaded for element in eachelement(dg, cache)
+    factor = cache.elements.inverse_jacobian[element]
+
+    for j in eachnode(dg), i in eachnode(dg)
+      for v in eachvariable(equations)
+        du[v, i, j, element] *= factor
+      end
+    end
+  end
+
+  return nothing
+end
\ No newline at end of file

From 68ac9e1e10107187c396fd588705d9f671f0e85f Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Thu, 25 May 2023 15:56:44 -0500
Subject: [PATCH 03/46] add normals

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 0f0a0134b25..bc897845a6c 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -189,12 +189,12 @@ function calc_gradient!(gradients, u_transformed, t,
       for l in eachnode(dg)
         for v in eachvariable(equations_parabolic)
           # surface at -x
-          normal_direction_x, _ = get_normal_direction(1, contravariant_vectors, l, 1, element)
+          normal_direction_x, _ = get_normal_direction(1, contravariant_vectors, 1, l, element)
           gradients_x[v, 1,          l, element] = (
             gradients_x[v, 1,          l, element] + surface_flux_values[v, l, 1, element] * factor_1 * normal_direction_x)
 
           # surface at +x
-          normal_direction_x, _ = get_normal_direction(2, contravariant_vectors, l, 1, element)
+          normal_direction_x, _ = get_normal_direction(2, contravariant_vectors, nnodes(dg), l, element)
           gradients_x[v, nnodes(dg), l, element] = (
             gradients_x[v, nnodes(dg), l, element] + surface_flux_values[v, l, 2, element] * factor_2 * normal_direction_x)
 
@@ -204,7 +204,7 @@ function calc_gradient!(gradients, u_transformed, t,
             gradients_y[v, l, 1,          element] + surface_flux_values[v, l, 3, element] * factor_1 * normal_direction_y)
 
           # surface at +y
-          _, normal_direction_y = get_normal_direction(4, contravariant_vectors, l, 1, element)
+          _, normal_direction_y = get_normal_direction(4, contravariant_vectors, l, nnodes(dg), element)
           gradients_y[v, l, nnodes(dg), element] = (
             gradients_y[v, l, nnodes(dg), element] + surface_flux_values[v, l, 4, element] * factor_2 * normal_direction_y)
         end

From 4064f89a8281cf5c11029934adfb8fe166aa552f Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Thu, 25 May 2023 15:57:30 -0500
Subject: [PATCH 04/46] add surface integrals

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 126 ++++++++-------------
 1 file changed, 49 insertions(+), 77 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index bc897845a6c..271b3b6c419 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -79,86 +79,58 @@ function calc_gradient!(gradients, u_transformed, t,
   @trixi_timeit timer() "prolong2interfaces" prolong2interfaces!(
     cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg)
 
-  # # Calculate interface fluxes
-  # @trixi_timeit timer() "interface flux" begin
-  #   (; surface_flux_values) = cache_parabolic.elements
-  #   (; neighbor_ids, orientations) = cache_parabolic.interfaces
-
-  #   @threaded for interface in eachinterface(dg, cache_parabolic)
-  #     # Get neighboring elements
-  #     left_id  = neighbor_ids[1, interface]
-  #     right_id = neighbor_ids[2, interface]
-
-  #     # Determine interface direction with respect to elements:
-  #     # orientation = 1: left -> 2, right -> 1
-  #     # orientation = 2: left -> 4, right -> 3
-  #     left_direction  = 2 * orientations[interface]
-  #     right_direction = 2 * orientations[interface] - 1
-
-  #     for i in eachnode(dg)
-  #       # Call pointwise Riemann solver
-  #       u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u,
-  #                                          equations_parabolic, dg, i, interface)
-  #       flux = 0.5 * (u_ll + u_rr)
-
-  #       # Copy flux to left and right element storage
-  #       for v in eachvariable(equations_parabolic)
-  #         surface_flux_values[v, i, left_direction,  left_id]  = flux[v]
-  #         surface_flux_values[v, i, right_direction, right_id] = flux[v]
-  #       end
-  #     end
-  #   end
-  # end
-
-  (; surface_flux_values) = cache_parabolic.elements
-  (; neighbor_ids, node_indices) = cache.interfaces
-  index_range = eachnode(dg)
-  index_end = last(index_range)
-
-  @threaded for interface in eachinterface(dg, cache)
-    # Get element and side index information on the primary element
-    primary_element = neighbor_ids[1, interface]
-    primary_indices = node_indices[1, interface]
-    primary_direction = indices2direction(primary_indices)
-
-    # Create the local i,j indexing on the primary element used to pull normal direction information
-    i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
-    j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)
-
-    i_primary = i_primary_start
-    j_primary = j_primary_start
-
-    # Get element and side index information on the secondary element
-    secondary_element = neighbor_ids[2, interface]
-    secondary_indices = node_indices[2, interface]
-    secondary_direction = indices2direction(secondary_indices)
-
-    # Initiate the secondary index to be used in the surface for loop.
-    # This index on the primary side will always run forward but
-    # the secondary index might need to run backwards for flipped sides.
-    if :i_backward in secondary_indices
-      node_secondary = index_end
-      node_secondary_step = -1
-    else
-      node_secondary = 1
-      node_secondary_step = 1
-    end
+  # Calculate interface fluxes
+  @trixi_timeit timer() "interface flux" begin
+    (; surface_flux_values) = cache_parabolic.elements
+    (; neighbor_ids, node_indices) = cache.interfaces
+    index_range = eachnode(dg)
+    index_end = last(index_range)
+
+    @threaded for interface in eachinterface(dg, cache)
+      # Get element and side index information on the primary element
+      primary_element = neighbor_ids[1, interface]
+      primary_indices = node_indices[1, interface]
+      primary_direction = indices2direction(primary_indices)
+
+      # Create the local i,j indexing on the primary element used to pull normal direction information
+      i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
+      j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)
+
+      i_primary = i_primary_start
+      j_primary = j_primary_start
+
+      # Get element and side index information on the secondary element
+      secondary_element = neighbor_ids[2, interface]
+      secondary_indices = node_indices[2, interface]
+      secondary_direction = indices2direction(secondary_indices)
+
+      # Initiate the secondary index to be used in the surface for loop.
+      # This index on the primary side will always run forward but
+      # the secondary index might need to run backwards for flipped sides.
+      if :i_backward in secondary_indices
+        node_secondary = index_end
+        node_secondary_step = -1
+      else
+        node_secondary = 1
+        node_secondary_step = 1
+      end
 
-    for node in eachnode(dg)
-      u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u,
-                                         equations_parabolic, dg, node, interface)                                              
-      flux = 0.5 * (u_ll + u_rr)
+      for node in eachnode(dg)
+        u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u,
+                                          equations_parabolic, dg, node, interface)                                              
+        flux = 0.5 * (u_ll + u_rr)
 
-      for v in eachvariable(equations_parabolic)
-        surface_flux_values[v, node, primary_direction, primary_element] = flux[v]
-        surface_flux_values[v, node_secondary, secondary_direction, secondary_element] = flux[v]
+        for v in eachvariable(equations_parabolic)
+          surface_flux_values[v, node, primary_direction, primary_element] = flux[v]
+          surface_flux_values[v, node_secondary, secondary_direction, secondary_element] = flux[v]
+        end
+    
+        # Increment primary element indices to pull the normal direction
+        i_primary += i_primary_step
+        j_primary += j_primary_step
+        # Increment the surface node index along the secondary element
+        node_secondary += node_secondary_step
       end
-  
-      # Increment primary element indices to pull the normal direction
-      i_primary += i_primary_step
-      j_primary += j_primary_step
-      # Increment the surface node index along the secondary element
-      node_secondary += node_secondary_step
     end
   end
 

From 32cb2c13a042e58e3f7d1f419b31865a1c53acfe Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Thu, 25 May 2023 16:01:55 -0500
Subject: [PATCH 05/46] fix type ambiguity

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

diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl
index a8972dfe766..6239a6572fe 100644
--- a/src/solvers/dgsem_structured/dg_2d.jl
+++ b/src/solvers/dgsem_structured/dg_2d.jl
@@ -544,7 +544,7 @@ end
 
 function apply_jacobian!(du,
                          mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}},
-                         equations, dg::DG, cache)
+                         equations::AbstractEquations, dg::DG, cache)
   @unpack inverse_jacobian = cache.elements
 
   @threaded for element in eachelement(dg, cache)

From ce00a0c7a83cee18206e9c259d7441cc0be59fe9 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Thu, 25 May 2023 16:02:35 -0500
Subject: [PATCH 06/46] generalizing `apply_jacobian!` to P4estMesh

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

diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
index b7cb05856f4..7decee18fd2 100644
--- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
@@ -588,7 +588,7 @@ end
 # This is because the parabolic fluxes are assumed to be of the form
 #   `du/dt + df/dx = dg/dx + source(x,t)`,
 # where f(u) is the inviscid flux and g(u) is the viscous flux.
-function apply_jacobian!(du, mesh::TreeMesh{2},
+function apply_jacobian!(du, mesh::Union{TreeMesh{2}, P4estMesh{2}}, 
                          equations::AbstractEquationsParabolic, dg::DG, cache)
 
   @threaded for element in eachelement(dg, cache)

From ca9d485a1e604d1b0c34b68ab6bd38d216da3b44 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Thu, 25 May 2023 16:21:22 -0500
Subject: [PATCH 07/46] resolving type ambiguity with apply_jacobian!

d
---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 156 +++++++++++++++++++--
 src/solvers/dgsem_structured/dg_2d.jl      |   2 +-
 src/solvers/dgsem_tree/dg_2d_parabolic.jl  |   8 +-
 3 files changed, 148 insertions(+), 18 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 271b3b6c419..a8c6f69004b 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -186,29 +186,159 @@ function calc_gradient!(gradients, u_transformed, t,
 
   # Apply Jacobian from mapping to reference element
   @trixi_timeit timer() "Jacobian" begin
-    apply_jacobian!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic)
-    apply_jacobian!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic)
   end
 
   return nothing
 end
 
-# Needed to *not* flip the sign of the inverse Jacobian.
-# This is because the parabolic fluxes are assumed to be of the form
-#   `du/dt + df/dx = dg/dx + source(x,t)`,
-# where f(u) is the inviscid flux and g(u) is the viscous flux.
-function apply_jacobian!(du, mesh::TreeMesh{2},
-                         equations::AbstractEquationsParabolic, dg::DG, cache)
+# This is the version used when calculating the divergence of the viscous fluxes
+function calc_volume_integral!(du, flux_viscous,
+                               mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic,
+                               dg::DGSEM, cache)
+  (; derivative_dhat) = dg.basis
+  (; contravariant_vectors) = cache.elements
+  flux_viscous_x, flux_viscous_y = flux_viscous
 
   @threaded for element in eachelement(dg, cache)
-    factor = cache.elements.inverse_jacobian[element]
-
+    # Calculate volume terms in one element
     for j in eachnode(dg), i in eachnode(dg)
-      for v in eachvariable(equations)
-        du[v, i, j, element] *= factor
+      flux1 = get_node_vars(flux_viscous_x, equations_parabolic, dg, i, j, element)
+      flux2 = get_node_vars(flux_viscous_y, equations_parabolic, dg, i, j, element)
+  
+      # 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, derivative_dhat[ii, i], contravariant_flux1, equations_parabolic, 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, derivative_dhat[jj, j], contravariant_flux2, equations_parabolic, dg, i, jj, element)
+      end
+    end
+  end
+
+  return nothing
+end
+
+
+# This is the version used when calculating the divergence of the viscous fluxes
+# We pass the `surface_integral` argument solely for dispatch
+function prolong2interfaces!(cache_parabolic, flux_viscous,
+                             mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic,
+                             surface_integral, dg::DG, cache)
+  (; interfaces) = cache
+  index_range = eachnode(dg)
+  flux_viscous_x, flux_viscous_y = flux_viscous
+
+  @threaded for interface in eachinterface(dg, cache)
+    # Copy solution data from the primary element using "delayed indexing" with
+    # a start value and a step size to get the correct face and orientation.
+    # Note that in the current implementation, the interface will be
+    # "aligned at the primary element", i.e., the index of the primary side
+    # will always run forwards.
+    primary_element = interfaces.neighbor_ids[1, interface]
+    primary_indices = interfaces.node_indices[1, interface]
+
+    i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
+    j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)
+
+    i_primary = i_primary_start
+    j_primary = j_primary_start
+    for i in eachnode(dg)
+      for v in eachvariable(equations_parabolic)
+        # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
+        interfaces.u[1, v, i, interface] = flux_viscous_x[v, i_primary, j_primary, primary_element]
+      end
+      i_primary += i_primary_step
+      j_primary += j_primary_step
+    end
+
+    # Copy solution data from the secondary element using "delayed indexing" with
+    # a start value and a step size to get the correct face and orientation.
+    secondary_element = interfaces.neighbor_ids[2, interface]
+    secondary_indices = interfaces.node_indices[2, interface]
+
+    i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], index_range)
+    j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], index_range)
+
+    i_secondary = i_secondary_start
+    j_secondary = j_secondary_start
+    for i in eachnode(dg)
+      for v in eachvariable(equations_parabolic)
+        # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
+        interfaces.u[2, v, i, interface] = flux_viscous_y[v, i_secondary, j_secondary, secondary_element]
       end
+      i_secondary += i_secondary_step
+      j_secondary += j_secondary_step
     end
   end
 
   return nothing
-end
\ No newline at end of file
+end
+
+# TODO: Deal with orientations for the BR1 divergence flux = {σ} ⋅ normal
+function calc_interface_flux!(surface_flux_values,
+                              mesh::P4estMesh{2}, equations_parabolic,
+                              dg::DG, cache_parabolic)
+
+  (; neighbor_ids, node_indices) = cache.interfaces
+  (; contravariant_vectors) = cache.elements
+  index_range = eachnode(dg)
+  index_end = last(index_range)
+
+  @threaded for interface in eachinterface(dg, cache)
+    # Get element and side index information on the primary element
+    primary_element = neighbor_ids[1, interface]
+    primary_indices = node_indices[1, interface]
+    primary_direction = indices2direction(primary_indices)
+
+    # Create the local i,j indexing on the primary element used to pull normal direction information
+    i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
+    j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)
+
+    i_primary = i_primary_start
+    j_primary = j_primary_start
+
+    # Get element and side index information on the secondary element
+    secondary_element = neighbor_ids[2, interface]
+    secondary_indices = node_indices[2, interface]
+    secondary_direction = indices2direction(secondary_indices)
+
+    # Initiate the secondary index to be used in the surface for loop.
+    # This index on the primary side will always run forward but
+    # the secondary index might need to run backwards for flipped sides.
+    if :i_backward in secondary_indices
+      node_secondary = index_end
+      node_secondary_step = -1
+    else
+      node_secondary = 1
+      node_secondary_step = 1
+    end
+
+    for node in eachnode(dg)
+      # Get the normal direction on the primary element.
+      # Contravariant vectors at interfaces in negative coordinate direction
+      # are pointing inwards. This is handled by `get_normal_direction`.
+      normal_direction = get_normal_direction(primary_direction, contravariant_vectors,
+                                              i_primary, j_primary, primary_element)
+
+      
+
+      # Increment primary element indices to pull the normal direction
+      i_primary += i_primary_step
+      j_primary += j_primary_step
+      # Increment the surface node index along the secondary element
+      node_secondary += node_secondary_step
+    end
+  end
+
+  return nothing
+end
diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl
index 6239a6572fe..a8972dfe766 100644
--- a/src/solvers/dgsem_structured/dg_2d.jl
+++ b/src/solvers/dgsem_structured/dg_2d.jl
@@ -544,7 +544,7 @@ end
 
 function apply_jacobian!(du,
                          mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}},
-                         equations::AbstractEquations, dg::DG, cache)
+                         equations, dg::DG, cache)
   @unpack inverse_jacobian = cache.elements
 
   @threaded for element in eachelement(dg, cache)
diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
index 7decee18fd2..690e5cbd297 100644
--- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
@@ -542,8 +542,8 @@ function calc_gradient!(gradients, u_transformed, t,
 
   # Apply Jacobian from mapping to reference element
   @trixi_timeit timer() "Jacobian" begin
-    apply_jacobian!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic)
-    apply_jacobian!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic)
   end
 
   return nothing
@@ -588,8 +588,8 @@ end
 # This is because the parabolic fluxes are assumed to be of the form
 #   `du/dt + df/dx = dg/dx + source(x,t)`,
 # where f(u) is the inviscid flux and g(u) is the viscous flux.
-function apply_jacobian!(du, mesh::Union{TreeMesh{2}, P4estMesh{2}}, 
-                         equations::AbstractEquationsParabolic, dg::DG, cache)
+function apply_jacobian_parabolic!(du, mesh::Union{TreeMesh{2}, P4estMesh{2}}, 
+                                   equations::AbstractEquationsParabolic, dg::DG, cache)
 
   @threaded for element in eachelement(dg, cache)
     factor = cache.elements.inverse_jacobian[element]

From 70cda42554c7657faa14cd26d1f52a6573744f58 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 09:33:06 -0500
Subject: [PATCH 08/46] `apply_jacobian!` -> `apply_jacobian_parabolic!`

---
 src/solvers/dgsem_tree/dg_1d_parabolic.jl |  4 ++--
 src/solvers/dgsem_tree/dg_3d_parabolic.jl | 12 ++++++------
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl
index be4235c627b..c53d70f2a48 100644
--- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl
@@ -73,7 +73,7 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, equations_parabolic::Abstra
     du, u, mesh, equations_parabolic, dg.surface_integral, dg, cache_parabolic)
 
   # Apply Jacobian from mapping to reference element
-  @trixi_timeit timer() "Jacobian" apply_jacobian!(
+  @trixi_timeit timer() "Jacobian" apply_jacobian_parabolic!(
     du, mesh, equations_parabolic, dg, cache_parabolic)
 
   return nothing
@@ -483,7 +483,7 @@ end
 # This is because the parabolic fluxes are assumed to be of the form
 #   `du/dt + df/dx = dg/dx + source(x,t)`,
 # where f(u) is the inviscid flux and g(u) is the viscous flux.
-function apply_jacobian!(du, mesh::TreeMesh{1},
+function apply_jacobian_parabolic!(du, mesh::TreeMesh{1},
                          equations::AbstractEquationsParabolic, dg::DG, cache)
 
   @threaded for element in eachelement(dg, cache)
diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl
index d3a47cb06be..d4a197de172 100644
--- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl
@@ -76,7 +76,7 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{3}, equations_parabolic::Abstra
     du, u, mesh, equations_parabolic, dg.surface_integral, dg, cache_parabolic)
 
   # Apply Jacobian from mapping to reference element
-  @trixi_timeit timer() "Jacobian" apply_jacobian!(
+  @trixi_timeit timer() "Jacobian" apply_jacobian_parabolic!(
     du, mesh, equations_parabolic, dg, cache_parabolic)
 
   return nothing
@@ -601,9 +601,9 @@ function calc_gradient!(gradients, u_transformed, t,
 
   # Apply Jacobian from mapping to reference element
   @trixi_timeit timer() "Jacobian" begin
-    apply_jacobian!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic)
-    apply_jacobian!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic)
-    apply_jacobian!(gradients_z, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_z, mesh, equations_parabolic, dg, cache_parabolic)
   end
 
   return nothing
@@ -648,8 +648,8 @@ end
 # This is because the parabolic fluxes are assumed to be of the form
 #   `du/dt + df/dx = dg/dx + source(x,t)`,
 # where f(u) is the inviscid flux and g(u) is the viscous flux.
-function apply_jacobian!(du, mesh::TreeMesh{3},
-                         equations::AbstractEquationsParabolic, dg::DG, cache)
+function apply_jacobian_parabolic!(du, mesh::TreeMesh{3},
+                                   equations::AbstractEquationsParabolic, dg::DG, cache)
 
   @threaded for element in eachelement(dg, cache)
     factor = cache.elements.inverse_jacobian[element]

From 82b4b5e33a44abddf2f48f89d23dbb7795e8ea47 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 09:33:06 -0500
Subject: [PATCH 09/46] `apply_jacobian!` -> `apply_jacobian_parabolic!`

---
 src/solvers/dgsem_tree/dg_1d_parabolic.jl |  4 ++--
 src/solvers/dgsem_tree/dg_3d_parabolic.jl | 12 ++++++------
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl
index be4235c627b..c53d70f2a48 100644
--- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl
@@ -73,7 +73,7 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, equations_parabolic::Abstra
     du, u, mesh, equations_parabolic, dg.surface_integral, dg, cache_parabolic)
 
   # Apply Jacobian from mapping to reference element
-  @trixi_timeit timer() "Jacobian" apply_jacobian!(
+  @trixi_timeit timer() "Jacobian" apply_jacobian_parabolic!(
     du, mesh, equations_parabolic, dg, cache_parabolic)
 
   return nothing
@@ -483,7 +483,7 @@ end
 # This is because the parabolic fluxes are assumed to be of the form
 #   `du/dt + df/dx = dg/dx + source(x,t)`,
 # where f(u) is the inviscid flux and g(u) is the viscous flux.
-function apply_jacobian!(du, mesh::TreeMesh{1},
+function apply_jacobian_parabolic!(du, mesh::TreeMesh{1},
                          equations::AbstractEquationsParabolic, dg::DG, cache)
 
   @threaded for element in eachelement(dg, cache)
diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl
index d3a47cb06be..d4a197de172 100644
--- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl
@@ -76,7 +76,7 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{3}, equations_parabolic::Abstra
     du, u, mesh, equations_parabolic, dg.surface_integral, dg, cache_parabolic)
 
   # Apply Jacobian from mapping to reference element
-  @trixi_timeit timer() "Jacobian" apply_jacobian!(
+  @trixi_timeit timer() "Jacobian" apply_jacobian_parabolic!(
     du, mesh, equations_parabolic, dg, cache_parabolic)
 
   return nothing
@@ -601,9 +601,9 @@ function calc_gradient!(gradients, u_transformed, t,
 
   # Apply Jacobian from mapping to reference element
   @trixi_timeit timer() "Jacobian" begin
-    apply_jacobian!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic)
-    apply_jacobian!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic)
-    apply_jacobian!(gradients_z, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_z, mesh, equations_parabolic, dg, cache_parabolic)
   end
 
   return nothing
@@ -648,8 +648,8 @@ end
 # This is because the parabolic fluxes are assumed to be of the form
 #   `du/dt + df/dx = dg/dx + source(x,t)`,
 # where f(u) is the inviscid flux and g(u) is the viscous flux.
-function apply_jacobian!(du, mesh::TreeMesh{3},
-                         equations::AbstractEquationsParabolic, dg::DG, cache)
+function apply_jacobian_parabolic!(du, mesh::TreeMesh{3},
+                                   equations::AbstractEquationsParabolic, dg::DG, cache)
 
   @threaded for element in eachelement(dg, cache)
     factor = cache.elements.inverse_jacobian[element]

From dd36650a69658f909cc4b2374d8b27852b11d721 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 09:40:01 -0500
Subject: [PATCH 10/46] switch to `apply_jacobian_parabolic!`

---
 src/solvers/dgsem_tree/dg_2d_parabolic.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
index ca6394172ad..97bbc7e2633 100644
--- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
@@ -76,7 +76,7 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{2}, equations_parabolic::Abstra
     du, u, mesh, equations_parabolic, dg.surface_integral, dg, cache_parabolic)
 
   # Apply Jacobian from mapping to reference element
-  @trixi_timeit timer() "Jacobian" apply_jacobian!(
+  @trixi_timeit timer() "Jacobian" apply_jacobian_parabolic!(
     du, mesh, equations_parabolic, dg, cache_parabolic)
 
   return nothing
@@ -540,8 +540,8 @@ function calc_gradient!(gradients, u_transformed, t,
 
   # Apply Jacobian from mapping to reference element
   @trixi_timeit timer() "Jacobian" begin
-    apply_jacobian!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic)
-    apply_jacobian!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_x, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic)
   end
 
   return nothing
@@ -586,8 +586,8 @@ end
 # This is because the parabolic fluxes are assumed to be of the form
 #   `du/dt + df/dx = dg/dx + source(x,t)`,
 # where f(u) is the inviscid flux and g(u) is the viscous flux.
-function apply_jacobian!(du, mesh::TreeMesh{2},
-                         equations::AbstractEquationsParabolic, dg::DG, cache)
+function apply_jacobian_parabolic!(du, mesh::TreeMesh{2},
+                                   equations::AbstractEquationsParabolic, dg::DG, cache)
 
   @threaded for element in eachelement(dg, cache)
     factor = cache.elements.inverse_jacobian[element]

From 7d0bf071c92a622b5b21242fcff3406861ce4420 Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Fri, 26 May 2023 10:27:09 -0500
Subject: [PATCH 11/46] Update src/solvers/dgsem_tree/dg_1d_parabolic.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/solvers/dgsem_tree/dg_1d_parabolic.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl
index c53d70f2a48..0e044b3239c 100644
--- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl
@@ -484,7 +484,7 @@ end
 #   `du/dt + df/dx = dg/dx + source(x,t)`,
 # where f(u) is the inviscid flux and g(u) is the viscous flux.
 function apply_jacobian_parabolic!(du, mesh::TreeMesh{1},
-                         equations::AbstractEquationsParabolic, dg::DG, cache)
+                                   equations::AbstractEquationsParabolic, dg::DG, cache)
 
   @threaded for element in eachelement(dg, cache)
     factor = cache.elements.inverse_jacobian[element]

From b74ec449b37502e3fd77f0036b2c3c69a3400d8f Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 12:03:01 -0500
Subject: [PATCH 12/46] missed one

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

diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl
index 0e044b3239c..1bec34568d8 100644
--- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl
@@ -444,7 +444,7 @@ function calc_gradient!(gradients, u_transformed, t,
 
   # Apply Jacobian from mapping to reference element
   @trixi_timeit timer() "Jacobian" begin
-    apply_jacobian!(gradients, mesh, equations_parabolic, dg, cache_parabolic)
+    apply_jacobian_parabolic!(gradients, mesh, equations_parabolic, dg, cache_parabolic)
   end
 
   return nothing

From 3eb7bf5017fdfbb39f383d386eb98824be7ac36d Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 12:35:03 -0500
Subject: [PATCH 13/46] draft of prolong2interfaces and calc_interface_flux

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 43 +++++++++++++++++-----
 1 file changed, 34 insertions(+), 9 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index a8c6f69004b..5204aa1fe8b 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -246,6 +246,7 @@ function prolong2interfaces!(cache_parabolic, flux_viscous,
     # will always run forwards.
     primary_element = interfaces.neighbor_ids[1, interface]
     primary_indices = interfaces.node_indices[1, interface]
+    primary_direction = indices2direction(primary_indices)
 
     i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
     j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)
@@ -253,9 +254,17 @@ function prolong2interfaces!(cache_parabolic, flux_viscous,
     i_primary = i_primary_start
     j_primary = j_primary_start
     for i in eachnode(dg)
+
+      # this is the outward normal direction on the primary element
+      normal_direction = get_normal_direction(primary_direction, contravariant_vectors,
+                                              i_primary, j_primary, primary_element)
+
       for v in eachvariable(equations_parabolic)
         # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
-        interfaces.u[1, v, i, interface] = flux_viscous_x[v, i_primary, j_primary, primary_element]
+        flux_viscous = SVector(flux_viscous_x[v, i_primary, j_primary, primary_element], 
+                               flux_viscous_y[v, i_primary, j_primary, primary_element])
+
+        interfaces.u[1, v, i, interface] = dot(flux_viscous, normal_direction)
       end
       i_primary += i_primary_step
       j_primary += j_primary_step
@@ -265,6 +274,7 @@ function prolong2interfaces!(cache_parabolic, flux_viscous,
     # a start value and a step size to get the correct face and orientation.
     secondary_element = interfaces.neighbor_ids[2, interface]
     secondary_indices = interfaces.node_indices[2, interface]
+    secondary_direction = indices2direction(secondary_indices)
 
     i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], index_range)
     j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], index_range)
@@ -272,9 +282,18 @@ function prolong2interfaces!(cache_parabolic, flux_viscous,
     i_secondary = i_secondary_start
     j_secondary = j_secondary_start
     for i in eachnode(dg)
+      # This is the outward normal direction on the secondary element.
+      # Here, we assume that normal_direction on the secondary element is 
+      # the negative of normal_direction on the primary element.  
+      normal_direction = get_normal_direction(secondary_direction, contravariant_vectors,
+                                              i_secondary, j_secondary, secondary_element)
+
       for v in eachvariable(equations_parabolic)
         # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
-        interfaces.u[2, v, i, interface] = flux_viscous_y[v, i_secondary, j_secondary, secondary_element]
+        flux_viscous = SVector(flux_viscous_x[v, i_secondary, j_secondary, secondary_element], 
+                               flux_viscous_y[v, i_secondary, j_secondary, secondary_element])
+        # store the normal flux with respect to the primary normal direction
+        interfaces.u[2, v, i, interface] = -dot(flux_viscous, normal_direction)
       end
       i_secondary += i_secondary_step
       j_secondary += j_secondary_step
@@ -323,14 +342,20 @@ function calc_interface_flux!(surface_flux_values,
       node_secondary_step = 1
     end
 
-    for node in eachnode(dg)
-      # Get the normal direction on the primary element.
-      # Contravariant vectors at interfaces in negative coordinate direction
-      # are pointing inwards. This is handled by `get_normal_direction`.
-      normal_direction = get_normal_direction(primary_direction, contravariant_vectors,
-                                              i_primary, j_primary, primary_element)
+    for i in eachnode(dg)
+            
+      flux_viscous_normal_ll, flux_viscous_normal_rr = 
+        get_surface_node_vars(cache_parabolic.interfaces.u, 
+                              equations_parabolic, dg, i, interface)
 
-      
+      # TODO: parabolic; only BR1 at the moment
+      flux = 0.5 * (flux_viscous_normal_ll + flux_viscous_normal_rr)
+
+      # Copy flux to left and right element storage
+      for v in eachvariable(equations_parabolic)
+        surface_flux_values[v, i, left_direction,  left_id]  = flux[v]
+        surface_flux_values[v, i, right_direction, right_id] = -flux[v]
+      end
 
       # Increment primary element indices to pull the normal direction
       i_primary += i_primary_step

From 984fde1c91ac6426ac17a58306e732f6565c985f Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 13:58:19 -0500
Subject: [PATCH 14/46] cache -> cache_parabolic

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 5204aa1fe8b..7ced23b0fc1 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -234,7 +234,8 @@ end
 function prolong2interfaces!(cache_parabolic, flux_viscous,
                              mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic,
                              surface_integral, dg::DG, cache)
-  (; interfaces) = cache
+  (; interfaces) = cache_parabolic
+  (; contravariant_vectors) = cache_parabolic.elements
   index_range = eachnode(dg)
   flux_viscous_x, flux_viscous_y = flux_viscous
 
@@ -308,12 +309,12 @@ function calc_interface_flux!(surface_flux_values,
                               mesh::P4estMesh{2}, equations_parabolic,
                               dg::DG, cache_parabolic)
 
-  (; neighbor_ids, node_indices) = cache.interfaces
-  (; contravariant_vectors) = cache.elements
+  (; neighbor_ids, node_indices) = cache_parabolic.interfaces
+  (; contravariant_vectors) = cache_parabolic.elements
   index_range = eachnode(dg)
   index_end = last(index_range)
 
-  @threaded for interface in eachinterface(dg, cache)
+  @threaded for interface in eachinterface(dg, cache_parabolic)
     # Get element and side index information on the primary element
     primary_element = neighbor_ids[1, interface]
     primary_indices = node_indices[1, interface]

From 0ad7cf75b3bb965655e3c55ba76f14d16a8602ab Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 13:59:01 -0500
Subject: [PATCH 15/46] adding prolong2boundaries! and
 calc_boundary_flux_gradients! back

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

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 7ced23b0fc1..b15e63f6cfb 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -134,14 +134,14 @@ function calc_gradient!(gradients, u_transformed, t,
     end
   end
 
-  # # Prolong solution to boundaries
-  # @trixi_timeit timer() "prolong2boundaries" prolong2boundaries!(
-  #   cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg)
-
-  # # Calculate boundary fluxes
-  # @trixi_timeit timer() "boundary flux" calc_boundary_flux_gradients!(
-  #   cache_parabolic, t, boundary_conditions_parabolic, mesh, equations_parabolic,
-  #   dg.surface_integral, dg)
+  # Prolong solution to boundaries
+  @trixi_timeit timer() "prolong2boundaries" prolong2boundaries!(
+    cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg)
+
+  # Calculate boundary fluxes
+  @trixi_timeit timer() "boundary flux" calc_boundary_flux_gradients!(
+    cache_parabolic, t, boundary_conditions_parabolic, mesh, equations_parabolic,
+    dg.surface_integral, dg)
 
   # TODO: parabolic; mortars
 

From 56dd5775614d026881ec87a847a851ace121d926 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 13:59:09 -0500
Subject: [PATCH 16/46] remove todo

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index b15e63f6cfb..7121e32d90a 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -304,7 +304,6 @@ function prolong2interfaces!(cache_parabolic, flux_viscous,
   return nothing
 end
 
-# TODO: Deal with orientations for the BR1 divergence flux = {σ} ⋅ normal
 function calc_interface_flux!(surface_flux_values,
                               mesh::P4estMesh{2}, equations_parabolic,
                               dg::DG, cache_parabolic)

From 194a6953f2ed9cdb924e7ccdda6c3535a890342e Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 13:59:23 -0500
Subject: [PATCH 17/46] variable renaming

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 23 +++++++++++-----------
 1 file changed, 11 insertions(+), 12 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 7121e32d90a..6563a7f9475 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -317,7 +317,7 @@ function calc_interface_flux!(surface_flux_values,
     # Get element and side index information on the primary element
     primary_element = neighbor_ids[1, interface]
     primary_indices = node_indices[1, interface]
-    primary_direction = indices2direction(primary_indices)
+    primary_direction_index = indices2direction(primary_indices)
 
     # Create the local i,j indexing on the primary element used to pull normal direction information
     i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
@@ -329,7 +329,7 @@ function calc_interface_flux!(surface_flux_values,
     # Get element and side index information on the secondary element
     secondary_element = neighbor_ids[2, interface]
     secondary_indices = node_indices[2, interface]
-    secondary_direction = indices2direction(secondary_indices)
+    secondary_direction_index = indices2direction(secondary_indices)
 
     # Initiate the secondary index to be used in the surface for loop.
     # This index on the primary side will always run forward but
@@ -342,19 +342,18 @@ function calc_interface_flux!(surface_flux_values,
       node_secondary_step = 1
     end
 
-    for i in eachnode(dg)
-            
-      flux_viscous_normal_ll, flux_viscous_normal_rr = 
-        get_surface_node_vars(cache_parabolic.interfaces.u, 
-                              equations_parabolic, dg, i, interface)
+    for node in eachnode(dg)
+
+      # We prolong the viscous flux dotted with the primary normal to the boundaries. 
+      # Assuming a BR-1 type of flux, we 
+      viscous_flux_normal_ll, viscous_flux_normal_rr = 
+        get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, dg, node, interface)
 
-      # TODO: parabolic; only BR1 at the moment
-      flux = 0.5 * (flux_viscous_normal_ll + flux_viscous_normal_rr)
+      flux = 0.5 * (viscous_flux_normal_ll + viscous_flux_normal_rr)
 
-      # Copy flux to left and right element storage
       for v in eachvariable(equations_parabolic)
-        surface_flux_values[v, i, left_direction,  left_id]  = flux[v]
-        surface_flux_values[v, i, right_direction, right_id] = -flux[v]
+        surface_flux_values[v, node, primary_direction_index, primary_element] = flux[v]
+        surface_flux_values[v, node_secondary, secondary_direction_index, secondary_element] = -flux[v]
       end
 
       # Increment primary element indices to pull the normal direction

From ff23f39e90dfa2696331e1583878b8f64aafc73b Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 14:33:58 -0500
Subject: [PATCH 18/46] extending TreeMesh parabolic functions to P4estMesh

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

diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
index 0cb9c482196..07146c8d79e 100644
--- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
@@ -283,20 +283,20 @@ function get_unsigned_normal_vector_2d(direction)
 end
 
 function calc_boundary_flux_gradients!(cache, t, boundary_conditions_parabolic::BoundaryConditionPeriodic,
-                                      mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic,
+                                      mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations_parabolic::AbstractEquationsParabolic,
                                       surface_integral, dg::DG)
   return nothing
 end
 
 function calc_boundary_flux_divergence!(cache, t, boundary_conditions_parabolic::BoundaryConditionPeriodic,
-                                        mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic,
+                                        mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations_parabolic::AbstractEquationsParabolic,
                                         surface_integral, dg::DG)
   return nothing
 end
 
 function calc_boundary_flux_gradients!(cache, t, boundary_conditions_parabolic::NamedTuple,
-                                      mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic,
-                                      surface_integral, dg::DG)
+                                       mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic,
+                                       surface_integral, dg::DG)
   @unpack surface_flux_values = cache.elements
   @unpack n_boundaries_per_direction = cache.boundaries
 

From 2a16098887fa3e40840795e8f1eecdd6d1a6df2e Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 15:33:30 -0500
Subject: [PATCH 19/46] adding elixir

---
 .../elixir_advection_diffusion_periodic.jl    | 66 +++++++++++++++++++
 1 file changed, 66 insertions(+)
 create mode 100644 examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl

diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl
new file mode 100644
index 00000000000..6bb1f958260
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl
@@ -0,0 +1,66 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+diffusivity() = 5.0e-2
+advection_velocity = (1.0, 0.0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
+
+function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D)
+  return SVector(exp(-100 * (x[1]^2 + x[2]^2)))
+end
+initial_condition = initial_condition_sharp_gaussian
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg=3, initial_refinement_level=2,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 periodicity=true)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+# The AliveCallback prints short status information in regular intervals
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+time_int_tol = 1.0e-11
+sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
+            ode_default_options()..., callback=callbacks)
+
+# Print the timer summary
+summary_callback()

From de2513d0fe6dab6b2573bc8a3a72bc89dc049523 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 15:34:50 -0500
Subject: [PATCH 20/46] comments

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 6563a7f9475..c5ebf81d191 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -235,7 +235,7 @@ function prolong2interfaces!(cache_parabolic, flux_viscous,
                              mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic,
                              surface_integral, dg::DG, cache)
   (; interfaces) = cache_parabolic
-  (; contravariant_vectors) = cache_parabolic.elements
+  (; contravariant_vectors) = cache_parabolic.elements 
   index_range = eachnode(dg)
   flux_viscous_x, flux_viscous_y = flux_viscous
 
@@ -344,8 +344,8 @@ function calc_interface_flux!(surface_flux_values,
 
     for node in eachnode(dg)
 
-      # We prolong the viscous flux dotted with the primary normal to the boundaries. 
-      # Assuming a BR-1 type of flux, we 
+      # We prolong the viscous flux dotted with respect the outward normal on the 
+      # primary element. We assume a BR-1 type of flux.
       viscous_flux_normal_ll, viscous_flux_normal_rr = 
         get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, dg, node, interface)
 
@@ -366,3 +366,5 @@ function calc_interface_flux!(surface_flux_values,
 
   return nothing
 end
+
+# TODO: parabolic, finish implementing `prolong2boundaries!`, `calc_boundary_flux_gradients!` and `calc_boundary_flux_divergence!`
\ No newline at end of file

From f482e68ba72aab54a6cce3a85f4d39d6df05d712 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 15:38:55 -0500
Subject: [PATCH 21/46] add prolong2boundaries! (untested)

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 40 +++++++++++++++++++++-
 1 file changed, 39 insertions(+), 1 deletion(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index c5ebf81d191..81ad93aef90 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -367,4 +367,42 @@ function calc_interface_flux!(surface_flux_values,
   return nothing
 end
 
-# TODO: parabolic, finish implementing `prolong2boundaries!`, `calc_boundary_flux_gradients!` and `calc_boundary_flux_divergence!`
\ No newline at end of file
+# TODO: parabolic, finish implementing `calc_boundary_flux_gradients!` and `calc_boundary_flux_divergence!`
+function prolong2boundaries!(cache_parabolic, flux_viscous,
+                             mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic,
+                             surface_integral, dg::DG, cache)
+  (; boundaries) = cache_parabolic
+  (; contravariant_vectors) = cache_parabolic.elements
+  index_range = eachnode(dg)
+
+  flux_viscous_x, flux_viscous_y = flux_viscous
+
+  @threaded for boundary in eachboundary(dg, cache_parabolic)
+    # Copy solution data from the element using "delayed indexing" with
+    # a start value and a step size to get the correct face and orientation.
+    element       = boundaries.neighbor_ids[boundary]
+    node_indices  = boundaries.node_indices[boundary]
+
+    i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+    j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+    i_node = i_node_start
+    j_node = j_node_start
+    for i in eachnode(dg)
+      # this is the outward normal direction on the primary element
+      normal_direction = get_normal_direction(primary_direction, contravariant_vectors,
+                        i_node, j_node, primary_element)
+
+      for v in eachvariable(equations_parabolic)
+        flux_viscous = SVector(flux_viscous_x[v, i_primary, j_primary, primary_element], 
+            flux_viscous_y[v, i_primary, j_primary, primary_element])
+
+        boundaries.u[v, i, boundary] = dot(flux_viscous, normal_direction)
+      end
+      i_node += i_node_step
+      j_node += j_node_step
+    end
+  end
+
+  return nothing
+end

From 55278272144d358a8bc8d64463b8e0babd065ff1 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 15:44:18 -0500
Subject: [PATCH 22/46] update test

---
 .../elixir_advection_diffusion_periodic.jl    | 27 +++++++++++++++----
 test/test_parabolic_2d.jl                     |  9 +++++++
 2 files changed, 31 insertions(+), 5 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl
index 6bb1f958260..a13486ff3d6 100644
--- a/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl
+++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl
@@ -9,16 +9,33 @@ advection_velocity = (1.0, 0.0)
 equations = LinearScalarAdvectionEquation2D(advection_velocity)
 equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
 
-function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D)
-  return SVector(exp(-100 * (x[1]^2 + x[2]^2)))
+function x_trans_periodic(x, domain_length = SVector(2 * pi), center = SVector(0.0))
+    x_normalized = x .- center
+    x_shifted = x_normalized .% domain_length
+    x_offset = ((x_shifted .< -0.5 * domain_length) - (x_shifted .> 0.5 * domain_length)) .* domain_length
+    return center + x_shifted + x_offset
 end
-initial_condition = initial_condition_sharp_gaussian
+
+# Define initial condition (copied from "examples/tree_1d_dgsem/elixir_advection_diffusion.jl")
+function initial_condition_diffusive_convergence_test(x, t, equation::LinearScalarAdvectionEquation2D)
+    # Store translated coordinate for easy use of exact solution
+    # Assumes that advection_velocity[2] = 0 (effectively that we are solving a 1D equation)
+    x_trans = x_trans_periodic(x[1] - equation.advection_velocity[1] * t)
+    
+    nu = diffusivity()
+    c = 0.0
+    A = 1.0
+    omega = 1.0
+    scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t)
+    return SVector(scalar)
+end
+initial_condition = initial_condition_diffusive_convergence_test
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
 solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
 
-coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
-coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
+coordinates_min = (-pi, -pi) # minimum coordinates (min(x), min(y))
+coordinates_max = ( pi,  pi) # maximum coordinates (max(x), max(y))
 
 trees_per_dimension = (4, 4)
 mesh = P4estMesh(trees_per_dimension,
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 588f43e4543..e00bf728b76 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -184,6 +184,15 @@ isdir(outdir) && rm(outdir, recursive=true)
     )
   end
 
+  @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic.jl"),
+      trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5),
+      l2 = [0.002633614394267485], 
+      linf = [0.013444960515762605]
+    )
+  end
+
+
 end
 
 # Clean up afterwards: delete Trixi.jl output directory

From e5b91b594ad54934af1b3f6c2036e0b79aed27a2 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 15:48:58 -0500
Subject: [PATCH 23/46] initial commit

---
 .../elixir_navierstokes_convergence.jl        | 220 ++++++++++++++++++
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl    | 106 +++++++++
 2 files changed, 326 insertions(+)
 create mode 100644 examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
new file mode 100644
index 00000000000..0b49dc9197a
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
@@ -0,0 +1,220 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+prandtl_number() = 0.72
+mu() = 0.01
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(),
+                                                          gradient_variables=GradientVariablesPrimitive())
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
+               volume_integral=VolumeIntegralWeakForm())
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                periodicity=(true, true),
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg=3, initial_refinement_level=2,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 periodicity=true)
+
+# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D`
+#       since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`)
+#       and by the initial condition (which passes in `CompressibleEulerEquations2D`).
+# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000)
+function initial_condition_navier_stokes_convergence_test(x, t, equations)
+  # Amplitude and shift
+  A = 0.5
+  c = 2.0
+
+  # convenience values for trig. functions
+  pi_x = pi * x[1]
+  pi_y = pi * x[2]
+  pi_t = pi * t
+
+  rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+  v1  = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t)
+  v2  = v1
+  p   = rho^2
+
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations)
+  y = x[2]
+
+  # TODO: parabolic
+  # we currently need to hardcode these parameters until we fix the "combined equation" issue
+  # see also https://github.com/trixi-framework/Trixi.jl/pull/1160
+  inv_gamma_minus_one = inv(equations.gamma - 1)
+  Pr = prandtl_number()
+  mu_ = mu()
+
+  # Same settings as in `initial_condition`
+  # Amplitude and shift
+  A = 0.5
+  c = 2.0
+
+  # convenience values for trig. functions
+  pi_x = pi * x[1]
+  pi_y = pi * x[2]
+  pi_t = pi * t
+
+  # compute the manufactured solution and all necessary derivatives
+  rho    =  c  + A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+  rho_t  = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t)
+  rho_x  =  pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t)
+  rho_y  = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t)
+  rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+  rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+
+  v1    =       sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
+  v1_t  = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t)
+  v1_x  =  pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
+  v1_y  =       sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t)
+  v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
+  v1_xy =  pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t)
+  v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0)
+                         - A * A * log(y + 2.0) * exp(-A * (y - 1.0))
+                         - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t))
+  v2    = v1
+  v2_t  = v1_t
+  v2_x  = v1_x
+  v2_y  = v1_y
+  v2_xx = v1_xx
+  v2_xy = v1_xy
+  v2_yy = v1_yy
+
+  p    = rho * rho
+  p_t  = 2.0 * rho * rho_t
+  p_x  = 2.0 * rho * rho_x
+  p_y  = 2.0 * rho * rho_y
+  p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x
+  p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y
+
+  # Note this simplifies slightly because the ansatz assumes that v1 = v2
+  E   = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2)
+  E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t
+  E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x
+  E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y
+
+  # Some convenience constants
+  T_const = equations.gamma * inv_gamma_minus_one / Pr
+  inv_rho_cubed = 1.0 / (rho^3)
+
+  # compute the source terms
+  # density equation
+  du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y
+
+  # x-momentum equation
+  du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2
+                                        + 2.0   * rho  * v1 * v1_x
+                                        + rho_y * v1   * v2
+                                        + rho   * v1_y * v2
+                                        + rho   * v1   * v2_y
+    # stress tensor from x-direction
+                      - 4.0 / 3.0 * v1_xx * mu_
+                      + 2.0 / 3.0 * v2_xy * mu_
+                      - v1_yy             * mu_
+                      - v2_xy             * mu_ )
+  # y-momentum equation
+  du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1    * v2
+                                        + rho   * v1_x  * v2
+                                        + rho   * v1    * v2_x
+                                        +         rho_y * v2^2
+                                        + 2.0   * rho   * v2 * v2_y
+    # stress tensor from y-direction
+                      - v1_xy             * mu_
+                      - v2_xx             * mu_
+                      - 4.0 / 3.0 * v2_yy * mu_
+                      + 2.0 / 3.0 * v1_xy * mu_ )
+  # total energy equation
+  du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x)
+              + v2_y * (E + p) + v2 * (E_y + p_y)
+    # stress tensor and temperature gradient terms from x-direction
+                                - 4.0 / 3.0 * v1_xx * v1   * mu_
+                                + 2.0 / 3.0 * v2_xy * v1   * mu_
+                                - 4.0 / 3.0 * v1_x  * v1_x * mu_
+                                + 2.0 / 3.0 * v2_y  * v1_x * mu_
+                                - v1_xy     * v2           * mu_
+                                - v2_xx     * v2           * mu_
+                                - v1_y      * v2_x         * mu_
+                                - v2_x      * v2_x         * mu_
+         - T_const * inv_rho_cubed * (        p_xx * rho   * rho
+                                      - 2.0 * p_x  * rho   * rho_x
+                                      + 2.0 * p    * rho_x * rho_x
+                                      -       p    * rho   * rho_xx ) * mu_
+    # stress tensor and temperature gradient terms from y-direction
+                                - v1_yy     * v1           * mu_
+                                - v2_xy     * v1           * mu_
+                                - v1_y      * v1_y         * mu_
+                                - v2_x      * v1_y         * mu_
+                                - 4.0 / 3.0 * v2_yy * v2   * mu_
+                                + 2.0 / 3.0 * v1_xy * v2   * mu_
+                                - 4.0 / 3.0 * v2_y  * v2_y * mu_
+                                + 2.0 / 3.0 * v1_x  * v2_y * mu_
+         - T_const * inv_rho_cubed * (        p_yy * rho   * rho
+                                      - 2.0 * p_y  * rho   * rho_y
+                                      + 2.0 * p    * rho_y * rho_y
+                                      -       p    * rho   * rho_yy ) * mu_ )
+
+  return SVector(du1, du2, du3, du4)
+end
+
+initial_condition = initial_condition_navier_stokes_convergence_test
+
+# BC types
+# TODO: parabolic. This is wrong! We got lucky because the velocity is actually zero at the boundaries. 
+velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3])
+heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
+boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom)
+
+# define inviscid boundary conditions
+boundary_conditions = (; x_neg = boundary_condition_periodic,
+                         x_pos = boundary_condition_periodic,
+                         y_neg = boundary_condition_slip_wall,
+                         y_pos = boundary_condition_slip_wall)
+
+# define viscous boundary conditions
+boundary_conditions_parabolic = (; x_neg = boundary_condition_periodic,
+                                   x_pos = boundary_condition_periodic,
+                                   y_neg = boundary_condition_top_bottom,
+                                   y_pos = boundary_condition_top_bottom)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver;
+                                             boundary_conditions=(boundary_conditions, boundary_conditions_parabolic),
+                                             source_terms=source_terms_navier_stokes_convergence_test)
+
+# ###############################################################################
+# # ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-8
+sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5,
+            ode_default_options()..., callback=callbacks)
+summary_callback() # print the timer summary
+
diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 81ad93aef90..4764a228485 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -406,3 +406,109 @@ function prolong2boundaries!(cache_parabolic, flux_viscous,
 
   return nothing
 end
+
+# function calc_boundary_flux_gradients!(cache, t, boundary_conditions_parabolic::NamedTuple,
+#                                        mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic,
+#                                        surface_integral, dg::DG)
+#   (; boundaries) = cache
+#   (; node_coordinates, surface_flux_values) = cache.elements
+#   (; boundary_condition_types, boundary_indices = boundary_conditions) = boundary_conditions_parabolic
+#   index_range = eachnode(dg)
+
+#   @threaded for local_index in eachindex(boundary_indexing)
+#     # Use the local index to get the global boundary index from the pre-sorted list
+#     boundary = boundary_indexing[local_index]
+
+#     # Get information on the adjacent element, compute the surface fluxes,
+#     # and store them
+#     element         = boundaries.neighbor_ids[boundary]
+#     node_indices    = boundaries.node_indices[boundary]
+#     direction_index = indices2direction(node_indices)
+
+#     i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+#     j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+#     i_node = i_node_start
+#     j_node = j_node_start
+#     for node in eachnode(dg)
+#       # Extract solution data from boundary container
+#       u_inner = get_node_vars(boundaries.u, equations_parabolic, dg, node, boundary)
+
+#       # Outward-pointing normal direction (not normalized)
+#       normal_direction = get_normal_direction(direction_index, contravariant_vectors,
+#                                               i_node, j_node, element)
+
+#       # TODO: revisit if we want more general boundary treatments.
+#       # This assumes the gradient numerical flux at the boundary is the gradient variable,
+#       # which is consistent with BR1, LDG.
+#       flux_inner = u_inner
+
+#       # Coordinates at boundary node
+#       x = get_node_coords(node_coordinates, equations_parabolic, dg, i_node, j_node, element)
+
+#       flux_ = boundary_condition(flux_inner, u_inner, normal_direction, 
+#                                  x, t, Gradient(), equations_parabolic)
+
+#       # Copy flux to element storage in the correct orientation
+#       for v in eachvariable(equations_parabolic)
+#         surface_flux_values[v, node_index, direction_index, element_index] = flux_[v]
+#       end
+
+#       i_node += i_node_step
+#       j_node += j_node_step
+#     end
+#   end
+# end
+
+# function calc_boundary_flux_divergence!(cache, t, boundary_conditions_parabolic::NamedTuple,
+#                                         mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic,
+#                                         surface_integral, dg::DG)
+#   (; boundaries) = cache
+#   (; node_coordinates, surface_flux_values) = cache.elements
+#   index_range = eachnode(dg)
+
+#   @threaded for local_index in eachindex(boundary_indexing)
+#     # Use the local index to get the global boundary index from the pre-sorted list
+#     boundary = boundary_indexing[local_index]
+
+#     # Get information on the adjacent element, compute the surface fluxes,
+#     # and store them
+#     element         = boundaries.neighbor_ids[boundary]
+#     node_indices    = boundaries.node_indices[boundary]
+#     direction_index = indices2direction(node_indices)
+
+#     i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+#     j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+#     i_node = i_node_start
+#     j_node = j_node_start
+#     for node in eachnode(dg)
+#       # Extract solution data from boundary container
+#       flux_inner = get_node_vars(boundaries.u, equations_parabolic, dg, node, boundary)
+
+#       # Outward-pointing normal direction (not normalized)
+#       normal_direction = get_normal_direction(direction_index, contravariant_vectors,
+#                                               i_node, j_node, element)
+
+#       # Coordinates at boundary node
+#       x = get_node_coords(node_coordinates, equations_parabolic, dg, i_node, j_node, element)
+
+#       # TODO: add a field in `cache.boundaries` for gradient information.
+#       # Here, we pass in `u_inner = nothing` since we overwrite cache.boundaries.u with gradient information.
+#       # This currently works with Dirichlet/Neuman boundary conditions for LaplaceDiffusion2D and
+#       # NoSlipWall/Adiabatic boundary conditions for CompressibleNavierStokesDiffusion2D as of 2022-6-27.
+#       # It will not work with implementations which utilize `u_inner` to impose boundary conditions.
+#       flux_ = boundary_condition(flux_inner, nothing, normal_direction,
+#                                  x, t, Divergence(), equations_parabolic)
+
+
+#       # Copy flux to element storage in the correct orientation
+#       for v in eachvariable(equations_parabolic)
+#         surface_flux_values[v, node_index, direction_index, element_index] = flux_[v]
+#       end
+
+#       i_node += i_node_step
+#       j_node += j_node_step
+#     end
+#   end                                        
+# end
\ No newline at end of file

From 4d5776f1f27c27bc71138360025e68085353b7e4 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 26 May 2023 16:36:48 -0500
Subject: [PATCH 24/46] fix CI

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

diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index e00bf728b76..09d84361dd9 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -187,8 +187,8 @@ isdir(outdir) && rm(outdir, recursive=true)
   @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin
     @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic.jl"),
       trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5),
-      l2 = [0.002633614394267485], 
-      linf = [0.013444960515762605]
+      l2 = [0.0023754695605828443], 
+      linf = [0.008154128363741964]
     )
   end
 

From 2d5804983f1ff6be9d95c3ee85a6b66ec4be6547 Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Sat, 27 May 2023 05:29:09 -0500
Subject: [PATCH 25/46] Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 81ad93aef90..a5cc9a5615f 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -1,4 +1,4 @@
-# This method is called when a SemidiscretizationHyperbolic is constructed.
+# This method is called when a SemidiscretizationHyperbolicParabolic is constructed.
 # It constructs the basic `cache` used throughout the simulation to compute
 # the RHS etc.
 function create_cache_parabolic(mesh::P4estMesh, equations_hyperbolic::AbstractEquations,

From 55368dcb9a67042a89de10c925736b6ad1b4eacd Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Sat, 27 May 2023 05:29:18 -0500
Subject: [PATCH 26/46] Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index a5cc9a5615f..e274d6b6ab6 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -117,7 +117,7 @@ function calc_gradient!(gradients, u_transformed, t,
 
       for node in eachnode(dg)
         u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u,
-                                          equations_parabolic, dg, node, interface)                                              
+                                           equations_parabolic, dg, node, interface)                                              
         flux = 0.5 * (u_ll + u_rr)
 
         for v in eachvariable(equations_parabolic)

From f3bb821f001bc9e6bd7f5cb40220c5bfc81a3302 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Sat, 27 May 2023 06:18:11 -0500
Subject: [PATCH 27/46] add "no mortars" check

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index e274d6b6ab6..f9d1372b35a 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -144,6 +144,7 @@ function calc_gradient!(gradients, u_transformed, t,
     dg.surface_integral, dg)
 
   # TODO: parabolic; mortars
+  @assert nmortars(dg, cache) == 0
 
   # Calculate surface integrals
   @trixi_timeit timer() "surface integral" begin

From 1315428788e510dc0772563965337dd35eb4e6c6 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Sat, 27 May 2023 06:49:03 -0500
Subject: [PATCH 28/46] add curved elixir

---
 ...xir_advection_diffusion_periodic_curved.jl | 86 +++++++++++++++++++
 1 file changed, 86 insertions(+)
 create mode 100644 examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl

diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl
new file mode 100644
index 00000000000..30039c2f69e
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl
@@ -0,0 +1,86 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+diffusivity() = 5.0e-2
+advection_velocity = (1.0, 0.0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
+
+function x_trans_periodic(x, domain_length = SVector(2 * pi), center = SVector(0.0))
+    x_normalized = x .- center
+    x_shifted = x_normalized .% domain_length
+    x_offset = ((x_shifted .< -0.5 * domain_length) - (x_shifted .> 0.5 * domain_length)) .* domain_length
+    return center + x_shifted + x_offset
+end
+
+# Define initial condition (copied from "examples/tree_1d_dgsem/elixir_advection_diffusion.jl")
+function initial_condition_diffusive_convergence_test(x, t, equation::LinearScalarAdvectionEquation2D)
+    # Store translated coordinate for easy use of exact solution
+    # Assumes that advection_velocity[2] = 0 (effectively that we are solving a 1D equation)
+    x_trans = x_trans_periodic(x[1] - equation.advection_velocity[1] * t)
+    
+    nu = diffusivity()
+    c = 0.0
+    A = 1.0
+    omega = 1.0
+    scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t)
+    return SVector(scalar)
+end
+initial_condition = initial_condition_diffusive_convergence_test
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+function mapping(xi, eta)
+  x = xi  + 0.1 * sin(pi * xi) * sin(pi * eta)
+  y = eta + 0.1 * sin(pi * xi) * sin(pi * eta)
+  return pi * SVector(x, y)
+end
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg=3, initial_refinement_level=2,                 
+                 mapping=mapping,
+                 periodicity=true)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+# The AliveCallback prints short status information in regular intervals
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+time_int_tol = 1.0e-11
+sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
+            ode_default_options()..., callback=callbacks)
+
+# Print the timer summary
+summary_callback()

From 9ced755cb6011a6c50f9a0a7e644dce0c13e9026 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Wed, 31 May 2023 12:01:04 -0500
Subject: [PATCH 29/46] fix gradient bug

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 26 +++++++++++++++++++++-
 1 file changed, 25 insertions(+), 1 deletion(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index f9d1372b35a..5d9d962f229 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -152,7 +152,6 @@ function calc_gradient!(gradients, u_transformed, t,
     (; surface_flux_values) = cache_parabolic.elements
     (; contravariant_vectors) = cache.elements
 
-    # Note that all fluxes have been computed with outward-pointing normal vectors.
     # Access the factors only once before beginning the loop to increase performance.
     # We also use explicit assignments instead of `+=` to let `@muladd` turn these
     # into FMAs (see comment at the top of the file).
@@ -161,6 +160,9 @@ function calc_gradient!(gradients, u_transformed, t,
     @threaded for element in eachelement(dg, cache)
       for l in eachnode(dg)
         for v in eachvariable(equations_parabolic)
+
+          # Compute x-component of gradients
+
           # surface at -x
           normal_direction_x, _ = get_normal_direction(1, contravariant_vectors, 1, l, element)
           gradients_x[v, 1,          l, element] = (
@@ -171,6 +173,28 @@ function calc_gradient!(gradients, u_transformed, t,
           gradients_x[v, nnodes(dg), l, element] = (
             gradients_x[v, nnodes(dg), l, element] + surface_flux_values[v, l, 2, element] * factor_2 * normal_direction_x)
 
+          # surface at -y
+          normal_direction_x, _ = get_normal_direction(3, contravariant_vectors, l, 1, element)
+          gradients_x[v, l, 1,          element] = (
+            gradients_x[v, l, 1,          element] + surface_flux_values[v, l, 3, element] * factor_1 * normal_direction_x)
+
+          # surface at +y
+          normal_direction_x, _ = get_normal_direction(4, contravariant_vectors, l, nnodes(dg), element)
+          gradients_x[v, l, nnodes(dg), element] = (
+            gradients_x[v, l, nnodes(dg), element] + surface_flux_values[v, l, 4, element] * factor_2 * normal_direction_x)
+
+          # Compute y-component of gradients
+
+          # surface at -x
+          _, normal_direction_y = get_normal_direction(1, contravariant_vectors, 1, l, element)
+          gradients_y[v, 1,          l, element] = (
+            gradients_y[v, 1,          l, element] + surface_flux_values[v, l, 1, element] * factor_1 * normal_direction_y)
+
+          # surface at +x
+          _, normal_direction_y = get_normal_direction(2, contravariant_vectors, nnodes(dg), l, element)
+          gradients_y[v, nnodes(dg), l, element] = (
+            gradients_y[v, nnodes(dg), l, element] + surface_flux_values[v, l, 2, element] * factor_2 * normal_direction_y)
+
           # surface at -y
           _, normal_direction_y = get_normal_direction(3, contravariant_vectors, l, 1, element)
           gradients_y[v, l, 1,          element] = (

From 1b3c1d24738ba7ab3c6a271f8a880010d1fc550c Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Wed, 31 May 2023 12:03:16 -0500
Subject: [PATCH 30/46] add curved test

---
 test/test_parabolic_2d.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 09d84361dd9..b0ac63d4ce9 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -192,6 +192,13 @@ isdir(outdir) && rm(outdir, recursive=true)
     )
   end
 
+  @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic_curved.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic_curved.jl"),
+      trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5),
+      l2 = [0.012380458938507371], 
+      linf = [0.10860506906472567]
+    )
+  end
 
 end
 

From 1cce0879e6d4b377ff4ab752692a682003368452 Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Thu, 1 Jun 2023 09:08:39 -0500
Subject: [PATCH 31/46] Apply suggestions from code review

Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 .../p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl     | 2 +-
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl                    | 4 +---
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl
index a13486ff3d6..1cd075e84ea 100644
--- a/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl
+++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic.jl
@@ -9,7 +9,7 @@ advection_velocity = (1.0, 0.0)
 equations = LinearScalarAdvectionEquation2D(advection_velocity)
 equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
 
-function x_trans_periodic(x, domain_length = SVector(2 * pi), center = SVector(0.0))
+function x_trans_periodic(x, domain_length=SVector(2 * pi), center=SVector(0.0))
     x_normalized = x .- center
     x_shifted = x_normalized .% domain_length
     x_offset = ((x_shifted .< -0.5 * domain_length) - (x_shifted .> 0.5 * domain_length)) .* domain_length
diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 5d9d962f229..8723a22547f 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -165,8 +165,7 @@ function calc_gradient!(gradients, u_transformed, t,
 
           # surface at -x
           normal_direction_x, _ = get_normal_direction(1, contravariant_vectors, 1, l, element)
-          gradients_x[v, 1,          l, element] = (
-            gradients_x[v, 1,          l, element] + surface_flux_values[v, l, 1, element] * factor_1 * normal_direction_x)
+          gradients_x[v, 1,          l, element] += surface_flux_values[v, l, 1, element] * factor_1 * normal_direction_x
 
           # surface at +x
           normal_direction_x, _ = get_normal_direction(2, contravariant_vectors, nnodes(dg), l, element)
@@ -368,7 +367,6 @@ function calc_interface_flux!(surface_flux_values,
     end
 
     for node in eachnode(dg)
-
       # We prolong the viscous flux dotted with respect the outward normal on the 
       # primary element. We assume a BR-1 type of flux.
       viscous_flux_normal_ll, viscous_flux_normal_rr = 

From 3886a14d9df7e025960e1a2c6dae183979c1f6e6 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Thu, 1 Jun 2023 09:49:06 -0500
Subject: [PATCH 32/46] add comment on mapping

---
 .../elixir_advection_diffusion_periodic_curved.jl               | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl
index 30039c2f69e..e63f7947952 100644
--- a/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl
+++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl
@@ -34,6 +34,8 @@ initial_condition = initial_condition_diffusive_convergence_test
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
 solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
 
+# This maps the domain [-1, 1]^2 to [-pi, pi]^2 while also 
+# introducing a curved warping to interior nodes. 
 function mapping(xi, eta)
   x = xi  + 0.1 * sin(pi * xi) * sin(pi * eta)
   y = eta + 0.1 * sin(pi * xi) * sin(pi * eta)

From 8976bc400069410f02944070ed33cbd8b3266d32 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Thu, 1 Jun 2023 09:49:17 -0500
Subject: [PATCH 33/46] reuse P4estMesh{2} code

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 83 ++++++++--------------
 1 file changed, 29 insertions(+), 54 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 8723a22547f..f62e9d17e61 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -79,60 +79,12 @@ function calc_gradient!(gradients, u_transformed, t,
   @trixi_timeit timer() "prolong2interfaces" prolong2interfaces!(
     cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg)
 
-  # Calculate interface fluxes
-  @trixi_timeit timer() "interface flux" begin
-    (; surface_flux_values) = cache_parabolic.elements
-    (; neighbor_ids, node_indices) = cache.interfaces
-    index_range = eachnode(dg)
-    index_end = last(index_range)
-
-    @threaded for interface in eachinterface(dg, cache)
-      # Get element and side index information on the primary element
-      primary_element = neighbor_ids[1, interface]
-      primary_indices = node_indices[1, interface]
-      primary_direction = indices2direction(primary_indices)
-
-      # Create the local i,j indexing on the primary element used to pull normal direction information
-      i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
-      j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)
-
-      i_primary = i_primary_start
-      j_primary = j_primary_start
-
-      # Get element and side index information on the secondary element
-      secondary_element = neighbor_ids[2, interface]
-      secondary_indices = node_indices[2, interface]
-      secondary_direction = indices2direction(secondary_indices)
-
-      # Initiate the secondary index to be used in the surface for loop.
-      # This index on the primary side will always run forward but
-      # the secondary index might need to run backwards for flipped sides.
-      if :i_backward in secondary_indices
-        node_secondary = index_end
-        node_secondary_step = -1
-      else
-        node_secondary = 1
-        node_secondary_step = 1
-      end
-
-      for node in eachnode(dg)
-        u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u,
-                                           equations_parabolic, dg, node, interface)                                              
-        flux = 0.5 * (u_ll + u_rr)
-
-        for v in eachvariable(equations_parabolic)
-          surface_flux_values[v, node, primary_direction, primary_element] = flux[v]
-          surface_flux_values[v, node_secondary, secondary_direction, secondary_element] = flux[v]
-        end
-    
-        # Increment primary element indices to pull the normal direction
-        i_primary += i_primary_step
-        j_primary += j_primary_step
-        # Increment the surface node index along the secondary element
-        node_secondary += node_secondary_step
-      end
-    end
-  end
+  # Calculate interface fluxes for the gradient. This reuses P4est `calc_interface_flux!` along with a 
+  # specialization for AbstractEquationsParabolic.
+  @trixi_timeit timer() "interface flux" calc_interface_flux!(cache_parabolic.elements.surface_flux_values, 
+                                                              mesh, False(), # False() = no nonconservative terms
+                                                              equations_parabolic, dg.surface_integral, dg, 
+                                                              cache_parabolic)
 
   # Prolong solution to boundaries
   @trixi_timeit timer() "prolong2boundaries" prolong2boundaries!(
@@ -217,6 +169,29 @@ function calc_gradient!(gradients, u_transformed, t,
   return nothing
 end
 
+# This version is used for parabolic gradient computations
+@inline function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2},
+                                      nonconservative_terms::False, 
+                                      equations::AbstractEquationsParabolic,
+                                      surface_integral, dg::DG, cache,
+                                      interface_index, normal_direction,
+                                      primary_node_index, primary_direction_index, primary_element_index,
+                                      secondary_node_index, secondary_direction_index, secondary_element_index)
+  @unpack u = cache.interfaces
+  @unpack surface_flux = surface_integral
+
+  u_ll, u_rr = get_surface_node_vars(u, equations, dg, primary_node_index, interface_index)
+
+  flux_ = 0.5 * (u_ll + u_rr) # we assume that the gradient computations utilize a central flux
+
+  # Note that we don't flip the sign on the secondondary flux. This is because for parabolic terms,
+  # the normals are not embedded in `flux_` for the parabolic gradient computations. 
+  for v in eachvariable(equations)
+    surface_flux_values[v, primary_node_index, primary_direction_index, primary_element_index] = flux_[v]
+    surface_flux_values[v, secondary_node_index, secondary_direction_index, secondary_element_index] = flux_[v]
+  end
+end
+
 # This is the version used when calculating the divergence of the viscous fluxes
 function calc_volume_integral!(du, flux_viscous,
                                mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic,

From 0a3f1fb7bfa0f75680d42ea861c2fe30b9df5880 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Thu, 1 Jun 2023 10:07:37 -0500
Subject: [PATCH 34/46] fix += for muladd

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index f62e9d17e61..e9c3ee31853 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -117,7 +117,8 @@ function calc_gradient!(gradients, u_transformed, t,
 
           # surface at -x
           normal_direction_x, _ = get_normal_direction(1, contravariant_vectors, 1, l, element)
-          gradients_x[v, 1,          l, element] += surface_flux_values[v, l, 1, element] * factor_1 * normal_direction_x
+          gradients_x[v, 1,          l, element] = (
+            gradients_x[v, 1,          l, element] + surface_flux_values[v, l, 1, element] * factor_1 * normal_direction_x)
 
           # surface at +x
           normal_direction_x, _ = get_normal_direction(2, contravariant_vectors, nnodes(dg), l, element)

From 193ec00c535aca443063aa3576adc54ef3bec15f Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Thu, 1 Jun 2023 10:09:36 -0500
Subject: [PATCH 35/46] Update
 examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl

Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
---
 .../elixir_advection_diffusion_periodic_curved.jl               | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl
index e63f7947952..b438fb8a29c 100644
--- a/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl
+++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl
@@ -9,7 +9,7 @@ advection_velocity = (1.0, 0.0)
 equations = LinearScalarAdvectionEquation2D(advection_velocity)
 equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
 
-function x_trans_periodic(x, domain_length = SVector(2 * pi), center = SVector(0.0))
+function x_trans_periodic(x, domain_length=SVector(2 * pi), center=SVector(0.0))
     x_normalized = x .- center
     x_shifted = x_normalized .% domain_length
     x_offset = ((x_shifted .< -0.5 * domain_length) - (x_shifted .> 0.5 * domain_length)) .* domain_length

From 895658722777d15a2c8a26e360b6804f5ff882fc Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Thu, 1 Jun 2023 10:57:36 -0500
Subject: [PATCH 36/46] comment

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index e9c3ee31853..7ddb83f97db 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -65,6 +65,10 @@ function calc_gradient!(gradients, u_transformed, t,
         gradients_reference_1 = get_node_vars(gradients_x, equations_parabolic, dg, i, j, element)
         gradients_reference_2 = get_node_vars(gradients_y, equations_parabolic, dg, i, j, element)
 
+        # note that the contravariant vectors are transposed compared with computations of flux
+        # divergences in `calc_volume_integral!`. See 
+        # https://github.com/trixi-framework/Trixi.jl/pull/1490#discussion_r1213345190 
+        # for a more detailed discussion. 
         gradient_x_node = Ja11 * gradients_reference_1 + Ja21 * gradients_reference_2
         gradient_y_node = Ja12 * gradients_reference_1 + Ja22 * gradients_reference_2
 

From 2fe0aa10fabe8de8dd0a4daa5263ddb30ba12628 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 2 Jun 2023 09:37:22 -0500
Subject: [PATCH 37/46] comments + remove cruft

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

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 7b3213f2e1b..25922fbe53e 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -308,12 +308,12 @@ function prolong2interfaces!(cache_parabolic, flux_viscous,
   return nothing
 end
 
+# This version is used for divergence flux computations 
 function calc_interface_flux!(surface_flux_values,
                               mesh::P4estMesh{2}, equations_parabolic,
                               dg::DG, cache_parabolic)
 
   (; neighbor_ids, node_indices) = cache_parabolic.interfaces
-  (; contravariant_vectors) = cache_parabolic.elements
   index_range = eachnode(dg)
   index_end = last(index_range)
 

From 1ffda402eea8394ed65443205fd0ebc8724f30d2 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Mon, 5 Jun 2023 15:24:44 -0500
Subject: [PATCH 38/46] add BCs for parabolic P43st

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 230 +++++++++++----------
 1 file changed, 121 insertions(+), 109 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index 25922fbe53e..671205daf0f 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -385,6 +385,7 @@ function prolong2boundaries!(cache_parabolic, flux_viscous,
     # a start value and a step size to get the correct face and orientation.
     element       = boundaries.neighbor_ids[boundary]
     node_indices  = boundaries.node_indices[boundary]
+    direction     = indices2direction(node_indices)
 
     i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
     j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
@@ -393,12 +394,12 @@ function prolong2boundaries!(cache_parabolic, flux_viscous,
     j_node = j_node_start
     for i in eachnode(dg)
       # this is the outward normal direction on the primary element
-      normal_direction = get_normal_direction(primary_direction, contravariant_vectors,
-                        i_node, j_node, primary_element)
+      normal_direction = get_normal_direction(direction, contravariant_vectors,
+                                              i_node, j_node, element)
 
       for v in eachvariable(equations_parabolic)
-        flux_viscous = SVector(flux_viscous_x[v, i_primary, j_primary, primary_element], 
-            flux_viscous_y[v, i_primary, j_primary, primary_element])
+        flux_viscous = SVector(flux_viscous_x[v, i_node, j_node, element], 
+                               flux_viscous_y[v, i_node, j_node, element])
 
         boundaries.u[v, i, boundary] = dot(flux_viscous, normal_direction)
       end
@@ -410,108 +411,119 @@ function prolong2boundaries!(cache_parabolic, flux_viscous,
   return nothing
 end
 
-# function calc_boundary_flux_gradients!(cache, t, boundary_conditions_parabolic::NamedTuple,
-#                                        mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic,
-#                                        surface_integral, dg::DG)
-#   (; boundaries) = cache
-#   (; node_coordinates, surface_flux_values) = cache.elements
-#   (; boundary_condition_types, boundary_indices = boundary_conditions) = boundary_conditions_parabolic
-#   index_range = eachnode(dg)
-
-#   @threaded for local_index in eachindex(boundary_indexing)
-#     # Use the local index to get the global boundary index from the pre-sorted list
-#     boundary = boundary_indexing[local_index]
-
-#     # Get information on the adjacent element, compute the surface fluxes,
-#     # and store them
-#     element         = boundaries.neighbor_ids[boundary]
-#     node_indices    = boundaries.node_indices[boundary]
-#     direction_index = indices2direction(node_indices)
-
-#     i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
-#     j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
-
-#     i_node = i_node_start
-#     j_node = j_node_start
-#     for node in eachnode(dg)
-#       # Extract solution data from boundary container
-#       u_inner = get_node_vars(boundaries.u, equations_parabolic, dg, node, boundary)
-
-#       # Outward-pointing normal direction (not normalized)
-#       normal_direction = get_normal_direction(direction_index, contravariant_vectors,
-#                                               i_node, j_node, element)
-
-#       # TODO: revisit if we want more general boundary treatments.
-#       # This assumes the gradient numerical flux at the boundary is the gradient variable,
-#       # which is consistent with BR1, LDG.
-#       flux_inner = u_inner
-
-#       # Coordinates at boundary node
-#       x = get_node_coords(node_coordinates, equations_parabolic, dg, i_node, j_node, element)
-
-#       flux_ = boundary_condition(flux_inner, u_inner, normal_direction, 
-#                                  x, t, Gradient(), equations_parabolic)
-
-#       # Copy flux to element storage in the correct orientation
-#       for v in eachvariable(equations_parabolic)
-#         surface_flux_values[v, node_index, direction_index, element_index] = flux_[v]
-#       end
-
-#       i_node += i_node_step
-#       j_node += j_node_step
-#     end
-#   end
-# end
-
-# function calc_boundary_flux_divergence!(cache, t, boundary_conditions_parabolic::NamedTuple,
-#                                         mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic,
-#                                         surface_integral, dg::DG)
-#   (; boundaries) = cache
-#   (; node_coordinates, surface_flux_values) = cache.elements
-#   index_range = eachnode(dg)
-
-#   @threaded for local_index in eachindex(boundary_indexing)
-#     # Use the local index to get the global boundary index from the pre-sorted list
-#     boundary = boundary_indexing[local_index]
-
-#     # Get information on the adjacent element, compute the surface fluxes,
-#     # and store them
-#     element         = boundaries.neighbor_ids[boundary]
-#     node_indices    = boundaries.node_indices[boundary]
-#     direction_index = indices2direction(node_indices)
-
-#     i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
-#     j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
-
-#     i_node = i_node_start
-#     j_node = j_node_start
-#     for node in eachnode(dg)
-#       # Extract solution data from boundary container
-#       flux_inner = get_node_vars(boundaries.u, equations_parabolic, dg, node, boundary)
-
-#       # Outward-pointing normal direction (not normalized)
-#       normal_direction = get_normal_direction(direction_index, contravariant_vectors,
-#                                               i_node, j_node, element)
-
-#       # Coordinates at boundary node
-#       x = get_node_coords(node_coordinates, equations_parabolic, dg, i_node, j_node, element)
-
-#       # TODO: add a field in `cache.boundaries` for gradient information.
-#       # Here, we pass in `u_inner = nothing` since we overwrite cache.boundaries.u with gradient information.
-#       # This currently works with Dirichlet/Neuman boundary conditions for LaplaceDiffusion2D and
-#       # NoSlipWall/Adiabatic boundary conditions for CompressibleNavierStokesDiffusion2D as of 2022-6-27.
-#       # It will not work with implementations which utilize `u_inner` to impose boundary conditions.
-#       flux_ = boundary_condition(flux_inner, nothing, normal_direction,
-#                                  x, t, Divergence(), equations_parabolic)
-
-
-#       # Copy flux to element storage in the correct orientation
-#       for v in eachvariable(equations_parabolic)
-#         surface_flux_values[v, node_index, direction_index, element_index] = flux_[v]
-#       end
-
-#       i_node += i_node_step
-#       j_node += j_node_step
-#     end
-#   end                                        
-# end
\ No newline at end of file
+function calc_boundary_flux_gradients!(cache, t, 
+                                       boundary_condition::Union{BoundaryConditionPeriodic, BoundaryConditionDoNothing}, 
+                                       mesh::P4estMesh, equations, surface_integral, dg::DG)
+  @assert isempty(eachboundary(dg, cache))
+end
+
+# Function barrier for type stability
+function calc_boundary_flux_gradients!(cache, t, boundary_conditions, mesh::P4estMesh,
+                                       equations, surface_integral, dg::DG)
+  (; boundary_condition_types, boundary_indices) = boundary_conditions
+
+  calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices,
+                              Gradient(), mesh, equations, surface_integral, dg)
+  return nothing
+end
+
+function calc_boundary_flux_divergence!(cache, t, boundary_conditions, mesh::P4estMesh,
+                                        equations, surface_integral, dg::DG)
+  (; boundary_condition_types, boundary_indices) = boundary_conditions
+
+  calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices,
+                              Divergence(), mesh, equations, surface_integral, dg)
+  return nothing
+end
+
+# Iterate over tuples of boundary condition types and associated indices
+# in a type-stable way using "lispy tuple programming".
+function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N, Any},
+                                     BC_indices::NTuple{N, Vector{Int}},
+                                     operator_type, 
+                                     mesh::P4estMesh,
+                                     equations, surface_integral, dg::DG) where {N}
+  # Extract the boundary condition type and index vector
+  boundary_condition = first(BCs)
+  boundary_condition_indices = first(BC_indices)
+  # Extract the remaining types and indices to be processed later
+  remaining_boundary_conditions = Base.tail(BCs)
+  remaining_boundary_condition_indices = Base.tail(BC_indices)
+
+  # process the first boundary condition type
+  calc_boundary_flux!(cache, t, boundary_condition, boundary_condition_indices,
+                      operator_type, mesh, equations, surface_integral, dg)
+
+  # recursively call this method with the unprocessed boundary types
+  calc_boundary_flux_by_type!(cache, t, remaining_boundary_conditions,
+                              remaining_boundary_condition_indices,
+                              operator_type, 
+                              mesh, equations, surface_integral, dg)
+
+  return nothing
+end
+
+# terminate the type-stable iteration over tuples
+function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{},
+                                     operator_type, mesh::P4estMesh, equations, 
+                                     surface_integral, dg::DG)
+  nothing
+end
+
+
+function calc_boundary_flux!(cache, t, 
+                             boundary_condition_parabolic, # works with Dict types
+                             boundary_condition_indices, 
+                             operator_type,  mesh::P4estMesh{2}, 
+                             equations_parabolic::AbstractEquationsParabolic,
+                             surface_integral, dg::DG)
+  (; boundaries) = cache
+  (; node_coordinates, surface_flux_values) = cache.elements
+  (; contravariant_vectors) = cache.elements
+  index_range = eachnode(dg)
+
+  @threaded for local_index in eachindex(boundary_condition_indices)
+    # Use the local index to get the global boundary index from the pre-sorted list
+    boundary_index = boundary_condition_indices[local_index]
+
+    # Get information on the adjacent element, compute the surface fluxes,
+    # and store them
+    element         = boundaries.neighbor_ids[boundary_index]
+    node_indices    = boundaries.node_indices[boundary_index]
+    direction_index = indices2direction(node_indices)
+
+    i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+    j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+    i_node = i_node_start
+    j_node = j_node_start
+    for node_index in eachnode(dg)
+      # Extract solution data from boundary container
+      u_inner = get_node_vars(boundaries.u, equations_parabolic, dg, node_index, boundary_index)
+
+      # Outward-pointing normal direction (not normalized)
+      normal_direction = get_normal_direction(direction_index, contravariant_vectors,
+                                              i_node, j_node, element)
+
+      # TODO: revisit if we want more general boundary treatments.
+      # This assumes the gradient numerical flux at the boundary is the gradient variable,
+      # which is consistent with BR1, LDG.
+      flux_inner = u_inner
+
+      # Coordinates at boundary node
+      x = get_node_coords(node_coordinates, equations_parabolic, dg, i_node, j_node, element)
+
+      flux_ = boundary_condition_parabolic(flux_inner, u_inner, normal_direction, 
+                                           x, t, operator_type, equations_parabolic)
+
+      # Copy flux to element storage in the correct orientation
+      for v in eachvariable(equations_parabolic)
+        surface_flux_values[v, node_index, direction_index, element] = flux_[v]
+      end
+
+      i_node += i_node_step
+      j_node += j_node_step
+    end
+  end
+end
+

From 4ec387ee218b21d067e61fbdc7d915e1d04b3a00 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Mon, 5 Jun 2023 15:29:56 -0500
Subject: [PATCH 39/46] add tests

---
 ..._advection_diffusion_nonperiodic_curved.jl | 96 +++++++++++++++++++
 .../elixir_navierstokes_lid_driven_cavity.jl  | 82 ++++++++++++++++
 test/test_parabolic_2d.jl                     | 24 +++++
 3 files changed, 202 insertions(+)
 create mode 100644 examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_curved.jl
 create mode 100644 examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl

diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_curved.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_curved.jl
new file mode 100644
index 00000000000..55682f73fce
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_curved.jl
@@ -0,0 +1,96 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+diffusivity() = 5.0e-2
+advection_velocity = (1.0, 0.0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
+
+# Example setup taken from
+# - Truman Ellis, Jesse Chan, and Leszek Demkowicz (2016).
+#   Robust DPG methods for transient convection-diffusion.
+#   In: Building bridges: connections and challenges in modern approaches
+#   to numerical partial differential equations.
+#   [DOI](https://doi.org/10.1007/978-3-319-41640-3_6).
+function initial_condition_eriksson_johnson(x, t, equations)
+  l = 4
+  epsilon = diffusivity() # TODO: this requires epsilon < .6 due to sqrt
+  lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
+  lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
+  r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
+  s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
+  u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) +
+      cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1))
+  return SVector{1}(u)
+end
+initial_condition = initial_condition_eriksson_johnson
+
+boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition),
+                           :y_neg => BoundaryConditionDirichlet(initial_condition),
+                           :y_pos => BoundaryConditionDirichlet(initial_condition),
+                           :x_pos => boundary_condition_do_nothing)
+
+boundary_conditions_parabolic = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition), 
+                                     :x_pos => BoundaryConditionDirichlet(initial_condition), 
+                                     :y_neg => BoundaryConditionDirichlet(initial_condition), 
+                                     :y_pos => BoundaryConditionDirichlet(initial_condition))
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -0.5) 
+coordinates_max = ( 0.0,  0.5) 
+
+# This maps the domain [-1, 1]^2 to [-1, 0] x [-0.5, 0.5] while also 
+# introducing a curved warping to interior nodes. 
+function mapping(xi, eta)
+    x = xi  + 0.1 * sin(pi * xi) * sin(pi * eta)
+    y = eta + 0.1 * sin(pi * xi) * sin(pi * eta)
+    return SVector(0.5 * (1 + x) - 1, 0.5 * y)
+end 
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg=3, initial_refinement_level=2,
+                 mapping=mapping, periodicity=(false, false))
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver, 
+                                             boundary_conditions = (boundary_conditions, boundary_conditions_parabolic))
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+# The AliveCallback prints short status information in regular intervals
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+time_int_tol = 1.0e-11
+sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
+            ode_default_options()..., callback=callbacks)
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl
new file mode 100644
index 00000000000..051f4defe54
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl
@@ -0,0 +1,82 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+# TODO: parabolic; unify names of these accessor functions
+prandtl_number() = 0.72
+mu() = 0.001
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(),
+                                                          Prandtl=prandtl_number())
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
+
+# Create a uniformly refined mesh
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg=3, initial_refinement_level=2,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 periodicity=(false, false))
+
+function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D)
+  Ma = 0.1
+  rho = 1.0
+  u, v = 0.0, 0.0
+  p = 1.0 / (Ma^2 * equations.gamma)
+  return prim2cons(SVector(rho, u, v, p), equations)
+end
+initial_condition = initial_condition_cavity
+
+# BC types
+velocity_bc_lid = NoSlip((x, t, equations) -> SVector(1.0, 0.0))
+velocity_bc_cavity = NoSlip((x, t, equations) -> SVector(0.0, 0.0))
+heat_bc = Adiabatic((x, t, equations) -> 0.0)
+boundary_condition_lid = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc)
+boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc)
+
+# define periodic boundary conditions everywhere
+boundary_conditions = Dict( :x_neg => boundary_condition_slip_wall,
+                            :y_neg => boundary_condition_slip_wall,
+                            :y_pos => boundary_condition_slip_wall,
+                            :x_pos => boundary_condition_slip_wall)
+
+boundary_conditions_parabolic = Dict( :x_neg => boundary_condition_cavity,
+                                      :y_neg => boundary_condition_cavity,
+                                      :y_pos => boundary_condition_lid,
+                                      :x_pos => boundary_condition_cavity)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions=(boundary_conditions,
+                                                                  boundary_conditions_parabolic))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 25.0)
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=100)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+callbacks = CallbackSet(summary_callback, alive_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-8
+sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
+            ode_default_options()..., callback=callbacks)
+summary_callback() # print the timer summary
+
+
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index b0ac63d4ce9..8ab7bc7556b 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -200,6 +200,30 @@ isdir(outdir) && rm(outdir, recursive=true)
     )
   end
 
+  @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic_curved.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic_curved.jl"),
+      trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5),
+      l2 = [0.012380458938507371], 
+      linf = [0.10860506906472567]
+    )
+  end
+
+  @trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_curved.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_nonperiodic_curved.jl"),
+      trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5),
+      l2 = [0.04933902988507035], 
+      linf = [0.2550261714590271]
+    )
+  end
+
+  @trixi_testset "P4estMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_lid_driven_cavity.jl"),
+      initial_refinement_level = 2, tspan=(0.0, 0.5),
+      l2 = [0.00028716166408816073, 0.08101204560401647, 0.02099595625377768, 0.05008149754143295], 
+      linf = [0.014804500261322406, 0.9513271652357098, 0.7223919625994717, 1.4846907331004786]
+    )
+  end
+
 end
 
 # Clean up afterwards: delete Trixi.jl output directory

From 93df8686dd2f53231c6e9fe6db88337258ec7ff5 Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Tue, 6 Jun 2023 09:45:17 -0500
Subject: [PATCH 40/46] Update
 examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl

---
 examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
index 0b49dc9197a..1375f387cc1 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
@@ -176,7 +176,6 @@ end
 initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
-# TODO: parabolic. This is wrong! We got lucky because the velocity is actually zero at the boundaries. 
 velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3])
 heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom)

From b9bbd3971e3f0bb3f30271884f6ca10fe45f63d4 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Fri, 16 Jun 2023 11:16:42 -0500
Subject: [PATCH 41/46] formatting

---
 src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 299 ++++++++++-----------
 1 file changed, 135 insertions(+), 164 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
index e84e16dcbb2..73ac47ed1e3 100644
--- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl
@@ -405,57 +405,27 @@ function calc_interface_flux!(surface_flux_values,
             node_secondary_step = 1
         end
 
-  (; neighbor_ids, node_indices) = cache_parabolic.interfaces
-  index_range = eachnode(dg)
-  index_end = last(index_range)
-
-  @threaded for interface in eachinterface(dg, cache_parabolic)
-    # Get element and side index information on the primary element
-    primary_element = neighbor_ids[1, interface]
-    primary_indices = node_indices[1, interface]
-    primary_direction_index = indices2direction(primary_indices)
-
-    # Create the local i,j indexing on the primary element used to pull normal direction information
-    i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
-    j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)
-
-    i_primary = i_primary_start
-    j_primary = j_primary_start
-
-    # Get element and side index information on the secondary element
-    secondary_element = neighbor_ids[2, interface]
-    secondary_indices = node_indices[2, interface]
-    secondary_direction_index = indices2direction(secondary_indices)
-
-    # Initiate the secondary index to be used in the surface for loop.
-    # This index on the primary side will always run forward but
-    # the secondary index might need to run backwards for flipped sides.
-    if :i_backward in secondary_indices
-      node_secondary = index_end
-      node_secondary_step = -1
-    else
-      node_secondary = 1
-      node_secondary_step = 1
-    end
-
-    for node in eachnode(dg)
-      # We prolong the viscous flux dotted with respect the outward normal on the 
-      # primary element. We assume a BR-1 type of flux.
-      viscous_flux_normal_ll, viscous_flux_normal_rr = 
-        get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, dg, node, interface)
+        for node in eachnode(dg)
+            # We prolong the viscous flux dotted with respect the outward normal on the 
+            # primary element. We assume a BR-1 type of flux.
+            viscous_flux_normal_ll, viscous_flux_normal_rr = get_surface_node_vars(cache_parabolic.interfaces.u,
+                                                                                   equations_parabolic,
+                                                                                   dg, node,
+                                                                                   interface)
 
-      flux = 0.5 * (viscous_flux_normal_ll + viscous_flux_normal_rr)
+            flux = 0.5 * (viscous_flux_normal_ll + viscous_flux_normal_rr)
 
-      for v in eachvariable(equations_parabolic)
-        surface_flux_values[v, node, primary_direction_index, primary_element] = flux[v]
-        surface_flux_values[v, node_secondary, secondary_direction_index, secondary_element] = -flux[v]
-      end
+            for v in eachvariable(equations_parabolic)
+                surface_flux_values[v, node, primary_direction_index, primary_element] = flux[v]
+                surface_flux_values[v, node_secondary, secondary_direction_index, secondary_element] = -flux[v]
+            end
 
-      # Increment primary element indices to pull the normal direction
-      i_primary += i_primary_step
-      j_primary += j_primary_step
-      # Increment the surface node index along the secondary element
-      node_secondary += node_secondary_step
+            # Increment primary element indices to pull the normal direction
+            i_primary += i_primary_step
+            j_primary += j_primary_step
+            # Increment the surface node index along the secondary element
+            node_secondary += node_secondary_step
+        end
     end
 
     return nothing
@@ -466,156 +436,157 @@ function prolong2boundaries!(cache_parabolic, flux_viscous,
                              mesh::P4estMesh{2},
                              equations_parabolic::AbstractEquationsParabolic,
                              surface_integral, dg::DG, cache)
-  (; boundaries) = cache_parabolic
-  (; contravariant_vectors) = cache_parabolic.elements
-  index_range = eachnode(dg)
-
-  flux_viscous_x, flux_viscous_y = flux_viscous
-
-  @threaded for boundary in eachboundary(dg, cache_parabolic)
-    # Copy solution data from the element using "delayed indexing" with
-    # a start value and a step size to get the correct face and orientation.
-    element       = boundaries.neighbor_ids[boundary]
-    node_indices  = boundaries.node_indices[boundary]
-    direction     = indices2direction(node_indices)
-
-    i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
-    j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
-
-    i_node = i_node_start
-    j_node = j_node_start
-    for i in eachnode(dg)
-      # this is the outward normal direction on the primary element
-      normal_direction = get_normal_direction(direction, contravariant_vectors,
-                                              i_node, j_node, element)
-
-      for v in eachvariable(equations_parabolic)
-        flux_viscous = SVector(flux_viscous_x[v, i_node, j_node, element], 
-                               flux_viscous_y[v, i_node, j_node, element])
-
-        boundaries.u[v, i, boundary] = dot(flux_viscous, normal_direction)
-      end
-      i_node += i_node_step
-      j_node += j_node_step
+    (; boundaries) = cache_parabolic
+    (; contravariant_vectors) = cache_parabolic.elements
+    index_range = eachnode(dg)
 
-    end
+    flux_viscous_x, flux_viscous_y = flux_viscous
 
+    @threaded for boundary in eachboundary(dg, cache_parabolic)
+        # Copy solution data from the element using "delayed indexing" with
+        # a start value and a step size to get the correct face and orientation.
+        element = boundaries.neighbor_ids[boundary]
+        node_indices = boundaries.node_indices[boundary]
+        direction = indices2direction(node_indices)
+
+        i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+        j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+        i_node = i_node_start
+        j_node = j_node_start
+        for i in eachnode(dg)
+            # this is the outward normal direction on the primary element
+            normal_direction = get_normal_direction(direction, contravariant_vectors,
+                                                    i_node, j_node, element)
+
+            for v in eachvariable(equations_parabolic)
+                flux_viscous = SVector(flux_viscous_x[v, i_node, j_node, element],
+                                       flux_viscous_y[v, i_node, j_node, element])
+
+                boundaries.u[v, i, boundary] = dot(flux_viscous, normal_direction)
+            end
+            i_node += i_node_step
+            j_node += j_node_step
+        end
+    end
     return nothing
 end
 
-function calc_boundary_flux_gradients!(cache, t, 
-                                       boundary_condition::Union{BoundaryConditionPeriodic, BoundaryConditionDoNothing}, 
+function calc_boundary_flux_gradients!(cache, t,
+                                       boundary_condition::Union{BoundaryConditionPeriodic,
+                                                                 BoundaryConditionDoNothing
+                                                                 },
                                        mesh::P4estMesh, equations, surface_integral, dg::DG)
-  @assert isempty(eachboundary(dg, cache))
+    @assert isempty(eachboundary(dg, cache))
 end
 
 # Function barrier for type stability
 function calc_boundary_flux_gradients!(cache, t, boundary_conditions, mesh::P4estMesh,
                                        equations, surface_integral, dg::DG)
-  (; boundary_condition_types, boundary_indices) = boundary_conditions
+    (; boundary_condition_types, boundary_indices) = boundary_conditions
 
-  calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices,
-                              Gradient(), mesh, equations, surface_integral, dg)
-  return nothing
+    calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices,
+                                Gradient(), mesh, equations, surface_integral, dg)
+    return nothing
 end
 
 function calc_boundary_flux_divergence!(cache, t, boundary_conditions, mesh::P4estMesh,
                                         equations, surface_integral, dg::DG)
-  (; boundary_condition_types, boundary_indices) = boundary_conditions
+    (; boundary_condition_types, boundary_indices) = boundary_conditions
 
-  calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices,
-                              Divergence(), mesh, equations, surface_integral, dg)
-  return nothing
+    calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices,
+                                Divergence(), mesh, equations, surface_integral, dg)
+    return nothing
 end
 
 # Iterate over tuples of boundary condition types and associated indices
 # in a type-stable way using "lispy tuple programming".
 function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N, Any},
                                      BC_indices::NTuple{N, Vector{Int}},
-                                     operator_type, 
+                                     operator_type,
                                      mesh::P4estMesh,
                                      equations, surface_integral, dg::DG) where {N}
-  # Extract the boundary condition type and index vector
-  boundary_condition = first(BCs)
-  boundary_condition_indices = first(BC_indices)
-  # Extract the remaining types and indices to be processed later
-  remaining_boundary_conditions = Base.tail(BCs)
-  remaining_boundary_condition_indices = Base.tail(BC_indices)
-
-  # process the first boundary condition type
-  calc_boundary_flux!(cache, t, boundary_condition, boundary_condition_indices,
-                      operator_type, mesh, equations, surface_integral, dg)
-
-  # recursively call this method with the unprocessed boundary types
-  calc_boundary_flux_by_type!(cache, t, remaining_boundary_conditions,
-                              remaining_boundary_condition_indices,
-                              operator_type, 
-                              mesh, equations, surface_integral, dg)
-
-  return nothing
+    # Extract the boundary condition type and index vector
+    boundary_condition = first(BCs)
+    boundary_condition_indices = first(BC_indices)
+    # Extract the remaining types and indices to be processed later
+    remaining_boundary_conditions = Base.tail(BCs)
+    remaining_boundary_condition_indices = Base.tail(BC_indices)
+
+    # process the first boundary condition type
+    calc_boundary_flux!(cache, t, boundary_condition, boundary_condition_indices,
+                        operator_type, mesh, equations, surface_integral, dg)
+
+    # recursively call this method with the unprocessed boundary types
+    calc_boundary_flux_by_type!(cache, t, remaining_boundary_conditions,
+                                remaining_boundary_condition_indices,
+                                operator_type,
+                                mesh, equations, surface_integral, dg)
+
+    return nothing
 end
 
 # terminate the type-stable iteration over tuples
 function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{},
-                                     operator_type, mesh::P4estMesh, equations, 
+                                     operator_type, mesh::P4estMesh, equations,
                                      surface_integral, dg::DG)
-  nothing
+    nothing
 end
 
-
-function calc_boundary_flux!(cache, t, 
+function calc_boundary_flux!(cache, t,
                              boundary_condition_parabolic, # works with Dict types
-                             boundary_condition_indices, 
-                             operator_type,  mesh::P4estMesh{2}, 
+                             boundary_condition_indices,
+                             operator_type, mesh::P4estMesh{2},
                              equations_parabolic::AbstractEquationsParabolic,
                              surface_integral, dg::DG)
-  (; boundaries) = cache
-  (; node_coordinates, surface_flux_values) = cache.elements
-  (; contravariant_vectors) = cache.elements
-  index_range = eachnode(dg)
-
-  @threaded for local_index in eachindex(boundary_condition_indices)
-    # Use the local index to get the global boundary index from the pre-sorted list
-    boundary_index = boundary_condition_indices[local_index]
-
-    # Get information on the adjacent element, compute the surface fluxes,
-    # and store them
-    element         = boundaries.neighbor_ids[boundary_index]
-    node_indices    = boundaries.node_indices[boundary_index]
-    direction_index = indices2direction(node_indices)
-
-    i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
-    j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
-
-    i_node = i_node_start
-    j_node = j_node_start
-    for node_index in eachnode(dg)
-      # Extract solution data from boundary container
-      u_inner = get_node_vars(boundaries.u, equations_parabolic, dg, node_index, boundary_index)
-
-      # Outward-pointing normal direction (not normalized)
-      normal_direction = get_normal_direction(direction_index, contravariant_vectors,
-                                              i_node, j_node, element)
-
-      # TODO: revisit if we want more general boundary treatments.
-      # This assumes the gradient numerical flux at the boundary is the gradient variable,
-      # which is consistent with BR1, LDG.
-      flux_inner = u_inner
-
-      # Coordinates at boundary node
-      x = get_node_coords(node_coordinates, equations_parabolic, dg, i_node, j_node, element)
-
-      flux_ = boundary_condition_parabolic(flux_inner, u_inner, normal_direction, 
-                                           x, t, operator_type, equations_parabolic)
-
-      # Copy flux to element storage in the correct orientation
-      for v in eachvariable(equations_parabolic)
-        surface_flux_values[v, node_index, direction_index, element] = flux_[v]
-      end
-
-      i_node += i_node_step
-      j_node += j_node_step
+    (; boundaries) = cache
+    (; node_coordinates, surface_flux_values) = cache.elements
+    (; contravariant_vectors) = cache.elements
+    index_range = eachnode(dg)
+
+    @threaded for local_index in eachindex(boundary_condition_indices)
+        # Use the local index to get the global boundary index from the pre-sorted list
+        boundary_index = boundary_condition_indices[local_index]
+
+        # Get information on the adjacent element, compute the surface fluxes,
+        # and store them
+        element = boundaries.neighbor_ids[boundary_index]
+        node_indices = boundaries.node_indices[boundary_index]
+        direction_index = indices2direction(node_indices)
+
+        i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+        j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+        i_node = i_node_start
+        j_node = j_node_start
+        for node_index in eachnode(dg)
+            # Extract solution data from boundary container
+            u_inner = get_node_vars(boundaries.u, equations_parabolic, dg, node_index,
+                                    boundary_index)
+
+            # Outward-pointing normal direction (not normalized)
+            normal_direction = get_normal_direction(direction_index, contravariant_vectors,
+                                                    i_node, j_node, element)
+
+            # TODO: revisit if we want more general boundary treatments.
+            # This assumes the gradient numerical flux at the boundary is the gradient variable,
+            # which is consistent with BR1, LDG.
+            flux_inner = u_inner
+
+            # Coordinates at boundary node
+            x = get_node_coords(node_coordinates, equations_parabolic, dg, i_node, j_node,
+                                element)
+
+            flux_ = boundary_condition_parabolic(flux_inner, u_inner, normal_direction,
+                                                 x, t, operator_type, equations_parabolic)
+
+            # Copy flux to element storage in the correct orientation
+            for v in eachvariable(equations_parabolic)
+                surface_flux_values[v, node_index, direction_index, element] = flux_[v]
+            end
+
+            i_node += i_node_step
+            j_node += j_node_step
+        end
     end
-  end
 end
-

From dbddf8586bec800e2d40906a586870cb5417f942 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Mon, 19 Jun 2023 09:17:01 -0500
Subject: [PATCH 42/46] fix CNS convergence elixir and add to tests

---
 .../elixir_navierstokes_convergence.jl        | 20 +++++--------------
 test/test_parabolic_2d.jl                     |  8 ++++++++
 2 files changed, 13 insertions(+), 15 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
index 1375f387cc1..8111df8251a 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
@@ -18,17 +18,11 @@ solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
 coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
 coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
 
-# Create a uniformly refined mesh with periodic boundaries
-mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level=4,
-                periodicity=(true, true),
-                n_cells_max=30_000) # set maximum capacity of tree data structure
-
 trees_per_dimension = (4, 4)
 mesh = P4estMesh(trees_per_dimension,
                  polydeg=3, initial_refinement_level=2,
                  coordinates_min=coordinates_min, coordinates_max=coordinates_max,
-                 periodicity=true)
+                 periodicity=(true, false))
 
 # Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D`
 #       since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`)
@@ -181,16 +175,12 @@ heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom)
 
 # define inviscid boundary conditions
-boundary_conditions = (; x_neg = boundary_condition_periodic,
-                         x_pos = boundary_condition_periodic,
-                         y_neg = boundary_condition_slip_wall,
-                         y_pos = boundary_condition_slip_wall)
+boundary_conditions = Dict(:y_neg => boundary_condition_slip_wall,
+                           :y_pos => boundary_condition_slip_wall)
 
 # define viscous boundary conditions
-boundary_conditions_parabolic = (; x_neg = boundary_condition_periodic,
-                                   x_pos = boundary_condition_periodic,
-                                   y_neg = boundary_condition_top_bottom,
-                                   y_pos = boundary_condition_top_bottom)
+boundary_conditions_parabolic = Dict(:y_neg => boundary_condition_top_bottom,
+                                     :y_pos => boundary_condition_top_bottom)
 
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver;
                                              boundary_conditions=(boundary_conditions, boundary_conditions_parabolic),
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 8ab7bc7556b..6d953f8372f 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -216,6 +216,14 @@ isdir(outdir) && rm(outdir, recursive=true)
     )
   end
 
+  @trixi_testset "P4estMesh2D: elixir_navierstokes_convergence.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_convergence.jl"),
+      initial_refinement_level = 2, tspan=(0.0, 0.5),
+      l2 = [4.2852610204634015e-6, 2.0679143814141558e-5, 1.4336134288439704e-5, 7.800744943199459e-6], 
+      linf = [3.570853695711307e-5, 0.00034481633921521794, 8.593188464759635e-5, 0.00010656913331352769]
+    )
+  end
+
   @trixi_testset "P4estMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin
     @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_lid_driven_cavity.jl"),
       initial_refinement_level = 2, tspan=(0.0, 0.5),

From 64488b2cb05a39177c1f99681685d0f1afeb19ca Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Mon, 19 Jun 2023 10:54:35 -0500
Subject: [PATCH 43/46] update test values

---
 test/test_parabolic_2d.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 6d953f8372f..471b976e990 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -218,9 +218,9 @@ isdir(outdir) && rm(outdir, recursive=true)
 
   @trixi_testset "P4estMesh2D: elixir_navierstokes_convergence.jl" begin
     @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_convergence.jl"),
-      initial_refinement_level = 2, tspan=(0.0, 0.5),
-      l2 = [4.2852610204634015e-6, 2.0679143814141558e-5, 1.4336134288439704e-5, 7.800744943199459e-6], 
-      linf = [3.570853695711307e-5, 0.00034481633921521794, 8.593188464759635e-5, 0.00010656913331352769]
+      initial_refinement_level = 1, tspan=(0.0, 0.2), 
+      l2 = [0.0003811978985836709, 0.0005874314969169538, 0.0009142898787923481, 0.0011613918899727263], 
+      linf = [0.0021633623982135752, 0.009484348274135372, 0.004231572066492217, 0.011661660275365193]
     )
   end
 

From 47d15cad7f67c06d511b67b68634e8b4be7a72ab Mon Sep 17 00:00:00 2001
From: Ahmad Peyvan <ahmad_peyvan@brown.edu>
Date: Thu, 22 Jun 2023 15:32:04 -0400
Subject: [PATCH 44/46] adding Dirichlet BCs for parabolic terms

---
 ...ir_navierstokes_convergence_nonperiodic.jl | 215 ++++++++++++++++++
 .../compressible_navier_stokes_2d.jl          |  28 +++
 2 files changed, 243 insertions(+)
 create mode 100644 examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl
new file mode 100644
index 00000000000..935f132ba4b
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl
@@ -0,0 +1,215 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+prandtl_number() = 0.72
+mu() = 0.01
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(),
+                                                          gradient_variables=GradientVariablesPrimitive())
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
+               volume_integral=VolumeIntegralWeakForm())
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg=3, initial_refinement_level=2,
+                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
+                 periodicity=(false, false))
+
+# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D`
+#       since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`)
+#       and by the initial condition (which passes in `CompressibleEulerEquations2D`).
+# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000)
+function initial_condition_navier_stokes_convergence_test(x, t, equations)
+  # Amplitude and shift
+  A = 0.5
+  c = 2.0
+
+  # convenience values for trig. functions
+  pi_x = pi * x[1]
+  pi_y = pi * x[2]
+  pi_t = pi * t
+
+  rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+  v1  = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t)
+  v2  = v1
+  p   = rho^2
+
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations)
+  y = x[2]
+
+  # TODO: parabolic
+  # we currently need to hardcode these parameters until we fix the "combined equation" issue
+  # see also https://github.com/trixi-framework/Trixi.jl/pull/1160
+  inv_gamma_minus_one = inv(equations.gamma - 1)
+  Pr = prandtl_number()
+  mu_ = mu()
+
+  # Same settings as in `initial_condition`
+  # Amplitude and shift
+  A = 0.5
+  c = 2.0
+
+  # convenience values for trig. functions
+  pi_x = pi * x[1]
+  pi_y = pi * x[2]
+  pi_t = pi * t
+
+  # compute the manufactured solution and all necessary derivatives
+  rho    =  c  + A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+  rho_t  = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t)
+  rho_x  =  pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t)
+  rho_y  = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t)
+  rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+  rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t)
+
+  v1    =       sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
+  v1_t  = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t)
+  v1_x  =  pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
+  v1_y  =       sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t)
+  v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
+  v1_xy =  pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t)
+  v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0)
+                         - A * A * log(y + 2.0) * exp(-A * (y - 1.0))
+                         - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t))
+  v2    = v1
+  v2_t  = v1_t
+  v2_x  = v1_x
+  v2_y  = v1_y
+  v2_xx = v1_xx
+  v2_xy = v1_xy
+  v2_yy = v1_yy
+
+  p    = rho * rho
+  p_t  = 2.0 * rho * rho_t
+  p_x  = 2.0 * rho * rho_x
+  p_y  = 2.0 * rho * rho_y
+  p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x
+  p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y
+
+  # Note this simplifies slightly because the ansatz assumes that v1 = v2
+  E   = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2)
+  E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t
+  E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x
+  E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y
+
+  # Some convenience constants
+  T_const = equations.gamma * inv_gamma_minus_one / Pr
+  inv_rho_cubed = 1.0 / (rho^3)
+
+  # compute the source terms
+  # density equation
+  du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y
+
+  # x-momentum equation
+  du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2
+                                        + 2.0   * rho  * v1 * v1_x
+                                        + rho_y * v1   * v2
+                                        + rho   * v1_y * v2
+                                        + rho   * v1   * v2_y
+    # stress tensor from x-direction
+                      - 4.0 / 3.0 * v1_xx * mu_
+                      + 2.0 / 3.0 * v2_xy * mu_
+                      - v1_yy             * mu_
+                      - v2_xy             * mu_ )
+  # y-momentum equation
+  du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1    * v2
+                                        + rho   * v1_x  * v2
+                                        + rho   * v1    * v2_x
+                                        +         rho_y * v2^2
+                                        + 2.0   * rho   * v2 * v2_y
+    # stress tensor from y-direction
+                      - v1_xy             * mu_
+                      - v2_xx             * mu_
+                      - 4.0 / 3.0 * v2_yy * mu_
+                      + 2.0 / 3.0 * v1_xy * mu_ )
+  # total energy equation
+  du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x)
+              + v2_y * (E + p) + v2 * (E_y + p_y)
+    # stress tensor and temperature gradient terms from x-direction
+                                - 4.0 / 3.0 * v1_xx * v1   * mu_
+                                + 2.0 / 3.0 * v2_xy * v1   * mu_
+                                - 4.0 / 3.0 * v1_x  * v1_x * mu_
+                                + 2.0 / 3.0 * v2_y  * v1_x * mu_
+                                - v1_xy     * v2           * mu_
+                                - v2_xx     * v2           * mu_
+                                - v1_y      * v2_x         * mu_
+                                - v2_x      * v2_x         * mu_
+         - T_const * inv_rho_cubed * (        p_xx * rho   * rho
+                                      - 2.0 * p_x  * rho   * rho_x
+                                      + 2.0 * p    * rho_x * rho_x
+                                      -       p    * rho   * rho_xx ) * mu_
+    # stress tensor and temperature gradient terms from y-direction
+                                - v1_yy     * v1           * mu_
+                                - v2_xy     * v1           * mu_
+                                - v1_y      * v1_y         * mu_
+                                - v2_x      * v1_y         * mu_
+                                - 4.0 / 3.0 * v2_yy * v2   * mu_
+                                + 2.0 / 3.0 * v1_xy * v2   * mu_
+                                - 4.0 / 3.0 * v2_y  * v2_y * mu_
+                                + 2.0 / 3.0 * v1_x  * v2_y * mu_
+         - T_const * inv_rho_cubed * (        p_yy * rho   * rho
+                                      - 2.0 * p_y  * rho   * rho_y
+                                      + 2.0 * p    * rho_y * rho_y
+                                      -       p    * rho   * rho_yy ) * mu_ )
+
+  return SVector(du1, du2, du3, du4)
+end
+
+initial_condition = initial_condition_navier_stokes_convergence_test
+
+# BC types
+velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3])
+heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
+boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom)
+
+boundary_condition_left_right = BoundaryConditionDirichlet(initial_condition_navier_stokes_convergence_test)
+
+# define inviscid boundary conditions
+boundary_conditions = Dict(:x_neg => boundary_condition_left_right,
+                           :x_pos => boundary_condition_left_right,
+                           :y_neg => boundary_condition_slip_wall,
+                           :y_pos => boundary_condition_slip_wall)
+
+# define viscous boundary conditions
+boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_left_right,
+                                     :x_pos => boundary_condition_left_right,
+                                     :y_neg => boundary_condition_top_bottom,
+                                     :y_pos => boundary_condition_top_bottom)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver;
+                                             boundary_conditions=(boundary_conditions, boundary_conditions_parabolic),
+                                             source_terms=source_terms_navier_stokes_convergence_test)
+
+# ###############################################################################
+# # ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-8
+sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5,
+            ode_default_options()..., callback=callbacks)
+summary_callback() # print the timer summary
+
diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl
index 9b06e0b5abf..3bbee1b7556 100644
--- a/src/equations/compressible_navier_stokes_2d.jl
+++ b/src/equations/compressible_navier_stokes_2d.jl
@@ -387,6 +387,8 @@ end
     return SVector(flux_inner[1], flux_inner[2], flux_inner[3], normal_energy_flux)
 end
 
+
+
 @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                         <:Isothermal})(flux_inner,
                                                                                        u_inner,
@@ -489,3 +491,29 @@ end
                                                                                                                                       })
     return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4])
 end
+
+# Diriclet Boundary Condition for P4est
+@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner,
+                                                                        u_inner,
+                                                                        normal::AbstractVector,
+                                                                        x, t,
+                                                                        operator_type::Gradient,
+                                                                        equations::CompressibleNavierStokesDiffusion2D{
+                                                                                                                        GradientVariablesPrimitive
+                                                                                                                        })
+   # We have to convert to primitive variables here since the "boundary_value_function" returns conservative variables
+   u_boundary = boundary_condition.boundary_value_function(x, t, equations)
+
+   return cons2prim(u_boundary,equations)
+end
+
+@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner,
+                                                                                      u_inner,
+                                                                                      normal::AbstractVector,
+                                                                                      x, t,
+                                                                                      operator_type::Divergence,
+                                                                                      equations::CompressibleNavierStokesDiffusion2D{
+                                                                                                                                     GradientVariablesPrimitive
+                                                                                                                                     })
+    return flux_inner
+end
\ No newline at end of file

From 5afbdb7cd1983f6063e8d53f737a06f270532474 Mon Sep 17 00:00:00 2001
From: Ahmad Peyvan <ahmad_peyvan@brown.edu>
Date: Thu, 22 Jun 2023 15:35:57 -0400
Subject: [PATCH 45/46] add test

---
 test/test_parabolic_2d.jl | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 471b976e990..ec5d7de2b7a 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -224,6 +224,14 @@ isdir(outdir) && rm(outdir, recursive=true)
     )
   end
 
+  @trixi_testset "P4estMesh2D: elixir_navierstokes_convergence_nonperiodic.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_convergence_nonperiodic.jl"),
+      initial_refinement_level = 1, tspan=(0.0, 0.2), 
+      l2 = [0.00040364962558511795, 0.0005869762481506936, 0.00091488537427274, 0.0011984191566376762], 
+      linf = [0.0024993634941723464, 0.009487866203944725, 0.004505829506628117, 0.011634902776245681]
+    )
+  end
+
   @trixi_testset "P4estMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin
     @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_lid_driven_cavity.jl"),
       initial_refinement_level = 2, tspan=(0.0, 0.5),

From 946c3ce8b9de313c6dc99793405ded98a39bae68 Mon Sep 17 00:00:00 2001
From: Ahmad Peyvan <ahmad_peyvan@brown.edu>
Date: Thu, 22 Jun 2023 15:53:56 -0400
Subject: [PATCH 46/46] Dirichlet BCs

---
 .../compressible_navier_stokes_2d.jl          | 40 +++++++++----------
 1 file changed, 19 insertions(+), 21 deletions(-)

diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl
index 3bbee1b7556..ac5013f6fc6 100644
--- a/src/equations/compressible_navier_stokes_2d.jl
+++ b/src/equations/compressible_navier_stokes_2d.jl
@@ -387,8 +387,6 @@ end
     return SVector(flux_inner[1], flux_inner[2], flux_inner[3], normal_energy_flux)
 end
 
-
-
 @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                         <:Isothermal})(flux_inner,
                                                                                        u_inner,
@@ -494,26 +492,26 @@ end
 
 # Diriclet Boundary Condition for P4est
 @inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner,
-                                                                        u_inner,
-                                                                        normal::AbstractVector,
-                                                                        x, t,
-                                                                        operator_type::Gradient,
-                                                                        equations::CompressibleNavierStokesDiffusion2D{
-                                                                                                                        GradientVariablesPrimitive
-                                                                                                                        })
-   # We have to convert to primitive variables here since the "boundary_value_function" returns conservative variables
-   u_boundary = boundary_condition.boundary_value_function(x, t, equations)
-
-   return cons2prim(u_boundary,equations)
+                                                                  u_inner,
+                                                                  normal::AbstractVector,
+                                                                  x, t,
+                                                                  operator_type::Gradient,
+                                                                  equations::CompressibleNavierStokesDiffusion2D{
+                                                                                                                 GradientVariablesPrimitive
+                                                                                                                 })
+    # We have to convert to primitive variables here since the "boundary_value_function" returns conservative variables
+    u_boundary = boundary_condition.boundary_value_function(x, t, equations)
+
+    return cons2prim(u_boundary, equations)
 end
 
 @inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner,
-                                                                                      u_inner,
-                                                                                      normal::AbstractVector,
-                                                                                      x, t,
-                                                                                      operator_type::Divergence,
-                                                                                      equations::CompressibleNavierStokesDiffusion2D{
-                                                                                                                                     GradientVariablesPrimitive
-                                                                                                                                     })
+                                                                  u_inner,
+                                                                  normal::AbstractVector,
+                                                                  x, t,
+                                                                  operator_type::Divergence,
+                                                                  equations::CompressibleNavierStokesDiffusion2D{
+                                                                                                                 GradientVariablesPrimitive
+                                                                                                                 })
     return flux_inner
-end
\ No newline at end of file
+end