From 1da692e2f8db707c9b811fe9aed5a42b06db74c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 27 Nov 2024 18:59:53 +0100 Subject: [PATCH 1/2] Added non-conservative terms to Cartesian SWE 3D --- ...wwater_cubed_sphere_shell_EC_correction.jl | 14 ++-- ...wwater_cubed_sphere_shell_EC_projection.jl | 12 ++- ...allowwater_cubed_sphere_shell_advection.jl | 9 +- ...hallowwater_cubed_sphere_shell_standard.jl | 7 +- src/equations/shallow_water_3d.jl | 84 ++++++++++++++++++- 5 files changed, 113 insertions(+), 13 deletions(-) diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl b/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl index 956c7ac..bd1b789 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl @@ -8,10 +8,11 @@ using TrixiAtmo equations = ShallowWaterEquations3D(gravity_constant = 9.81) -# Create DG solver with polynomial degree = 3 and Wintemeyer et al.'s flux as surface flux +# Create DG solver with polynomial degree = 3 and Fjordholm et al.'s flux as surface flux polydeg = 3 -volume_flux = flux_fjordholm_etal #flux_wintermeyer_etal -surface_flux = flux_fjordholm_etal #flux_wintermeyer_etal #flux_lax_friedrichs +volume_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) + solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) @@ -43,7 +44,7 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio # Compute the velocity of a rotating solid body # (alpha is the angle between the rotation axis and the polar axis of the spherical coordinate system) v0 = 1.0 # Velocity at the "equator" - alpha = 0.0 #pi / 4 + alpha = pi / 4 v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) vlat = -v0 * sin(phi) * sin(alpha) @@ -52,7 +53,10 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio v2 = -cos(colat) * sin(phi) * vlat + cos(phi) * v_long v3 = sin(colat) * vlat - return prim2cons(SVector(rho, v1, v2, v3, 0), equations) + # Non-constant topography + b = 0.1 * cos(10 * lat) + + return prim2cons(SVector(rho, v1, v2, v3, b), equations) end # Source term function to apply the entropy correction term and the Lagrange multiplier to the semi-discretization. diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl b/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl index 21934c5..bbb80e7 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl @@ -12,8 +12,9 @@ equations = ShallowWaterEquations3D(gravity_constant = 9.81) # Create DG solver with polynomial degree = 3 and Wintemeyer et al.'s flux as surface flux polydeg = 3 -volume_flux = flux_wintermeyer_etal # flux_fjordholm_etal -surface_flux = flux_wintermeyer_etal # flux_fjordholm_etal #flux_lax_friedrichs +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) + solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) @@ -45,7 +46,7 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio # Compute the velocity of a rotating solid body # (alpha is the angle between the rotation axis and the polar axis of the spherical coordinate system) v0 = 1.0 # Velocity at the "equator" - alpha = 0.0 #pi / 4 + alpha = pi / 4 v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) vlat = -v0 * sin(phi) * sin(alpha) @@ -54,7 +55,10 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio v2 = -cos(colat) * sin(phi) * vlat + cos(phi) * v_long v3 = sin(colat) * vlat - return prim2cons(SVector(rho, v1, v2, v3, 0), equations) + # Non-constant topography + b = 0.1 * cos(10 * lat) + + return prim2cons(SVector(rho, v1, v2, v3, b), equations) end initial_condition = initial_condition_advection_sphere diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl index 0749c83..7de23ce 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl @@ -8,13 +8,18 @@ using TrixiAtmo cells_per_dimension = 5 initial_condition = initial_condition_gaussian +polydeg = 3 # We use the ShallowWaterEquations3D equations structure but modify the rhs! function to # convert it to a variable-coefficient advection equation equations = ShallowWaterEquations3D(gravity_constant = 0.0) -# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) +volume_flux = (flux_central, flux_nonconservative_fjordholm_etal) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal) + +solver = DGSEM(polydeg = polydeg, + surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # Source term function to transform the Euler equations into a linear advection equation with variable advection velocity function source_terms_convert_to_linear_advection(u, du, x, t, diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl b/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl index ba1c3bf..71d597c 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl @@ -10,7 +10,12 @@ equations = ShallowWaterEquations3D(gravity_constant = 9.81) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs as surface flux polydeg = 3 -solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs) +volume_flux = (flux_central, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_wintermeyer_etal) + +solver = DGSEM(polydeg = polydeg, + surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # Initial condition for a Gaussian density profile with constant pressure # and the velocity of a rotating solid body diff --git a/src/equations/shallow_water_3d.jl b/src/equations/shallow_water_3d.jl index 8340868..897f467 100644 --- a/src/equations/shallow_water_3d.jl +++ b/src/equations/shallow_water_3d.jl @@ -61,7 +61,7 @@ function ShallowWaterEquations3D(; gravity_constant, H0 = zero(gravity_constant) ShallowWaterEquations3D(gravity_constant, H0) end -Trixi.have_nonconservative_terms(::ShallowWaterEquations3D) = False() # Deactivate non-conservative terms for the moment... +Trixi.have_nonconservative_terms(::ShallowWaterEquations3D) = True() # Deactivate non-conservative terms for the moment... Trixi.varnames(::typeof(cons2cons), ::ShallowWaterEquations3D) = ("h", "h_v1", "h_v2", "h_v3", "b") # Note, we use the total water height, H = h + b, as the first primitive variable for easier @@ -159,6 +159,88 @@ Further details are available in Theorem 1 of the paper: return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) end +""" + flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquations2D) + flux_nonconservative_wintermeyer_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterEquations2D) + +Non-symmetric two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref). + +For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can +be used to ensure well-balancedness and entropy conservation. + +Further details are available in the papers: +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) + An entropy stable nodal discontinuous Galerkin method for the two dimensional + shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry + [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) +""" +@inline function Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterEquations3D) + # Pull the necessary left and right state information + h_ll = waterheight(u_ll, equations) + b_jump = u_rr[5] - u_ll[5] + + # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) + return SVector(0, + normal_direction[1] * equations.gravity * h_ll * b_jump, + normal_direction[2] * equations.gravity * h_ll * b_jump, + normal_direction[3] * equations.gravity * h_ll * b_jump, + 0) +end + +""" + flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquations2D) + flux_nonconservative_fjordholm_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterEquations2D) + +Non-symmetric two-point surface flux discretizing the nonconservative (source) term of +that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref). + +This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy +conservation and well-balancedness. + +Further details for the original finite volume formulation are available in +- Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011) + Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography + [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042) +and for curvilinear 2D case in the paper: +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) + An entropy stable nodal discontinuous Galerkin method for the two dimensional + shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry + [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +""" +@inline function Trixi.flux_nonconservative_fjordholm_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterEquations3D) + # Pull the necessary left and right state information + h_ll, _, _, _, b_ll = u_ll + h_rr, _, _, _, b_rr = u_rr + + h_average = 0.5f0 * (h_ll + h_rr) + b_jump = b_rr - b_ll + + # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) + f2 = normal_direction[1] * equations.gravity * h_average * b_jump + f3 = normal_direction[2] * equations.gravity * h_average * b_jump + f4 = normal_direction[3] * equations.gravity * h_average * b_jump + + # First and last equations do not have a nonconservative flux + f1 = f5 = 0 + + return SVector(f1, f2, f3, f4, f5) +end + """ flux_fjordholm_etal(u_ll, u_rr, orientation, equations::ShallowWaterEquations1D) From 60094635be4ec53062892b770165b833af88fa9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 28 Nov 2024 10:43:57 +0100 Subject: [PATCH 2/2] Added well-balancing test and corrected reference results for EC tests (already tested that EC works) --- ...water_cubed_sphere_shell_well_balancing.jl | 96 +++++++++++++++++++ src/equations/shallow_water_3d.jl | 2 +- test/test_3d_shallow_water.jl | 67 +++++++++---- 3 files changed, 144 insertions(+), 21 deletions(-) create mode 100644 examples/elixir_shallowwater_cubed_sphere_shell_well_balancing.jl diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_well_balancing.jl b/examples/elixir_shallowwater_cubed_sphere_shell_well_balancing.jl new file mode 100644 index 0000000..f60a9f8 --- /dev/null +++ b/examples/elixir_shallowwater_cubed_sphere_shell_well_balancing.jl @@ -0,0 +1,96 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiAtmo + +############################################################################### +# Entropy conservation for the spherical shallow water equations in Cartesian +# form obtained through the projection of the momentum onto the divergence-free +# tangential contravariant vectors + +equations = ShallowWaterEquations3D(gravity_constant = 9.81) + +# Create DG solver with polynomial degree = 3 and Wintemeyer et al.'s flux as surface flux +polydeg = 3 +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) + +solver = DGSEM(polydeg = polydeg, + surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +# Initial condition for a Gaussian density profile with constant pressure +# and the velocity of a rotating solid body +function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquations3D) + # Constant water height + H = 1.0 + + # Non-constant topography + lat = asin(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2)) + b = 0.1 * cos(10 * lat) + + return prim2cons(SVector(H, 0, 0, 0, b), equations) +end + +initial_condition = initial_condition_advection_sphere + +mesh = P4estMeshCubedSphere2D(5, 2.0, polydeg = polydeg, + initial_refinement_level = 0) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + metric_terms = MetricTermsInvariantCurl(), + source_terms = source_terms_lagrange_multiplier) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to π +tspan = (0.0, pi) +ode = semidiscretize(semi, tspan) + +# Clean the initial condition +for element in eachelement(solver, semi.cache) + for j in eachnode(solver), i in eachnode(solver) + u0 = Trixi.wrap_array(ode.u0, semi) + + contravariant_normal_vector = Trixi.get_contravariant_vector(3, + semi.cache.elements.contravariant_vectors, + i, j, element) + clean_solution_lagrange_multiplier!(u0[:, i, j, element], equations, + contravariant_normal_vector) + end +end + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 10, + save_analysis = true, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight, energy_total), + analysis_polydeg = polydeg) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 10, + solution_variables = cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.7) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/equations/shallow_water_3d.jl b/src/equations/shallow_water_3d.jl index 897f467..fd269de 100644 --- a/src/equations/shallow_water_3d.jl +++ b/src/equations/shallow_water_3d.jl @@ -61,7 +61,7 @@ function ShallowWaterEquations3D(; gravity_constant, H0 = zero(gravity_constant) ShallowWaterEquations3D(gravity_constant, H0) end -Trixi.have_nonconservative_terms(::ShallowWaterEquations3D) = True() # Deactivate non-conservative terms for the moment... +Trixi.have_nonconservative_terms(::ShallowWaterEquations3D) = True() Trixi.varnames(::typeof(cons2cons), ::ShallowWaterEquations3D) = ("h", "h_v1", "h_v2", "h_v3", "b") # Note, we use the total water height, H = h + b, as the first primitive variable for easier diff --git a/test/test_3d_shallow_water.jl b/test/test_3d_shallow_water.jl index 0287a35..6560956 100644 --- a/test/test_3d_shallow_water.jl +++ b/test/test_3d_shallow_water.jl @@ -11,18 +11,18 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_cubed_sphere_shell_EC_correction.jl"), l2=[ - 7.23800458e-02, - 9.98871590e-02, - 4.55606969e-02, - 3.17422064e-02, - 0.00000000e+00 + 0.07271743465440743, + 0.07435870820933804, + 0.05059253372287277, + 0.07712245696955002, + 0.002166321001178188 ], linf=[ - 1.05686060e+00, - 1.04413842e+00, - 3.12374228e-01, - 3.19064636e-01, - 0.00000000e+00 + 1.057593830959098, + 0.7420535585104321, + 0.33165510628037276, + 0.6538582658635804, + 0.009155117525905546 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -38,18 +38,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_cubed_sphere_shell_EC_projection.jl"), l2=[ - 7.25131271e-02, - 1.01017589e-01, - 4.39697415e-02, - 3.08898940e-02, - 0.00000000e+00 + 0.07281376793171955, + 0.07425222671946745, + 0.05069528390329348, + 0.07741804767929857, + 0.002166321001178188 ], linf=[ - 1.06007536e+00, - 1.05786719e+00, - 3.93826726e-01, - 2.34129714e-01, - 0.00000000e+00 + 1.0576298065256564, + 0.7245020984321977, + 0.3273888680444672, + 0.6494576573261298, + 0.009155117525905546 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -88,4 +88,31 @@ end end end +@trixiatmo_testset "elixir_shallowwater_cubed_sphere_shell_well_balancing" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_cubed_sphere_shell_well_balancing.jl"), + l2=[ + 0.00000000e+00, + 0.00000000e+00, + 0.00000000e+00, + 0.00000000e+00, + 0.00000000e+00 + ], + linf=[ + 0.00000000e+00, + 0.00000000e+00, + 0.00000000e+00, + 0.00000000e+00, + 0.00000000e+00 + ], atol=8.0e-13) # Needs a slightly larger tolerance for linf + # 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 TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + end # module