From 8bc5469b3d883be0ce088386e1ec794e448003ef Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Wed, 11 Oct 2023 06:22:49 -0500
Subject: [PATCH] fix bug in `apply_jacobian_parabolic!` for `P4estMesh`
 (#1668)

* fix bug in parabolic jacobian computations

* formatting

* updating test results and removing duplicate test

---------

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/solvers/dgsem_tree/dg_2d_parabolic.jl | 20 +++++++++++++++++++-
 src/solvers/dgsem_tree/dg_3d_parabolic.jl | 23 +++++++++++++++++++++--
 test/test_parabolic_2d.jl                 | 16 ++++------------
 3 files changed, 44 insertions(+), 15 deletions(-)

diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
index 3dbc55412ad..1c32703c7c3 100644
--- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl
@@ -935,7 +935,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_parabolic!(du, mesh::Union{TreeMesh{2}, P4estMesh{2}},
+function apply_jacobian_parabolic!(du, mesh::TreeMesh{2},
                                    equations::AbstractEquationsParabolic, dg::DG, cache)
     @unpack inverse_jacobian = cache.elements
 
@@ -951,4 +951,22 @@ function apply_jacobian_parabolic!(du, mesh::Union{TreeMesh{2}, P4estMesh{2}},
 
     return nothing
 end
+
+function apply_jacobian_parabolic!(du, mesh::P4estMesh{2},
+                                   equations::AbstractEquationsParabolic,
+                                   dg::DG, cache)
+    @unpack inverse_jacobian = cache.elements
+
+    @threaded for element in eachelement(dg, cache)
+        for j in eachnode(dg), i in eachnode(dg)
+            factor = inverse_jacobian[i, j, element]
+
+            for v in eachvariable(equations)
+                du[v, i, j, element] *= factor
+            end
+        end
+    end
+
+    return nothing
+end
 end # @muladd
diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl
index 9817e0e5f0e..37492dbcb91 100644
--- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl
+++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl
@@ -1125,8 +1125,9 @@ 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_parabolic!(du, mesh::Union{TreeMesh{3}, P4estMesh{3}},
-                                   equations::AbstractEquationsParabolic, dg::DG, cache)
+function apply_jacobian_parabolic!(du, mesh::TreeMesh{3},
+                                   equations::AbstractEquationsParabolic,
+                                   dg::DG, cache)
     @unpack inverse_jacobian = cache.elements
 
     @threaded for element in eachelement(dg, cache)
@@ -1141,4 +1142,22 @@ function apply_jacobian_parabolic!(du, mesh::Union{TreeMesh{3}, P4estMesh{3}},
 
     return nothing
 end
+
+function apply_jacobian_parabolic!(du, mesh::P4estMesh{3},
+                                   equations::AbstractEquationsParabolic,
+                                   dg::DG, cache)
+    @unpack inverse_jacobian = cache.elements
+
+    @threaded for element in eachelement(dg, cache)
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            factor = inverse_jacobian[i, j, k, element]
+
+            for v in eachvariable(equations)
+                du[v, i, j, k, element] *= factor
+            end
+        end
+    end
+
+    return nothing
+end
 end # @muladd
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 967aa1069a0..98aa337afb4 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -260,24 +260,16 @@ isdir(outdir) && rm(outdir, recursive=true)
   @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_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]
+      l2 = [0.006708147442490916], 
+      linf = [0.04807038397976693]
     )
   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]
+      l2 = [0.00919917034843865], 
+      linf = [0.14186297438393505]
     )
   end