From 5589d2a4366ccba5cab5487676fd17ef42db3375 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 15 Jul 2023 12:38:45 +0200 Subject: [PATCH 1/3] enable fully periodic upwind SBP oeprators --- Project.toml | 2 +- .../tree_1d_fdsbp/elixir_advection_upwind.jl | 3 +- .../elixir_advection_upwind_periodic.jl | 57 +++++++++++++++++++ src/solvers/fdsbp_tree/fdsbp.jl | 2 + src/solvers/fdsbp_tree/fdsbp_1d.jl | 28 ++++++++- src/solvers/fdsbp_tree/fdsbp_2d.jl | 28 ++++++++- src/solvers/fdsbp_tree/fdsbp_3d.jl | 28 ++++++++- test/test_tree_1d_fdsbp.jl | 15 +++++ 8 files changed, 155 insertions(+), 8 deletions(-) create mode 100644 examples/tree_1d_fdsbp/elixir_advection_upwind_periodic.jl diff --git a/Project.toml b/Project.toml index 6c3c7fa020..a49cfb2e25 100644 --- 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 index 5c50e1a6c6..18dd818e3c 100644 --- 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 new file mode 100644 index 0000000000..3eb805095f --- /dev/null +++ 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 index cbb6fd1624..94d3a2397b 100644 --- a/src/solvers/fdsbp_tree/fdsbp.jl +++ b/src/solvers/fdsbp_tree/fdsbp.jl @@ -27,6 +27,8 @@ The other arguments have the same meaning as in [`DG`](@ref) or [`DGSEM`](@ref). """ const FDSBP = DG{Basis} where {Basis <: AbstractDerivativeOperator} +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 index c771207494..0de0cff485 100644 --- 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 index 241e0d9534..beff605629 100644 --- 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 index a4f69d3d48..0c3f18b6d6 100644 --- 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 index 118385c34b..ce0ca660d3 100644 --- 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 From 4da632c35e48d8b5eb4d0e93eaa77b2dfb62d22c Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 15 Jul 2023 12:49:33 +0200 Subject: [PATCH 2/3] 2D and 3D tests --- test/test_tree_2d_fdsbp.jl | 18 ++++++++++++++++++ test/test_tree_3d_fdsbp.jl | 23 ++++++++++++++++++++++- 2 files changed, 40 insertions(+), 1 deletion(-) diff --git a/test/test_tree_2d_fdsbp.jl b/test/test_tree_2d_fdsbp.jl index 7c58ef89a6..e81c82f3f3 100644 --- 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 index 9dceab3803..106dd007b0 100644 --- 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], From 2bfcb1af3740236f573fc379b3a9fdb1b833aa0b Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 15 Jul 2023 14:15:36 +0200 Subject: [PATCH 3/3] comment on PeriodicFDSBP --- src/solvers/fdsbp_tree/fdsbp.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/solvers/fdsbp_tree/fdsbp.jl b/src/solvers/fdsbp_tree/fdsbp.jl index 94d3a2397b..11b09c6df9 100644 --- a/src/solvers/fdsbp_tree/fdsbp.jl +++ b/src/solvers/fdsbp_tree/fdsbp.jl @@ -27,6 +27,7 @@ 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)