diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_convergence_sc_subcell_curved.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_convergence_sc_subcell_curved.jl new file mode 100644 index 00000000..cc198fdc --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_convergence_sc_subcell_curved.jl @@ -0,0 +1,134 @@ +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater +using Symbolics + +############################################################################### +# Semidiscretization of the multilayer shallow water equations with a single layer to test +# convergence. The initial condition and source term are created using the +# method of manufactured solutions (MMS) with the help of Symbolics.jl. + +equations = ShallowWaterMultiLayerEquations2D(gravity = 1.1, rhos = (1.0)) + +### Create manufactured solution for method of manufactured solutions (MMS) + +# Symbolic Variables +@variables x_sym[1:2], t_sym, g + +# Define Differentials +Dt, Dx, Dy = Differential(t_sym), Differential(x_sym[1]), Differential(x_sym[2]) + +## Initial condition +############################################################################### +# Primitive Variables +ω = pi * 1.0 +H = 4.0 + 0.2 * cos(ω * x_sym[1] + t_sym) + 0.2 * cos(ω * x_sym[2] + t_sym) +v1 = 0.5 +v2 = 0.5 +b = 1.0 + 0.2 * cos(ω * x_sym[1]) + 0.2 * cos(ω * x_sym[2]) +h = H - b + +init = [H, v1, v2, b] + +## PDE +############################################################################### +eqs = [ + Dt(h) + Dx(h * v1) + Dy(h * v2), + Dt(h * v1) + Dx(h * v1^2 + 0.5 * g * h^2) + Dy(h * v1 * v2) + g * h * Dx(b), + Dt(h * v2) + Dx(h * v1 * v2) + Dy(h * v2^2 + 0.5 * g * h^2) + g * h * Dy(b), + 0 +] + +## Create the functions for the manufactured solution +######################################################################################## +# Expand derivatives +du_exprs = expand_derivatives.(eqs) + +# Build functions +du_funcs = build_function.(du_exprs, Ref(x_sym), t_sym, g, expression = Val(false)) + +init_funcs = build_function.(init, Ref(x_sym), t_sym, expression = Val(false)) + +function initial_condition_convergence_mms(x, + t, + equations::ShallowWaterMultiLayerEquations2D) + prim = SVector{4, Float64}([f(x, t) for f in init_funcs]...) + return prim2cons(prim, equations) +end + +function source_terms_convergence_mms(u, + x, + t, + equations::ShallowWaterMultiLayerEquations2D) + g = equations.gravity + return SVector{4, Float64}([f(x, t, g) for f in du_funcs]...) +end + +initial_condition = initial_condition_convergence_mms + +############################################################################### +# Get the DG approximation space + +polydeg = 3 +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump) +surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump) + +basis = LobattoLegendreBasis(polydeg) +limiter_idp = SubcellLimiterIDP(equations, basis;) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Get the P4estMesh and setup a periodic mesh on the domain [-1,1]^2 with an affine type mapping to +# obtain a warped curvilinear mesh. +function mapping_twist(xi, eta) + y = eta + 0.1 * sin(pi * xi) * cos(0.5 * pi * eta) + x = xi + 0.1 * sin(pi * eta) * cos(0.5 * pi * xi) + return SVector(x, y) +end + +# Create warped P4estMesh with 8 x 8 elements +trees_per_dimension = (2, 2) +mesh = P4estMesh(trees_per_dimension, polydeg = 3, + initial_refinement_level = 2, + periodicity = true, + mapping = mapping_twist) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_mms, + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.1) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 0.7) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(),) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + ode_default_options()..., + callback = callbacks, adaptive = false); diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_ec_sc_subcell.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_ec_sc_subcell.jl new file mode 100644 index 00000000..4a4ee8a5 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_ec_sc_subcell.jl @@ -0,0 +1,110 @@ +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the multilayer shallow water equations with a single layer and a bottom +# topography function for a blast wave test with discontinuous initial conditions to test entropy +# conservation with the subcell limiting volume integral on a curvilinear mesh. + +equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, rhos = (1.0)) + +function initial_condition_weak_blast_wave(x, t, + equations::ShallowWaterMultiLayerEquations2D) + # Set up polar coordinates + RealT = eltype(x) + inicenter = SVector(convert(RealT, 0.0), convert(RealT, 0.0)) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Calculate primitive variables + H = r > 0.2f0 ? 3.8f0 : 4.0f0 + v1 = r > 0.2f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi + v2 = r > 0.2f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi + v1 = 0.0 + v2 = 0.0 + + # Setup a continuous bottom topography + r = 0.4 + b = (((x[1])^2 + (x[2])^2) < r^2 ? + 0.5 * (cos(1 / r * pi * sqrt((x[1])^2 + (x[2])^2)) + 1) : 0.0) + + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +############################################################################### +# Get the DG approximation space + +polydeg = 3 +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump) +surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump) + +basis = LobattoLegendreBasis(polydeg) +# For reproducibility set the limiting coefficients to pure DG, +# see https://github.com/trixi-framework/Trixi.jl/pull/2007 +limiter_idp = SubcellLimiterIDP(equations, basis;) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Get the P4estMesh and setup a periodic mesh on the domain [-1,1]^2 with an affine type mapping to +# obtain a warped curvilinear mesh. +function mapping_twist(xi, eta) + y = eta + 0.1 * sin(pi * xi) * cos(0.5 * pi * eta) + x = xi + 0.1 * sin(pi * eta) * cos(0.5 * pi * xi) + return SVector(x, y) +end + +# Create P4estMesh with 8 x 8 elements +trees_per_dimension = (2, 2) +mesh = P4estMesh(trees_per_dimension, polydeg = 2, + initial_refinement_level = 2, + periodicity = true, + mapping = mapping_twist) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solver + +tspan = (0.0, 0.2) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + analysis_polydeg = polydeg, + extra_analysis_errors = (:conservation_error,)) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.1, + save_initial_solution = true, + save_final_solution = true, + extra_node_variables = (:limiting_coefficient,)) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### + +# Setup the stage callbacks +stage_callbacks = (SubcellLimiterIDPCorrection(),) + +# run the simulation +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + ode_default_options()..., + callback = callbacks, adaptive = false); diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl index 5d5be01e..85337cb6 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl @@ -44,7 +44,7 @@ ode = semidiscretize(semi, tspan) ############################################################################### # Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. -# alternative version of the initial conditinon used to setup a truly discontinuous +# alternative version of the initial condition used to setup a truly discontinuous # bottom topography function and initial condition for this academic testcase of entropy conservation. # The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff # In contrast to the usual signature of initial conditions, this one get passed the diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index 71a24fcb..b8c0bc3c 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -624,6 +624,12 @@ or jump portion of the non-conservative flux based on the type of the nonconservative_type argument, employing multiple dispatch. They are used to compute the subcell fluxes in [`Trixi.VolumeIntegralSubcellLimiting`](@extref). +!!! note + While `flux_nonconservative_ersing_etal_local_jump` and [`flux_nonconservative_ersing_etal`](@ref) + are equivalent at interfaces, the volume formulation for [`Trixi.VolumeIntegralSubcellLimiting`](@extref) + differs as `flux_nonconservative_ersing_etal_local_jump` applies the local normal direction + instead of the averaged one. + ## References - Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf. @@ -653,7 +659,7 @@ end nonconservative_type::Trixi.NonConservativeLocal, nonconservative_term::Integer) flux_nonconservative_ersing_local_jump(u_ll, u_rr, - normal_direction::AbstractVector, + normal_direction_ll::AbstractVector, equations::ShallowWaterMultiLayerEquations2D, nonconservative_type::Trixi.NonConservativeLocal, nonconservative_term::Integer) @@ -695,7 +701,7 @@ Local part of the nonconservative term needed for the calculation of the non-con end @inline function flux_nonconservative_ersing_etal_local_jump(u_ll, - normal_direction::AbstractVector, + normal_direction_ll::AbstractVector, equations::ShallowWaterMultiLayerEquations2D, nonconservative_type::Trixi.NonConservativeLocal, nonconservative_term::Integer) @@ -714,7 +720,9 @@ end f_h = zero(real(equations)) f_hv = g * h_ll[i] - setlayer!(f, f_h, f_hv, f_hv, i, equations) + setlayer!(f, f_h, f_hv * normal_direction_ll[1], f_hv * normal_direction_ll[2], + i, + equations) end return SVector(f) @@ -812,8 +820,8 @@ end f_hv += h_jump[j] end end - setlayer!(f, f_h, f_hv * normal_direction[1], - f_hv * normal_direction[2], i, equations) + setlayer!(f, f_h, f_hv, + f_hv, i, equations) end return SVector(f) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index fbc2249b..ad148ec7 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -165,6 +165,25 @@ end # SWE @test_allocations(Trixi.rhs!, semi, sol, 15000) end + @trixi_testset "elixir_shallowwater_multilayer_convergence_sc_subcell_curved.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_convergence_sc_subcell_curved.jl"), + l2=[ + 0.00013275624153302434, + 0.0001586600356395913, + 0.000158660035639501, + 2.9583315272612922e-5 + ], + linf=[ + 0.0007544272991792944, + 0.0007250877164874936, + 0.00072508771648927, + 9.31948456051046e-5 + ]) + # Allocation testing is disabled as the symbolic source term computation is known to cause + # allocations. + end + @trixi_testset "elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl"), @@ -319,6 +338,30 @@ end # SWE # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test_allocations(Trixi.rhs!, semi, sol, 15000) end + + @trixi_testset "elixir_shallowwater_multilayer_ec_sc_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_ec_sc_subcell.jl"), + l2=[ + 0.04338581073333642, + 0.12068863355590975, + 0.12068863355590971, + 7.906244739074657e-18 + ], + linf=[ + 0.2550462084344076, + 0.42004261172048024, + 0.4200426117204863, + 1.1102230246251565e-16 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from + # OrdinaryDiffEq.jl + # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 + @test_allocations(Trixi.rhs!, semi, sol, 15000) + end end # MLSWE end # P4estMesh2D