Factorize matrix for SvaerdKalischEquations1D

src/DispersiveShallowWater.jl

@@ -3,7 +3,7 @@ module DispersiveShallowWater
 using BandedMatrices: BandedMatrix
 using DiffEqBase: DiffEqBase, SciMLBase, terminate!
 using Interpolations: Interpolations, linear_interpolation
-using LinearAlgebra: mul!, ldiv!, I, Diagonal, diag, lu
+using LinearAlgebra: mul!, ldiv!, I, Diagonal, Symmetric, diag, lu, cholesky, cholesky!
 using PolynomialBases: PolynomialBases
 using Printf: @printf, @sprintf
 using RecipesBase: RecipesBase, @recipe, @series

src/equations/svaerd_kalisch_1d.jl

@@ -177,7 +177,6 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D,
     beta_hat = equations.beta * D .^ 3
     gamma_hat = equations.gamma * sqrt.(equations.gravity * D) .* D .^ 3
     tmp2 = similar(h)
-    hmD1betaD1 = Array{RealT}(undef, nnodes(mesh), nnodes(mesh))
     if solver.D1 isa PeriodicDerivativeOperator ||
        solver.D1 isa UniformPeriodicCoupledOperator
         D1_central = solver.D1
@@ -189,9 +188,11 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D,
     else
         @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
     end
-    return (hmD1betaD1 = hmD1betaD1, D1betaD1 = D1betaD1, D = D, h = h, hv = hv,
-            alpha_hat = alpha_hat, beta_hat = beta_hat, gamma_hat = gamma_hat,
-            tmp2 = tmp2, D1_central = D1_central, D1 = solver.D1)
+    M = mass_matrix(solver.D1)
+    factorization = cholesky(Symmetric(M * (Diagonal(ones(nnodes(mesh))) - D1betaD1)))
+    return (factorization = factorization, D1betaD1 = D1betaD1, D = D, h = h, hv = hv,
+            alpha_hat = alpha_hat, gamma_hat = gamma_hat,
+            tmp2 = tmp2, D1_central = D1_central, M = M, D1 = solver.D1)
 end
 
 # Discretization that conserves the mass (for eta and for flat bottom hv) and the energy for periodic boundary conditions, see
@@ -201,7 +202,7 @@ end
 function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D,
               initial_condition, ::BoundaryConditionPeriodic, source_terms,
               solver, cache)
-    @unpack hmD1betaD1, D1betaD1, D, h, hv, alpha_hat, beta_hat, gamma_hat, tmp1, tmp2, D1_central = cache
+    @unpack factorization, D1betaD1, D, h, hv, alpha_hat, gamma_hat, tmp1, tmp2, D1_central, M = cache
     q = wrap_array(u_ode, mesh, equations, solver)
     dq = wrap_array(du_ode, mesh, equations, solver)
 
@@ -213,8 +214,8 @@ function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D,
     fill!(dD, zero(eltype(dD)))
 
     @timeit timer() "deta hyperbolic" begin
-        h = eta .+ D .- equations.eta0
-        hv = h .* v
+        @. h = eta + D - equations.eta0
+        @. hv = h * v
 
         if solver.D1 isa PeriodicDerivativeOperator ||
            solver.D1 isa UniformPeriodicCoupledOperator
@@ -241,11 +242,15 @@ function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D,
 
     # split form
     @timeit timer() "dv hyperbolic" begin
-        dv[:] = -(0.5 * (D1_central * (hv .* v) + hv .* D1v - v .* (D1_central * hv)) +
-                  equations.gravity * h .* D1eta +
+        D1_hv = D1_central * hv
+        D1_hv2 = D1_central * (hv .* v)
+        D1_gamma_hat_D2_v = D1_central * (gamma_hat .* (solver.D2 * v))
+        D2_gamma_hat_D1_v = solver.D2 * (gamma_hat .* D1v)
+        @. dv = -(0.5 * (D1_hv2 + hv * D1v - v * D1_hv) +
+                  equations.gravity * h * D1eta +
                   0.5 * (vD1y - D1vy - yD1v) -
-                  0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) -
-                  0.5 * solver.D2 * (gamma_hat .* D1v))
+                  0.5 * D1_gamma_hat_D2_v -
+                  0.5 * D2_gamma_hat_D1_v)
     end
 
     # no split form
@@ -257,8 +262,12 @@ function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D,
 
     @timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, solver)
     @timeit timer() "dv elliptic" begin
-        hmD1betaD1 = Diagonal(h) - D1betaD1
-        dv[:] = hmD1betaD1 \ dv
+        # decompose M * (h - D1betaD1) because it is guaranteed to be symmetric and pos. def.,
+        # while (h - D1betaD1) is not necessarily
+        hmD1betaD1 = Symmetric(M * (Diagonal(h) - D1betaD1))
+        cholesky!(factorization, hmD1betaD1)
+        tmp1[:] = M * dv
+        dv[:] = factorization \ tmp1
     end
 
     return nothing

test/test_visualization.jl

@@ -4,18 +4,16 @@ using Test
 using DispersiveShallowWater
 using Plots
 
-@testset "Unit tests" begin
-    @testset "Visualization" begin
-        trixi_include(@__MODULE__, default_example(), tspan = (0.0, 1.0))
-        @test_nowarn plot(semi => sol)
-        @test_nowarn plot!(semi => sol, plot_initial = true)
-        @test_nowarn plot(semi, sol, conversion = prim2cons, plot_bathymetry = false)
-        @test_nowarn plot(semi => sol, 0.0)
-        @test_nowarn plot(semi, sol, 0.0, conversion = prim2cons)
-        @test_nowarn plot(analysis_callback)
-        @test_nowarn plot(analysis_callback, what = (:errors,))
-        @test_nowarn plot(analysis_callback, what = (:integrals, :errors))
-    end
+@testset "Visualization" begin
+    trixi_include(@__MODULE__, default_example(), tspan = (0.0, 1.0))
+    @test_nowarn plot(semi => sol)
+    @test_nowarn plot!(semi => sol, plot_initial = true)
+    @test_nowarn plot(semi, sol, conversion = prim2cons, plot_bathymetry = false)
+    @test_nowarn plot(semi => sol, 0.0)
+    @test_nowarn plot(semi, sol, 0.0, conversion = prim2cons)
+    @test_nowarn plot(analysis_callback)
+    @test_nowarn plot(analysis_callback, what = (:errors,))
+    @test_nowarn plot(analysis_callback, what = (:integrals, :errors))
 end
 
 end # module