allow periodic FDSBP operators (#1570)

* enable fully periodic upwind SBP oeprators * 2D and 3D tests * comment on PeriodicFDSBP
trixi-framework · Jul 15, 2023 · 932f433 · 932f433
1 parent 07981c3
commit 932f433
Show file tree

Hide file tree

Showing 10 changed files with 196 additions and 9 deletions.
diff --git a/Project.toml b/Project.toml
@@ -79,7 +79,7 @@ StaticArrayInterface = "1.4"
 StaticArrays = "1"
 StrideArrays = "0.1.18"
 StructArrays = "0.6"
-SummationByPartsOperators = "0.5.25"
+SummationByPartsOperators = "0.5.41"
 TimerOutputs = "0.5"
 Triangulate = "2.0"
 TriplotBase = "0.1"

diff --git a/examples/tree_1d_fdsbp/elixir_advection_upwind.jl b/examples/tree_1d_fdsbp/elixir_advection_upwind.jl
@@ -27,7 +27,8 @@ coordinates_min = -1.0
 coordinates_max =  1.0
 mesh = TreeMesh(coordinates_min, coordinates_max,
                 initial_refinement_level = 4,
-                n_cells_max = 10_000)
+                n_cells_max = 10_000,
+                periodicity = true)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_sin, solver)
 

diff --git a/examples/tree_1d_fdsbp/elixir_advection_upwind_periodic.jl b/examples/tree_1d_fdsbp/elixir_advection_upwind_periodic.jl
@@ -0,0 +1,57 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear scalar advection equation equation
+
+equations = LinearScalarAdvectionEquation1D(1.0)
+
+function initial_condition_sin(x, t, equation::LinearScalarAdvectionEquation1D)
+    return SVector(sinpi(x[1] - equations.advection_velocity[1] * t))
+end
+
+D_upw = upwind_operators(SummationByPartsOperators.periodic_derivative_operator,
+                         accuracy_order = 4,
+                         xmin = -1.0, xmax = 1.0,
+                         N = 64)
+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 = 0,
+                n_cells_max = 10_000,
+                periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_sin, 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)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6,
+            ode_default_options()..., callback=callbacks);
+summary_callback() # print the timer summary
diff --git a/src/solvers/fdsbp_tree/fdsbp.jl b/src/solvers/fdsbp_tree/fdsbp.jl
@@ -27,6 +27,9 @@ The other arguments have the same meaning as in [`DG`](@ref) or [`DGSEM`](@ref).
 """
 const FDSBP = DG{Basis} where {Basis <: AbstractDerivativeOperator}
 
+# Internal abbreviation for easier-to-read dispatch (not exported)
+const PeriodicFDSBP = FDSBP{Basis} where {Basis <: AbstractPeriodicDerivativeOperator}
+
 function FDSBP(D_SBP::AbstractDerivativeOperator; surface_integral, volume_integral)
     # `nothing` is passed as `mortar`
     return DG(D_SBP, nothing, surface_integral, volume_integral)

diff --git a/src/solvers/fdsbp_tree/fdsbp_1d.jl b/src/solvers/fdsbp_tree/fdsbp_1d.jl
@@ -165,6 +165,14 @@ function calc_surface_integral!(du, u, mesh::TreeMesh{1},
     return nothing
 end
 
+# Periodic FDSBP operators need to use a single element without boundaries
+function calc_surface_integral!(du, u, mesh::TreeMesh1D,
+                                equations, surface_integral::SurfaceIntegralStrongForm,
+                                dg::PeriodicFDSBP, cache)
+    @assert nelements(dg, cache) == 1
+    return nothing
+end
+
 # Specialized interface flux computation because the upwind solver does
 # not require a standard numerical flux (Riemann solver). The flux splitting
 # already separates the solution information into right-traveling and
@@ -239,13 +247,25 @@ function calc_surface_integral!(du, u, mesh::TreeMesh{1},
     return nothing
 end
 
+# Periodic FDSBP operators need to use a single element without boundaries
+function calc_surface_integral!(du, u, mesh::TreeMesh1D,
+                                equations, surface_integral::SurfaceIntegralUpwind,
+                                dg::PeriodicFDSBP, cache)
+    @assert nelements(dg, cache) == 1
+    return nothing
+end
+
 # AnalysisCallback
 
 function integrate_via_indices(func::Func, u,
                                mesh::TreeMesh{1}, equations,
                                dg::FDSBP, cache, args...; normalize = true) where {Func}
     # TODO: FD. This is rather inefficient right now and allocates...
-    weights = diag(SummationByPartsOperators.mass_matrix(dg.basis))
+    M = SummationByPartsOperators.mass_matrix(dg.basis)
+    if M isa UniformScaling
+        M = M(nnodes(dg))
+    end
+    weights = diag(M)
 
     # Initialize integral with zeros of the right shape
     integral = zero(func(u, 1, 1, equations, dg, args...))
@@ -271,7 +291,11 @@ function calc_error_norms(func, u, t, analyzer,
                           mesh::TreeMesh{1}, equations, initial_condition,
                           dg::FDSBP, cache, cache_analysis)
     # TODO: FD. This is rather inefficient right now and allocates...
-    weights = diag(SummationByPartsOperators.mass_matrix(dg.basis))
+    M = SummationByPartsOperators.mass_matrix(dg.basis)
+    if M isa UniformScaling
+        M = M(nnodes(dg))
+    end
+    weights = diag(M)
     @unpack node_coordinates = cache.elements
 
     # Set up data structures

diff --git a/src/solvers/fdsbp_tree/fdsbp_2d.jl b/src/solvers/fdsbp_tree/fdsbp_2d.jl
@@ -201,6 +201,14 @@ function calc_surface_integral!(du, u, mesh::TreeMesh{2},
     return nothing
 end
 
+# Periodic FDSBP operators need to use a single element without boundaries
+function calc_surface_integral!(du, u, mesh::TreeMesh2D,
+                                equations, surface_integral::SurfaceIntegralStrongForm,
+                                dg::PeriodicFDSBP, cache)
+    @assert nelements(dg, cache) == 1
+    return nothing
+end
+
 # Specialized interface flux computation because the upwind solver does
 # not require a standard numerical flux (Riemann solver). The flux splitting
 # already separates the solution information into right-traveling and
@@ -295,12 +303,24 @@ function calc_surface_integral!(du, u, mesh::TreeMesh{2},
     return nothing
 end
 
+# Periodic FDSBP operators need to use a single element without boundaries
+function calc_surface_integral!(du, u, mesh::TreeMesh2D,
+                                equations, surface_integral::SurfaceIntegralUpwind,
+                                dg::PeriodicFDSBP, cache)
+    @assert nelements(dg, cache) == 1
+    return nothing
+end
+
 # AnalysisCallback
 function integrate_via_indices(func::Func, u,
                                mesh::TreeMesh{2}, equations,
                                dg::FDSBP, cache, args...; normalize = true) where {Func}
     # TODO: FD. This is rather inefficient right now and allocates...
-    weights = diag(SummationByPartsOperators.mass_matrix(dg.basis))
+    M = SummationByPartsOperators.mass_matrix(dg.basis)
+    if M isa UniformScaling
+        M = M(nnodes(dg))
+    end
+    weights = diag(M)
 
     # Initialize integral with zeros of the right shape
     integral = zero(func(u, 1, 1, 1, equations, dg, args...))
@@ -326,7 +346,11 @@ function calc_error_norms(func, u, t, analyzer,
                           mesh::TreeMesh{2}, equations, initial_condition,
                           dg::FDSBP, cache, cache_analysis)
     # TODO: FD. This is rather inefficient right now and allocates...
-    weights = diag(SummationByPartsOperators.mass_matrix(dg.basis))
+    M = SummationByPartsOperators.mass_matrix(dg.basis)
+    if M isa UniformScaling
+        M = M(nnodes(dg))
+    end
+    weights = diag(M)
     @unpack node_coordinates = cache.elements
 
     # Set up data structures

diff --git a/src/solvers/fdsbp_tree/fdsbp_3d.jl b/src/solvers/fdsbp_tree/fdsbp_3d.jl
@@ -237,6 +237,14 @@ function calc_surface_integral!(du, u, mesh::TreeMesh{3},
     return nothing
 end
 
+# Periodic FDSBP operators need to use a single element without boundaries
+function calc_surface_integral!(du, u, mesh::TreeMesh3D,
+                                equations, surface_integral::SurfaceIntegralStrongForm,
+                                dg::PeriodicFDSBP, cache)
+    @assert nelements(dg, cache) == 1
+    return nothing
+end
+
 # Specialized interface flux computation because the upwind solver does
 # not require a standard numerical flux (Riemann solver). The flux splitting
 # already separates the solution information into right-traveling and
@@ -346,13 +354,25 @@ function calc_surface_integral!(du, u, mesh::TreeMesh{3},
     return nothing
 end
 
+# Periodic FDSBP operators need to use a single element without boundaries
+function calc_surface_integral!(du, u, mesh::TreeMesh3D,
+                                equations, surface_integral::SurfaceIntegralUpwind,
+                                dg::PeriodicFDSBP, cache)
+    @assert nelements(dg, cache) == 1
+    return nothing
+end
+
 # AnalysisCallback
 
 function integrate_via_indices(func::Func, u,
                                mesh::TreeMesh{3}, equations,
                                dg::FDSBP, cache, args...; normalize = true) where {Func}
     # TODO: FD. This is rather inefficient right now and allocates...
-    weights = diag(SummationByPartsOperators.mass_matrix(dg.basis))
+    M = SummationByPartsOperators.mass_matrix(dg.basis)
+    if M isa UniformScaling
+        M = M(nnodes(dg))
+    end
+    weights = diag(M)
 
     # Initialize integral with zeros of the right shape
     integral = zero(func(u, 1, 1, 1, 1, equations, dg, args...))
@@ -378,7 +398,11 @@ function calc_error_norms(func, u, t, analyzer,
                           mesh::TreeMesh{3}, equations, initial_condition,
                           dg::FDSBP, cache, cache_analysis)
     # TODO: FD. This is rather inefficient right now and allocates...
-    weights = diag(SummationByPartsOperators.mass_matrix(dg.basis))
+    M = SummationByPartsOperators.mass_matrix(dg.basis)
+    if M isa UniformScaling
+        M = M(nnodes(dg))
+    end
+    weights = diag(M)
     @unpack node_coordinates = cache.elements
 
     # Set up data structures

diff --git a/test/test_tree_1d_fdsbp.jl b/test/test_tree_1d_fdsbp.jl
@@ -23,6 +23,21 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_fdsbp")
       @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
   end
+
+  @trixi_testset "elixir_advection_upwind_periodic.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_upwind_periodic.jl"),
+      l2   = [1.1672962783692568e-5],
+      linf = [1.650514414558435e-5])
+
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+      t = sol.t[end]
+      u_ode = sol.u[end]
+      du_ode = similar(u_ode)
+      @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+  end
 end
 
 @testset "Inviscid Burgers" begin

diff --git a/test/test_tree_2d_fdsbp.jl b/test/test_tree_2d_fdsbp.jl
@@ -23,6 +23,24 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_fdsbp")
       @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
   end
+
+  @trixi_testset "elixir_advection_extended.jl with periodic operators" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
+      l2   = [1.1239649404463432e-5],
+      linf = [1.5895264629195438e-5],
+      D_SBP = SummationByPartsOperators.periodic_derivative_operator(
+        derivative_order = 1, accuracy_order = 4, xmin = 0.0, xmax = 1.0, N = 40),
+      initial_refinement_level = 0)
+
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+      t = sol.t[end]
+      u_ode = sol.u[end]
+      du_ode = similar(u_ode)
+      @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+  end
 end
 
 @testset "Compressible Euler" begin

diff --git a/test/test_tree_3d_fdsbp.jl b/test/test_tree_3d_fdsbp.jl
@@ -7,7 +7,7 @@ include("test_trixi.jl")
 
 EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_3d_fdsbp")
 
-@testset "Compressible Euler" begin
+@testset "Linear scalar advection" begin
   @trixi_testset "elixir_advection_extended.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
       l2   = [0.005355755365412444],
@@ -23,6 +23,27 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_3d_fdsbp")
     end
   end
 
+  @trixi_testset "elixir_advection_extended.jl with periodic operators" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
+      l2   = [1.3819894522373702e-8],
+      linf = [3.381866298113323e-8],
+      D_SBP = SummationByPartsOperators.periodic_derivative_operator(
+        derivative_order = 1, accuracy_order = 4, xmin = 0.0, xmax = 1.0, N = 10),
+      initial_refinement_level = 0,
+      tspan = (0.0, 5.0))
+
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+      t = sol.t[end]
+      u_ode = sol.u[end]
+      du_ode = similar(u_ode)
+      @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+  end
+end
+
+@testset "Compressible Euler" begin
   @trixi_testset "elixir_euler_convergence.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence.jl"),
       l2   = [2.247522803543667e-5, 2.2499169224681058e-5, 2.24991692246826e-5, 2.2499169224684707e-5, 5.814121361417382e-5],