From 673a72d7cc7463a527082dba94508417672aeddc Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Thu, 12 Sep 2024 18:08:07 +0200
Subject: [PATCH] fix entropy_modified

---
 src/equations/bbm_1d.jl | 14 ++++++++------
 test/test_bbm_1d.jl     | 12 ++++++------
 2 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/src/equations/bbm_1d.jl b/src/equations/bbm_1d.jl
index 2ee56bbb..4c6ab528 100644
--- a/src/equations/bbm_1d.jl
+++ b/src/equations/bbm_1d.jl
@@ -96,7 +96,8 @@ function create_cache(mesh, equations::BBMEquation1D,
     eta2_x = zero(eta2)
     eta_x = zero(eta2)
     etaeta_x = zero(eta2)
-    cache = (; invImD2, eta2, eta2_x, eta_x, etaeta_x, solver.D1)
+    eta_xx = zero(eta2)
+    cache = (; invImD2, eta2, eta2_x, eta_x, etaeta_x, eta_xx, solver.D1, solver.D2)
     if solver.D1 isa PeriodicUpwindOperators
         eta_x_upwind = zero(eta2)
         cache = (; cache..., eta_x_upwind)
@@ -158,7 +159,7 @@ is a conserved quantity (for periodic boundary conditions).
 
 It is given by
 ```math
-\\frac{1}{2} \\eta^2 + \\frac{1}{2} \\eta_x^2.
+\\frac{1}{2} \\eta(\\eta - \\eta_{xx}).
 ```
 
 `q_global` is a vector of the primitive variables at ALL nodes.
@@ -169,14 +170,15 @@ See also [`energy_total_modified`](@ref).
 function energy_total_modified!(e, q_global, equations::BBMEquation1D, cache)
     eta, = q_global.x
 
-    (; D1, eta_x) = cache
+    (; D1, D2, eta_xx, tmp1) = cache
     if D1 isa PeriodicUpwindOperators
-        mul!(eta_x, D1.central, eta)
+        mul!(tmp1, D1.minus, eta)
+        mul!(eta_xx, D1.plus, tmp1)
     else
-        mul!(eta_x, D1, eta)
+        mul!(eta_xx, D2, eta)
     end
 
-    @.. e = 0.5 * (eta^2 + eta_x^2)
+    @.. e = 0.5 * eta * (eta - eta_xx)
     return nothing
 end
 
diff --git a/test/test_bbm_1d.jl b/test/test_bbm_1d.jl
index a5967649..e188207f 100644
--- a/test/test_bbm_1d.jl
+++ b/test/test_bbm_1d.jl
@@ -15,7 +15,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_1d")
                             linf=[9.756568022900591e-5],
                             cons_error=[7.105427357601002e-15],
                             change_waterheight=-7.105427357601002e-15,
-                            change_entropy_modified=-4.4198182447274803e-7,
+                            change_entropy_modified=-4.4274395039067826e-7,
                             change_invariant_cubic=-3.198066679033218e-6)
 
         @test_allocations(semi, sol, allocs=5_000)
@@ -55,12 +55,12 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_1d")
     @trixi_testset "bbm_1d_relaxation" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_1d_relaxation.jl"),
                             tspan=(0.0, 100.0),
-                            l2=[0.00025104961489025117],
-                            linf=[9.745906869124132e-5],
+                            l2=[0.0002510482044505331],
+                            linf=[9.745861147153478e-5],
                             cons_error=[5.329070518200751e-15],
-                            change_waterheight=-5.329070518200751e-15,
-                            change_entropy_modified=-2.220446049250313e-16,
-                            change_invariant_cubic=-1.5787975371495122e-8)
+                            change_waterheight=5.329070518200751e-15,
+                            change_entropy_modified=0.0,
+                            change_invariant_cubic=-1.0308752962373546e-8)
 
         @test_allocations(semi, sol, allocs=5_000)
     end