Add reflecting boundary conditions for Svärd-Kalisch equations (#166)

* add reflecting boundary conditions for Svärd-Kalisch equations

* Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* add back accidently removed FourierDerivativeOperator

* also implement upwind operators

* reformat

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/equations/svaerd_kalisch_1d.jl

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* set stricter time integration tolerances

* lower tolerance

* use Cholesky decomposition

* Update src/equations/svaerd_kalisch_1d.jl

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* temporarily add lu decomposition again

* fix scale_by_mass_matrix!

* remove LU decomposition again

* scale_by_mass_matrix! outside of solve_system_matrix!

* try CI without tolerances

* reorder imports

* fix variable name

* refactor to reduce code redundancy

* fix

* pass boundary_conditions

* format

* adjust tolerances

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/equations/equations.jl

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* begin + 1instead of 2

* set tolerances in tests

* add test for ArgumentError

---------

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: Hendrik Ranocha <[email protected]>

examples/svaerd_kalisch_1d/svaerd_kalisch_1d_basic_reflecting.jl

@@ -0,0 +1,46 @@
+using OrdinaryDiffEqTsit5
+using DispersiveShallowWater
+using SummationByPartsOperators: MattssonNordström2004, derivative_operator
+
+###############################################################################
+# Semidiscretization of the Svärd-Kalisch equations
+# For reflecting boundary conditions, alpha and gamma need to be 0
+equations = SvaerdKalischEquations1D(gravity_constant = 1.0, eta0 = 0.0,
+                                     alpha = 0.0, beta = 1 / 3, gamma = 0.0)
+
+initial_condition = initial_condition_manufactured_reflecting
+source_terms = source_terms_manufactured_reflecting
+boundary_conditions = boundary_condition_reflecting
+
+# create homogeneous mesh
+coordinates_min = 0.0
+coordinates_max = 1.0
+N = 512
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with SBP operators of accuracy order 4
+accuracy_order = 4
+D1 = derivative_operator(MattssonNordström2004(),
+                         derivative_order = 1, accuracy_order = accuracy_order,
+                         xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N)
+solver = Solver(D1, nothing)
+
+# semidiscretization holds all the necessary data structures for the spatial discretization
+semi = Semidiscretization(mesh, equations, initial_condition, solver,
+                          boundary_conditions = boundary_conditions,
+                          source_terms = source_terms)
+
+###############################################################################
+# Create `ODEProblem` and run the simulation
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 100,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 entropy_modified))
+callbacks = CallbackSet(analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 100)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+            save_everystep = false, callback = callbacks, saveat = saveat)

src/equations/bbm_bbm_1d.jl

@@ -338,8 +338,7 @@ function create_cache(mesh, equations::BBMBBMEquations1D{BathymetryFlat},
 
     # homogeneous Neumann boundary conditions
     if solver.D1 isa DerivativeOperator ||
-       solver.D1 isa UniformCoupledOperator ||
-       solver.D1 isa FourierDerivativeOperator
+       solver.D1 isa UniformCoupledOperator
         D1_b = BandedMatrix(solver.D1)
         invImD2n = lu(I + 1 / 6 * D^2 * inv(M) * D1_b' * PdM * D1_b)
     elseif solver.D1 isa UpwindOperators
@@ -382,8 +381,7 @@ function create_cache(mesh, equations::BBMBBMEquations1D,
 
     # homogeneous Neumann boundary conditions
     if solver.D1 isa DerivativeOperator ||
-       solver.D1 isa UniformCoupledOperator ||
-       solver.D1 isa FourierDerivativeOperator
+       solver.D1 isa UniformCoupledOperator
         D1_b = BandedMatrix(solver.D1)
         invImD2n = lu(I + 1 / 6 * inv(M) * D1_b' * PdM * K * D1_b)
     elseif solver.D1 isa UpwindOperators
@@ -501,8 +499,7 @@ function rhs!(dq, q, t, mesh, equations::BBMBBMEquations1D, initial_condition,
     end
     # energy and mass conservative semidiscretization
     if solver.D1 isa DerivativeOperator ||
-       solver.D1 isa UniformCoupledOperator ||
-       solver.D1 isa FourierDerivativeOperator
+       solver.D1 isa UniformCoupledOperator
         @trixi_timeit timer() "deta hyperbolic" begin
             mul!(deta, solver.D1, tmp1)
         end

src/equations/equations.jl

@@ -557,6 +557,23 @@ abstract type AbstractSerreGreenNaghdiEquations{NDIMS, NVARS} <:
 include("serre_green_naghdi_1d.jl")
 include("hyperbolic_serre_green_naghdi_1d.jl")
 
+function solve_system_matrix!(dv, system_matrix, rhs,
+                              equations::Union{SvaerdKalischEquations1D,
+                                               SerreGreenNaghdiEquations1D},
+                              D1, cache, ::BoundaryConditionPeriodic)
+    scale_by_mass_matrix!(rhs, D1)
+    solve_system_matrix!(dv, system_matrix, rhs, equations, D1, cache)
+end
+
+function solve_system_matrix!(dv, system_matrix, rhs,
+                              equations::Union{SvaerdKalischEquations1D,
+                                               SerreGreenNaghdiEquations1D},
+                              D1, cache, ::BoundaryConditionReflecting)
+    scale_by_mass_matrix!(rhs, D1)
+    solve_system_matrix!(dv, system_matrix, (@view rhs[(begin + 1):(end - 1)]), equations,
+                         D1, cache)
+end
+
 function solve_system_matrix!(dv, system_matrix, rhs,
                               ::Union{SvaerdKalischEquations1D,
                                       SerreGreenNaghdiEquations1D},
@@ -565,7 +582,6 @@ 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
         else
@@ -575,7 +591,6 @@ function solve_system_matrix!(dv, system_matrix, rhs,
         end
     else
         factorization = cholesky!(system_matrix)
-        scale_by_mass_matrix!(rhs, D1)
         ldiv!(dv, factorization, rhs)
     end
     return nothing

src/equations/serre_green_naghdi_1d.jl

@@ -314,19 +314,22 @@ end
 function rhs!(dq, q, t, mesh,
               equations::SerreGreenNaghdiEquations1D,
               initial_condition,
-              ::BoundaryConditionPeriodic,
+              boundary_conditions::BoundaryConditionPeriodic,
               source_terms::Nothing,
               solver, cache)
     if cache.D1 isa PeriodicUpwindOperators
-        rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry_type)
+        rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry_type,
+                        boundary_conditions)
     else
-        rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry_type)
+        rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry_type,
+                         boundary_conditions)
     end
 
     return nothing
 end
 
-function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFlat)
+function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFlat,
+                          boundary_conditions::BoundaryConditionPeriodic)
     # Unpack physical parameters and SBP operator `D1` as well as the
     # SBP operator in sparse matrix form `D1mat`
     g = gravity_constant(equations)
@@ -399,14 +402,15 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla
     end
 
     @trixi_timeit timer() "solving elliptic system" begin
-        solve_system_matrix!(dv, system_matrix, tmp,
-                             equations, D1, cache)
+        solve_system_matrix!(dv, system_matrix,
+                             tmp, equations, D1, cache, boundary_conditions)
     end
 
     return nothing
 end
 
-function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat)
+function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat,
+                         boundary_conditions::BoundaryConditionPeriodic)
     # Unpack physical parameters and SBP operator `D1` as well as the
     # SBP upwind operator in sparse matrix form `D1mat_minus`
     g = gravity_constant(equations)
@@ -484,15 +488,16 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat
     end
 
     @trixi_timeit timer() "solving elliptic system" begin
-        solve_system_matrix!(dv, system_matrix, tmp,
-                             equations, D1, cache)
+        solve_system_matrix!(dv, system_matrix,
+                             tmp, equations, D1, cache, boundary_conditions)
     end
 
     return nothing
 end
 
 function rhs_sgn_central!(dq, q, equations, source_terms, cache,
-                          ::Union{BathymetryMildSlope, BathymetryVariable})
+                          ::Union{BathymetryMildSlope, BathymetryVariable},
+                          boundary_conditions::BoundaryConditionPeriodic)
     # Unpack physical parameters and SBP operator `D1` as well as the
     # SBP operator in sparse matrix form `D1mat`
     g = gravity_constant(equations)
@@ -589,15 +594,16 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache,
     end
 
     @trixi_timeit timer() "solving elliptic system" begin
-        solve_system_matrix!(dv, system_matrix, tmp,
-                             equations, D1, cache)
+        solve_system_matrix!(dv, system_matrix,
+                             tmp, equations, D1, cache, boundary_conditions)
     end
 
     return nothing
 end
 
 function rhs_sgn_upwind!(dq, q, equations, source_terms, cache,
-                         ::Union{BathymetryMildSlope, BathymetryVariable})
+                         ::Union{BathymetryMildSlope, BathymetryVariable},
+                         boundary_conditions::BoundaryConditionPeriodic)
     # Unpack physical parameters and SBP operator `D1` as well as the
     # SBP operator in sparse matrix form `D1mat`
     g = gravity_constant(equations)
@@ -702,8 +708,8 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache,
     end
 
     @trixi_timeit timer() "solving elliptic system" begin
-        solve_system_matrix!(dv, system_matrix, tmp,
-                             equations, D1, cache)
+        solve_system_matrix!(dv, system_matrix,
+                             tmp, equations, D1, cache, boundary_conditions)
     end
 
     return nothing