scale_by_mass_matrix! outside of solve_system_matrix!

JoshuaLampert · Dec 11, 2024 · 3fa969f · 3fa969f
1 parent ca345a5
commit 3fa969f
Show file tree

Hide file tree

Showing 4 changed files with 18 additions and 14 deletions.
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
@@ -565,17 +565,15 @@ function solve_system_matrix!(dv, system_matrix, rhs,
         (; factorization) = cache
         cholesky!(factorization, system_matrix; check = false)
         if issuccess(factorization)
-            scale_by_mass_matrix!(rhs, D1)
             # see https://github.com/JoshuaLampert/DispersiveShallowWater.jl/issues/122
-            dv .= factorization \ rhs[2:(end - 1)]
+            dv .= factorization \ rhs
         else
             # The factorization may fail if the time step is too large
             # and h becomes negative.
             fill!(dv, NaN)
         end
     else
         factorization = cholesky!(system_matrix)
-        scale_by_mass_matrix!(rhs, D1)
         ldiv!(dv, factorization, rhs)
     end
     return nothing

diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl
@@ -399,6 +399,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla
     end
 
     @trixi_timeit timer() "solving elliptic system" begin
+        scale_by_mass_matrix!(tmp, D1)
         solve_system_matrix!(dv, system_matrix, tmp,
                              equations, D1, cache)
     end
@@ -484,6 +485,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat
     end
 
     @trixi_timeit timer() "solving elliptic system" begin
+        scale_by_mass_matrix!(tmp1, D1)
         solve_system_matrix!(dv, system_matrix, tmp,
                              equations, D1, cache)
     end
@@ -589,6 +591,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache,
     end
 
     @trixi_timeit timer() "solving elliptic system" begin
+        scale_by_mass_matrix!(tmp1, D1)
         solve_system_matrix!(dv, system_matrix, tmp,
                              equations, D1, cache)
     end
@@ -702,6 +705,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache,
     end
 
     @trixi_timeit timer() "solving elliptic system" begin
+        scale_by_mass_matrix!(tmp1, D1)
         solve_system_matrix!(dv, system_matrix, tmp,
                              equations, D1, cache)
     end

diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl
@@ -442,6 +442,7 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D,
     end
     @trixi_timeit timer() "solving elliptic system" begin
         tmp1 .= dv
+        scale_by_mass_matrix!(tmp1, D1)
         solve_system_matrix!(dv, system_matrix, tmp1,
                              equations, D1, cache)
     end
@@ -487,8 +488,9 @@ function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D,
     end
     @trixi_timeit timer() "solving elliptic system" begin
         tmp1 .= dv
+        scale_by_mass_matrix!(tmp1, D1)
         solve_system_matrix!((@view dv[2:(end - 1)]), system_matrix,
-                             tmp1, equations, D1,
+                             (@view tmp1[2:(end - 1)]), equations, D1,
                              cache)
         dv[1] = dv[end] = zero(eltype(dv))
     end

diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl
@@ -24,11 +24,11 @@ end
 ] begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"),
                         tspan=(0.0, 1.0),
-                        l2=[5.402315467757905e-6 6.705583190754449e-8 0.0],
-                        linf=[9.787584609100008e-5 1.4602668577112787e-7 0.0],
-                        cons_error=[1.0484865187415942e-9 0.5469460930247753 0.0],
-                        change_waterheight=1.0484865187415942e-9,
-                        change_entropy_modified=459.9037236233692,
+                        l2=[5.402315494643131e-6 6.70558305594986e-8 0.0],
+                        linf=[9.787585531206844e-5 1.460266869784954e-7 0.0],
+                        cons_error=[1.0484871583691058e-9 0.5469460930247998 0.0],
+                        change_waterheight=1.0484871583691058e-9,
+                        change_entropy_modified=459.90372362340514,
                         atol=1e-11) # to make CI pass
 
     @test_allocations(semi, sol, allocs=650_000)
@@ -42,11 +42,11 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_basic_reflecting.jl"),
                         tspan=(0.0, 1.0),
                         solver=solver,
-                        l2=[5.862278093002509e-6 4.111933117101271e-9 0.0],
-                        linf=[3.135229084794133e-5 8.797788287467911e-8 0.0],
-                        cons_error=[1.700425282207037e-9 0.5469460935004404 0.0],
-                        change_waterheight=-1.700425282207037e-9,
-                        change_entropy_modified=459.9037221442477,
+                        l2=[5.862278175937948e-6 4.11195454078554e-9 0.0],
+                        linf=[3.135228725170691e-5 8.797787950237668e-8 0.0],
+                        cons_error=[1.700425028441774e-9 0.5469460935005555 0.0],
+                        change_waterheight=-1.700425028441774e-9,
+                        change_entropy_modified=459.9037221442321,
                         atol=1e-11) # to make CI pass
 
     @test_allocations(semi, sol, allocs=650_000)