From 80237db76f0565223c7a0d063f74bb2247e42af6 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 30 Nov 2022 08:01:21 +0100 Subject: [PATCH 1/8] add NEWS --- NEWS.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/NEWS.md b/NEWS.md index 1cddb754d07..fbe6abc48c8 100644 --- 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 From 501466f5b76cf15ce019ec936ad01cdd3876f648 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 30 Nov 2022 08:59:30 +0100 Subject: [PATCH 2/8] WIP: tutorial on upwind FD SBP methods --- docs/literate/src/files/index.jl | 22 +++++++++++------ docs/literate/src/files/upwind_fdsbp.jl | 32 +++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 7 deletions(-) create mode 100644 docs/literate/src/files/upwind_fdsbp.jl diff --git a/docs/literate/src/files/index.jl b/docs/literate/src/files/index.jl index f4a8d83b845..899c18ff8d7 100644 --- 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 new file mode 100644 index 00000000000..d103c4ddc70 --- /dev/null +++ b/docs/literate/src/files/upwind_fdsbp.jl @@ -0,0 +1,32 @@ +#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. +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. +# 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. +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. From 47034469f1eb4b9c23dc00cb3854d0b4c07377f4 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 30 Nov 2022 09:32:29 +0100 Subject: [PATCH 3/8] first draft of tutorial --- docs/literate/src/files/upwind_fdsbp.jl | 28 +++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/docs/literate/src/files/upwind_fdsbp.jl b/docs/literate/src/files/upwind_fdsbp.jl index d103c4ddc70..712bd1dc8e8 100644 --- a/docs/literate/src/files/upwind_fdsbp.jl +++ b/docs/literate/src/files/upwind_fdsbp.jl @@ -30,3 +30,31 @@ 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 Pratical Introduction" of Eleuterio F. Toro (2009), +# [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761). 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) From 8308631960cc932ae4760aa04f16ca0b202a1c0a Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 30 Nov 2022 10:14:07 +0100 Subject: [PATCH 4/8] add to make.jl --- docs/make.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/make.jl b/docs/make.jl index d8d1298fba6..4e0980e6e10 100644 --- 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", From dd200635bc5b791b96d3395c68382f4269ad1fda Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 30 Nov 2022 12:28:59 +0100 Subject: [PATCH 5/8] fix tutorials? --- docs/literate/src/files/upwind_fdsbp.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/literate/src/files/upwind_fdsbp.jl b/docs/literate/src/files/upwind_fdsbp.jl index 712bd1dc8e8..63c1c1d9759 100644 --- a/docs/literate/src/files/upwind_fdsbp.jl +++ b/docs/literate/src/files/upwind_fdsbp.jl @@ -7,12 +7,13 @@ # 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. -# This yields an object representing the operator efficiently. If you want to +# 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) From fc3c5916358742be32a6a4e9028ff0918f61e27b Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 30 Nov 2022 13:18:03 +0100 Subject: [PATCH 6/8] Update upwind_fdsbp.jl --- docs/literate/src/files/upwind_fdsbp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/upwind_fdsbp.jl b/docs/literate/src/files/upwind_fdsbp.jl index 63c1c1d9759..ef8bc3b1db1 100644 --- a/docs/literate/src/files/upwind_fdsbp.jl +++ b/docs/literate/src/files/upwind_fdsbp.jl @@ -38,7 +38,7 @@ Matrix(D_upw.plus) # ```math # f(u) = f^-(u) + f^+(u) # ``` -# such that $f^-(u)$ is associated with left-going waves and `f^+(u)` with +# 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 From d24fa99d170b1504ad1ba4f8685564807c55ae97 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 1 Dec 2022 15:32:14 +0100 Subject: [PATCH 7/8] Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper --- docs/literate/src/files/upwind_fdsbp.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/literate/src/files/upwind_fdsbp.jl b/docs/literate/src/files/upwind_fdsbp.jl index ef8bc3b1db1..0ec0854e080 100644 --- a/docs/literate/src/files/upwind_fdsbp.jl +++ b/docs/literate/src/files/upwind_fdsbp.jl @@ -51,8 +51,8 @@ Matrix(D_upw.plus) # 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 Pratical Introduction" of Eleuterio F. Toro (2009), -# [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761). Such a well-known +# 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 From fe6390f898635b46d70c07a0ea1802b03f2e14f2 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 1 Dec 2022 15:35:08 +0100 Subject: [PATCH 8/8] explain bias --- docs/literate/src/files/upwind_fdsbp.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/docs/literate/src/files/upwind_fdsbp.jl b/docs/literate/src/files/upwind_fdsbp.jl index 0ec0854e080..36ca1b57404 100644 --- a/docs/literate/src/files/upwind_fdsbp.jl +++ b/docs/literate/src/files/upwind_fdsbp.jl @@ -24,7 +24,10 @@ 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. +# 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)