WIP: Add non-conservative terms to 3D Cartesian SWE

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.

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

examples/elixir_shallowwater_cubed_sphere_shell_advection.jl

@@ -17,14 +17,19 @@ using TrixiAtmo
 # semidiscretization of the linear advection equation
 
 initial_condition = initial_condition_gaussian
+polydeg = 3
 cells_per_dimension = (5, 5)
 
 # 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,

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

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()

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()
 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)

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