Merge pull request #1275 from trixi-framework/dev

Upwind finite difference SBP solver
trixi-framework · Dec 2, 2022 · d9b5c3c · d9b5c3c
2 parents 94703d9 + 9498920
commit d9b5c3c
Show file tree

Hide file tree

Showing 34 changed files with 2,766 additions and 79 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -9,6 +9,9 @@ for human readability.
 
 #### Added
 
+- Experimental support for upwind finite difference summation by parts (FDSBP)
+  has been added. The first implementation requires a `TreeMesh` and comes
+  with several examples in the `examples_dir()` of Trixi.jl.
 - Experimental support for 2D parabolic diffusion terms has been added.
   * `LaplaceDiffusion2D` and `CompressibleNavierStokesDiffusion2D` can be used to add
   diffusion to systems. `LaplaceDiffusion2D` can be used to add scalar diffusion to each

diff --git a/Project.toml b/Project.toml
@@ -67,7 +67,7 @@ Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8"
 StaticArrays = "1"
 StrideArrays = "0.1.18"
 StructArrays = "0.6"
-SummationByPartsOperators = "0.5.10"
+SummationByPartsOperators = "0.5.25"
 TimerOutputs = "0.5"
 Triangulate = "2.0"
 TriplotBase = "0.1"

diff --git a/docs/literate/src/files/index.jl b/docs/literate/src/files/index.jl
@@ -56,33 +56,41 @@
 # For instance, we show how to set up a finite differences (FD) scheme and a continuous Galerkin
 # (CGSEM) method.
 
-# ### [7 Adding a new scalar conservation law](@ref adding_new_scalar_equations)
+# ### [7 Upwind FD SBP schemes](@ref upwind_fdsbp)
+#-
+# General SBP schemes can not only be used via the [`DGMulti`](@ref) solver but
+# also with a general `DG` solver. In particular, upwind finite difference SBP
+# methods can be used together with the `TreeMesh`. Similar to general SBP
+# schemes in the `DGMulti` framework, the interface is based on the package
+# [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl).
+
+# ### [8 Adding a new scalar conservation law](@ref adding_new_scalar_equations)
 #-
 # This tutorial explains how to add a new physics model using the example of the cubic conservation
 # law. First, we define the equation using a `struct` `CubicEquation` and the physical flux. Then,
 # the corresponding standard setup in Trixi.jl (`mesh`, `solver`, `semi` and `ode`) is implemented
 # and the ODE problem is solved by OrdinaryDiffEq's `solve` method.
 
-# ### [8 Adding a non-conservative equation](@ref adding_nonconservative_equation)
+# ### [9 Adding a non-conservative equation](@ref adding_nonconservative_equation)
 #-
 # In this part, another physics model is implemented, the nonconservative linear advection equation.
 # We run two different simulations with different levels of refinement and compare the resulting errors.
 
-# ### [9 Adaptive mesh refinement](@ref adaptive_mesh_refinement)
+# ### [10 Adaptive mesh refinement](@ref adaptive_mesh_refinement)
 #-
 # Adaptive mesh refinement (AMR) helps to increase the accuracy in sensitive or turbolent regions while
 # not wasting ressources for less interesting parts of the domain. This leads to much more efficient
 # simulations. This tutorial presents the implementation strategy of AMR in Trixi, including the use of
 # different indicators and controllers.
 
-# ### [10 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping)
+# ### [11 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping)
 #-
 # In this tutorial, the use of Trixi's structured curved mesh type [`StructuredMesh`](@ref) is explained.
 # We present the two basic option to initialize such a mesh. First, the curved domain boundaries
 # of a circular cylinder are set by explicit boundary functions. Then, a fully curved mesh is
 # defined by passing the transformation mapping.
 
-# ### [11 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial)
+# ### [12 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial)
 #-
 # The purpose of this tutorial is to demonstrate how to use the [`UnstructuredMesh2D`](@ref)
 # functionality of Trixi.jl. This begins by running and visualizing an available unstructured
@@ -91,13 +99,13 @@
 # software in the Trixi.jl ecosystem, and then run a simulation using Trixi.jl on said mesh.
 # In the end, the tutorial briefly explains how to simulate an example using AMR via `P4estMesh`.
 
-# ### [12 Explicit time stepping](@ref time_stepping)
+# ### [13 Explicit time stepping](@ref time_stepping)
 # -
 # This tutorial is about time integration using [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl).
 # It explains how to use their algorithms and presents two types of time step choices - with error-based
 # and CFL-based adaptive step size control.
 
-# ### [13 Differentiable programming](@ref differentiable_programming)
+# ### [14 Differentiable programming](@ref differentiable_programming)
 #-
 # This part deals with some basic differentiable programming topics. For example, a Jacobian, its
 # eigenvalues and a curve of total energy (through the simulation) are calculated and plotted for

diff --git a/docs/literate/src/files/upwind_fdsbp.jl b/docs/literate/src/files/upwind_fdsbp.jl
@@ -0,0 +1,64 @@
+#src # Upwind FD SBP schemes
+
+# General tensor product SBP methods are supported via the `DGMulti` solver
+# in a reasonably complete way, see [the previous tutorial](@ref DGMulti_2).
+# Nevertheless, there is also experimental support for SBP methods with
+# other solver and mesh types.
+
+# The first step is to set up an SBP operator. A classical (central) SBP
+# operator can be created as follows.
+using Trixi
+D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
+                            derivative_order=1, accuracy_order=2,
+                            xmin=0.0, xmax=1.0, N=11)
+# Instead of prefixing the source of coefficients `MattssonNordström2004()`,
+# you can also load the package SummationByPartsOperators.jl. Either way,
+# this yields an object representing the operator efficiently. If you want to
+# compare it to coefficients presented in the literature, you can convert it
+# to a matrix.
+Matrix(D_SBP)
+
+# Upwind SBP operators are a concept introduced in 2017 by Ken Mattsson. You can
+# create them as follows.
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order=1, accuracy_order=2,
+                         xmin=0.0, xmax=1.0, N=11)
+# Upwind operators are derivative operators biased towards one direction.
+# The "minus" variants has a bias towards the left side, i.e., it uses values
+# from more nodes to the left than from the right to compute the discrete
+# derivative approximation at a given node (in the interior of the domain).
+# In matrix form, this means more non-zero entries are left from the diagonal.
+Matrix(D_upw.minus)
+# Analogously, the "plus" variant has a bias towards the right side.
+Matrix(D_upw.plus)
+# For more information on upwind SBP operators, please refer to the documentation
+# of [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl)
+# and references cited there.
+
+# The basic idea of upwind SBP schemes is to apply a flux vector splitting and
+# use appropriate upwind operators for both parts of the flux. In 1D, this means
+# to split the flux
+# ```math
+# f(u) = f^-(u) + f^+(u)
+# ```
+# such that $f^-(u)$ is associated with left-going waves and $f^+(u)$ with
+# right-going waves. Then, we apply upwind SBP operators $D^-, D^+$ with an
+# appropriate upwind bias, resulting in
+# ```math
+# \partial_x f(u) \approx D^+ f^-(u) + D^- f^+(u)
+# ```
+# Note that the established notations of upwind operators $D^\pm$ and flux
+# splittings $f^\pm$ clash. The right-going waves from $f^+$ need an operator
+# biased towards their upwind side, i.e., the left side. This upwind bias is
+# provided by the operator $D^-$.
+
+# Many classical flux vector splittings have been developed for finite volume
+# methods and are described in the book "Riemann Solvers and Numerical Methods
+# for Fluid Dynamics: A Practical Introduction" of Eleuterio F. Toro (2009),
+# [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761). One such a well-known
+# splitting provided by Trixi.jl is [`splitting_steger_warming`](@ref).
+
+# Trixi.jl comes with several example setups using upwind SBP methods with
+# flux vector splitting, e.g.,
+# - [`elixir_euler_vortex.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_fdsbp/elixir_euler_vortex.jl)
+# - [`elixir_euler_taylor_green_vortex.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_3d_fdsbp/elixir_euler_taylor_green_vortex.jl)
diff --git a/docs/make.jl b/docs/make.jl
@@ -36,6 +36,7 @@ files = [
     "Non-periodic boundaries" => "non_periodic_boundaries.jl",
     "DG schemes via `DGMulti` solver" => "DGMulti_1.jl",
     "Other SBP schemes (FD, CGSEM) via `DGMulti` solver" => "DGMulti_2.jl",
+    "Upwind FD SBP schemes" => "upwind_fdsbp.jl",
     # Topic: equations
     "Adding a new scalar conservation law" => "adding_new_scalar_equations.jl",
     "Adding a non-conservative equation" => "adding_nonconservative_equation.jl",

diff --git a/examples/tree_1d_fdsbp/elixir_burgers_basic.jl b/examples/tree_1d_fdsbp/elixir_burgers_basic.jl
@@ -0,0 +1,64 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the (inviscid) Burgers' equation
+
+equations = InviscidBurgersEquation1D()
+
+initial_condition = initial_condition_convergence_test
+
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order=1,
+                         accuracy_order=4,
+                         xmin=-1.0, xmax=1.0,
+                         N=32)
+flux_splitting = splitting_lax_friedrichs
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
+
+coordinates_min = 0.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=3,
+                n_cells_max=10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms=source_terms_convergence_test)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 200
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_errors=(:l2_error_primitive,
+                                                            :linf_error_primitive))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(), abstol=1.0e-9, reltol=1.0e-9,
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/tree_1d_fdsbp/elixir_burgers_linear_stability.jl b/examples/tree_1d_fdsbp/elixir_burgers_linear_stability.jl
@@ -0,0 +1,60 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the (inviscid) Burgers' equation
+
+equations = InviscidBurgersEquation1D()
+
+function initial_condition_linear_stability(x, t, equation::InviscidBurgersEquation1D)
+  k = 1
+  2 + sinpi(k * (x[1] - 0.7)) |> SVector
+end
+
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order=1,
+                         accuracy_order=4,
+                         xmin=-1.0, xmax=1.0,
+                         N=16)
+flux_splitting = splitting_lax_friedrichs
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
+
+coordinates_min = -1.0
+coordinates_max =  1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_linear_stability, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_errors=(:l2_error_primitive,
+                                                            :linf_error_primitive))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(), abstol=1.0e-6, reltol=1.0e-6,
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/tree_1d_fdsbp/elixir_euler_convergence.jl b/examples/tree_1d_fdsbp/elixir_euler_convergence.jl
@@ -0,0 +1,64 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                          derivative_order=1,
+                          accuracy_order=4,
+                          xmin=-1.0, xmax=1.0,
+                          N=32)
+flux_splitting = splitting_steger_warming
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralUpwind(flux_splitting),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
+
+coordinates_min = 0.0
+coordinates_max = 2.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=1,
+                n_cells_max=10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms=source_terms_convergence_test)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_errors=(:l2_error_primitive,
+                                                            :linf_error_primitive))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(), abstol=1.0e-6, reltol=1.0e-6,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary