Solvers

This chapter covers different solvers and how to use them in DispersiveShallowWater.jl.

To learn more about the analytical and mathematical background, go to the chapter about Summation-by-Parts Operators.

Introduction

DispersiveShallowWater.jl uses SummationByPartsOperators.jl as the foundation for all spatial discretizations. This package provides a comprehensive collection of summation-by-parts (SBP) operators that enable structure-preserving numerical methods with provable stability and conservation properties.

The Solver struct in DispersiveShallowWater.jl wraps the SBP operators needed for spatial discretization: typically a first-derivative operator D1, and optionally second- and third-derivative operators D2 and D3, depending on the equation system being solved.

Design Philosophy: Flexibility and Modularity

DispersiveShallowWater.jl follows a modular design that separates the choice of spatial discretization from the underlying equation systems. This design allows users to:

• Compare discretization methods: Test finite difference, discontinuous Galerkin, continuous Galerkin, and Fourier spectral approaches on identical problems
• Investigate operator properties: Analyze how different SBP operators preserve conservation laws and stability properties

This modularity is particularly relevant for research applications where the choice of discretization can significantly impact the preservation of physical invariants and long-time numerical stability. The unified interface enables systematic studies of numerical method performance without requiring separate implementations for each approach.

Basic Summation-by-Parts Solver

As we have already seen in the basic example, the easiest way to create a solver is to pass both a mesh and the desired accuracy order to the Solver function. This creates first-, second-, and third-derivative periodic summation-by-parts operators of the given accuracy order on the specified mesh:

accuracy_order = 4solver = Solver(mesh, accuracy_order)

This approach creates periodic SBP operators and only works with periodic boundary conditions. If a solver containing periodic SBP operators is passed to Semidiscretization with non-periodic boundary conditions, it will throw an error.

Customizing Solvers

For more advanced use cases, you can create custom solvers by explicitly constructing the required SBP operators. The general pattern is:

solver = Solver(D1, D2, D3)

where:

• D1 is always required and must be an AbstractDerivativeOperator
• D2 and D3 are optional and can be either AbstractDerivativeOperators, AbstractMatrixes, or nothing

Reflecting Boundary Conditions

For non-periodic boundary conditions, you need to use non-periodic SBP operators. Here's how to create a solver for reflecting boundary conditions:

using SummationByPartsOperators: MattssonNordström2004, derivative_operator# Create solver with SBP operators of accuracy order 4accuracy_order = 4D1 = derivative_operator(MattssonNordström2004(),                         derivative_order = 1, accuracy_order = accuracy_order,                         xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N)solver = Solver(D1)

The MattssonNordström2004() operator family provides SBP operators that satisfy the SBP property with non-periodic boundary conditions, making them suitable for reflecting boundary conditions. Other possible choices for sources can be found in the documentation of SummationByPartsOperators.jl.

Upwind Operators

For equations that benefit from upwind discretizations (such as the Serre-Green-Naghdi equations), you can use upwind SBP Operators:

For periodic boundary conditions:

using SummationByPartsOperators: upwind_operators, periodic_derivative_operatoraccuracy_order = 4D1 = upwind_operators(periodic_derivative_operator;                      derivative_order = 1,                      accuracy_order = accuracy_order,                      xmin = xmin(mesh), xmax = xmax(mesh),                      N = nnodes(mesh))solver = Solver(D1)

For reflecting boundary conditions:

using SummationByPartsOperators: Mattsson2017, upwind_operatorsaccuracy_order = 2D1 = upwind_operators(Mattsson2017, derivative_order = 1, accuracy_order = accuracy_order,                      xmin = xmin(mesh), xmax = xmax(mesh),                      N = nnodes(mesh))solver = Solver(D1)

Upwind operators provide additional dissipation and often achieve higher resolution for sharp features compared to central operators.

KdV Equations with Third-Derivative Operators

The KdV equation requires first- and third-derivative operators, but no second-derivative operator. You can specify this by passing nothing for the second argument:

using SummationByPartsOperators: fourier_derivative_operator# Create solver with Fourier SBP operatorsD1 = fourier_derivative_operator(xmin(mesh), xmax(mesh), nnodes(mesh))D3 = D1^3  # Third derivative via compositionsolver = Solver(D1, nothing, D3)

Using Sparse Matrices for Derivative Operators

While D1 must always be an SBP operator, D2 and D3 can be regular (sparse) matrices. This can be useful for creating custom discretizations:

using SummationByPartsOperators: upwind_operators, periodic_derivative_operatorusing SparseArrays: sparseaccuracy_order = 4D1 = upwind_operators(periodic_derivative_operator; derivative_order = 1,                      accuracy_order = accuracy_order, xmin = mesh.xmin, xmax = mesh.xmax,                      N = mesh.N)D2 = sparse(D1.plus) * sparse(D1.minus)  # Create D2 as sparse matrixsolver = Solver(D1, D2)

Discontinuous and continuous Galerkin Methods

DispersiveShallowWater.jl supports discontinuous Galerkin (DG) and continuous Galerkin (CG) methods through coupled Lobatto-Legendre operators:

using SummationByPartsOperators: legendre_derivative_operator,                                 UniformPeriodicMesh1D,                                 couple_discontinuouslyusing SparseArrays: sparse# Create solver with DG operatorsp = 3  # polynomial degree (N needs to be divisible by p + 1)D_legendre = legendre_derivative_operator(-1.0, 1.0, p + 1)uniform_mesh = UniformPeriodicMesh1D(coordinates_min, coordinates_max, div(N, p + 1))D1 = couple_discontinuously(D_legendre, uniform_mesh)D_pl = couple_discontinuously(D_legendre, uniform_mesh, Val(:plus))D_min = couple_discontinuously(D_legendre, uniform_mesh, Val(:minus))D2 = sparse(D_pl) * sparse(D_min)solver = Solver(D1, D2)

This approach creates a DG discretization with polynomial degree p within each element. Similarly, you can create a CG discretization using couple_continuously. Note that in this case, N needs to be divisible by p instead of p + 1.

Fourier Spectral Methods

Fourier collocation methods can be interpreted as periodic SBP operators, which can be constructed via fourier_derivative_operator:

using SummationByPartsOperators: fourier_derivative_operator# Create solver with Fourier pseudospectral collocation methodD1 = fourier_derivative_operator(xmin(mesh), xmax(mesh), nnodes(mesh))solver = Solver(D1)

Variable Coefficient Operators

Variable coefficient operators can be useful when solving PDEs that involve elliptic systems with variable coefficients. Using them can lead to fewer allocations.

using SummationByPartsOperators: Mattsson2012, derivative_operator,                                 var_coef_derivative_operatoraccuracy_order = 2D1 = derivative_operator(Mattsson2012();                         derivative_order = 1, accuracy_order,                         xmin = xmin(mesh), xmax = xmax(mesh), N = N)# Create a variable-coefficient second-derivative operator# Initialize with coefficient function `one` - coefficients will be set during simulationD2 = var_coef_derivative_operator(Mattsson2012(),                                  2, accuracy_order,                                  xmin(mesh), xmax(mesh), N, one)solver = Solver(D1, D2)

Using Your Custom Solver

Once you have created your solver object using any of the approaches above, you can use it in a Semidiscretization:

semi = Semidiscretization(mesh, equations, initial_condition, solver,                          boundary_conditions = boundary_conditions)

The solver you choose should be compatible with your boundary conditions:

• Periodic operators (created with periodic_derivative_operator or fourier_derivative_operator) require boundary_condition_periodic
• Non-periodic operators (created with MattssonNordström2004, etc.) are needed for boundary_condition_reflecting

Common Pitfalls

Mismatched Operators and Boundary Conditions

Problem: Using periodic operators with reflecting boundary conditions or vice versa.

# ❌ This will failD1 = periodic_derivative_operator(1, 4, mesh.xmin, mesh.xmax, mesh.N)solver = Solver(D1)semi = Semidiscretization(mesh, equations, initial_condition, solver,                          boundary_conditions = boundary_condition_reflecting)  # Error!

Solution: Match operator type to boundary conditions:

# ✅ Correct approachD1 = derivative_operator(MattssonNordström2004(), derivative_order = 1,                         accuracy_order = 4, xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N)solver = Solver(D1)semi = Semidiscretization(mesh, equations, initial_condition, solver,                          boundary_conditions = boundary_condition_reflecting)

Incorrect Grid Size for DG Methods

Problem: Grid size not divisible by polynomial degree + 1 for DG methods.

# ❌ N = 101, p = 3, but 101 is not divisible by (3+1) = 4N = 101p = 3mesh = Mesh1D(coordinates_min, coordinates_max, N)  # Error in DG setup!

Solution: Ensure N is divisible by p + 1:

# ✅ Adjust N to be divisible by p + 1p = 3N = 100  # 100 ÷ 4 = 25 elements exactly# Or N = 96, 104, etc. - any multiple of 4