Add bottom topography and Coriolis terms to 3D Cartesian SWE, and incorporate additional test cases (#53)

* Added non-conservative terms to Cartesian SWE 3D

* Added well-balancing test and corrected reference results for EC tests (already tested that EC works)

* Added DissipationLaxFriedrichsEntropyVariables for the 3D Cartesian SWE

* Added hack to enable the use of the weak form kernel for the spherical advection tests

* Moved non-conservative fluxes below conservative fluxes

* Renamed gravity_constant to gravity in elixir

* Correected docstrings

* Modify Cartesian well-balanced test to match with covariant test (Coriolis not yet implemented)

* Added Coriolis source terms

* Added geostrophic test for 3F SW equations

* Updated reference values for standard DGSEM discretization of 3D SWE

* Added isolated mountain test

* Updated EC and ES tests to use the reference initial conditions of the code (and match the covariant tests)

* Renamed elixirs for consistency

* Add rotation rate to well-balanced elixir

* Removed 'standard DGSEM' test (it was not adding anything)

* Improve docstring

* Minor improvements to docstring

---------

Co-authored-by: Tristan Montoya <[email protected]>

Author: amrueda

examples/elixir_shallowwater_cubed_sphere_shell_advection.jl renamed to examples/elixir_shallowwater_cartesian_advection_cubed_sphere.jl

@@ -17,14 +17,16 @@ 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 = 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)
+solver = DGSEM(polydeg = polydeg,
+               surface_flux = (flux_lax_friedrichs, flux_nonconservative_wintermeyer_etal))
 
 # 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,
@@ -41,6 +43,21 @@ function source_terms_convert_to_linear_advection(u, du, x, t,
     return SVector(0.0, s2, s3, s4, 0.0)
 end
 
+# Hack to use the weak form kernel with ShallowWaterEquations3D (a non-conservative equation).
+# This works only because we have a constant bottom topography, so the equations are effectively
+# conservative. Note that the weak form kernel is NOT equal to the flux differencing kernel
+# with central fluxes because of the curved geometry!
+@inline function Trixi.weak_form_kernel!(du, u,
+                                         element,
+                                         mesh::Union{StructuredMesh{2},
+                                                     UnstructuredMesh2D,
+                                                     P4estMesh{2}, T8codeMesh{2}},
+                                         nonconservative_terms::Trixi.True,
+                                         equations::Trixi.AbstractEquations{3},
+                                         dg::DGSEM, cache, alpha = true)
+    Trixi.weak_form_kernel!(du, u, element, mesh, Trixi.False(), equations, dg, cache)
+end
+
 # Create a 2D cubed sphere mesh the size of the Earth
 mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = 3,
                               #initial_refinement_level = 0, # Comment to use cells_per_dimension in the convergence test

examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl renamed to examples/elixir_shallowwater_cartesian_advection_quad_icosahedron.jl

@@ -24,7 +24,8 @@ cells_per_dimension = (2, 2)
 equations = ShallowWaterEquations3D(gravity = 0.0)
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
-solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs)
+solver = DGSEM(polydeg = polydeg,
+               surface_flux = (flux_lax_friedrichs, flux_nonconservative_wintermeyer_etal))
 
 # 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,
@@ -41,6 +42,21 @@ function source_terms_convert_to_linear_advection(u, du, x, t,
     return SVector(0.0, s2, s3, s4, 0.0)
 end
 
+# Hack to use the weak form kernel with ShallowWaterEquations3D (a non-conservative equation).
+# This works only because we have a constant bottom topography, so the equations are effectively
+# conservative. Note that the weak form kernel is NOT equal to the flux differencing kernel
+# with central fluxes because of the curved geometry!
+@inline function Trixi.weak_form_kernel!(du, u,
+                                         element,
+                                         mesh::Union{StructuredMesh{2},
+                                                     UnstructuredMesh2D,
+                                                     P4estMesh{2}, T8codeMesh{2}},
+                                         nonconservative_terms::Trixi.True,
+                                         equations::Trixi.AbstractEquations{3},
+                                         dg::DGSEM, cache, alpha = true)
+    Trixi.weak_form_kernel!(du, u, element, mesh, Trixi.False(), equations, dg, cache)
+end
+
 # Create a 2D quad-based icosahedral mesh the size of the Earth
 mesh = P4estMeshQuadIcosahedron2D(cells_per_dimension[1], EARTH_RADIUS,
                                   #initial_refinement_level = 0,

examples/elixir_shallowwater_cartesian_geostrophic_balance.jl

@@ -0,0 +1,89 @@
+###############################################################################
+# Entropy-stable DGSEM for the 3D shallow water equations in Cartesian form on 
+# the cubed sphere: Steady geostrophic balance (Case 2, Williamson et al., 1992)
+###############################################################################
+
+using OrdinaryDiffEq, Trixi, TrixiAtmo
+
+###############################################################################
+# Parameters
+
+initial_condition = initial_condition_geostrophic_balance
+polydeg = 3
+cells_per_dimension = (5, 5)
+n_saves = 10
+tspan = (0.0, 5.0 * SECONDS_PER_DAY)
+
+###############################################################################
+# Spatial discretization
+
+mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = polydeg,
+                              element_local_mapping = true)
+
+equations = ShallowWaterEquations3D(gravity = EARTH_GRAVITATIONAL_ACCELERATION,
+                                    rotation_rate = EARTH_ROTATION_RATE)
+
+# Create DG solver with polynomial degree = polydeg
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (FluxPlusDissipation(flux_wintermeyer_etal, DissipationLocalLaxFriedrichs()),
+                flux_nonconservative_wintermeyer_etal)
+
+solver = DGSEM(polydeg = polydeg,
+               surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+# Transform the initial condition to the proper set of conservative variables
+initial_condition_transformed = transform_initial_condition(initial_condition, equations)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver,
+                                    source_terms = source_terms_coriolis_lagrange_multiplier)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0 to T
+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 = 200,
+                                     save_analysis = true,
+                                     extra_analysis_errors = (:conservation_error,))
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves,
+                                     solution_variables = cons2cons)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 0.4)
+
+# 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 = 100.0, save_everystep = false, callback = callbacks)

examples/elixir_shallowwater_cartesian_isolated_mountain.jl

@@ -0,0 +1,83 @@
+###############################################################################
+# Entropy-stable DGSEM for the 3D shallow water equations in Cartesian form on 
+# the cubed sphere: Zonal flow over an isolated mountain (Case 5, Williamson et 
+# al., 1992)
+###############################################################################
+
+using OrdinaryDiffEq, Trixi, TrixiAtmo
+
+###############################################################################
+# Parameters
+
+initial_condition = initial_condition_isolated_mountain
+polydeg = 5
+cells_per_dimension = (30, 30)
+n_saves = 10
+tspan = (0.0, 15.0 * SECONDS_PER_DAY)
+
+###############################################################################
+# Spatial discretization
+
+mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = polydeg,
+                              element_local_mapping = true)
+
+equations = ShallowWaterEquations3D(gravity = EARTH_GRAVITATIONAL_ACCELERATION,
+                                    rotation_rate = EARTH_ROTATION_RATE)
+
+# Use entropy-conservative two-point flux for volume terms, dissipative surface flux with 
+# simplification for continuous bottom topography
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (FluxPlusDissipation(flux_wintermeyer_etal,
+                                    DissipationLaxFriedrichsEntropyVariables()),
+                flux_nonconservative_wintermeyer_etal)
+
+# Create DG solver with polynomial degree = polydeg
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+# Transform the initial condition to the proper set of conservative variables
+initial_condition_transformed = transform_initial_condition(initial_condition, equations)
+
+# A semidiscretization collects data structures and functions for the spatial 
+# discretization. Here, we pass in the additional keyword argument "auxiliary_field" to 
+# specify the bottom topography.
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver,
+                                    source_terms = source_terms_coriolis_lagrange_multiplier)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0 to T
+ode = semidiscretize(semi, tspan)
+
+# 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. Note that entropy should be conserved at the semi-discrete level.
+analysis_callback = AnalysisCallback(semi, interval = 200,
+                                     save_analysis = true,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (entropy,),
+                                     analysis_polydeg = polydeg)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves,
+                                     solution_variables = cons2prim_and_vorticity)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 0.4)
+
+# 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 = 100.0, save_everystep = false, callback = callbacks)

examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl renamed to examples/elixir_shallowwater_cartesian_unsteady_solid_body_rotation_EC_correction.jl

@@ -6,55 +6,28 @@ using TrixiAtmo
 # Entropy conservation for the spherical shallow water equations in Cartesian
 # form obtained through an entropy correction term
 
-equations = ShallowWaterEquations3D(gravity = 9.81)
+###############################################################################
+# Parameters
 
-# Create DG solver with polynomial degree = 3 and Wintemeyer et al.'s flux as surface flux
+initial_condition = initial_condition_unsteady_solid_body_rotation
 polydeg = 3
-volume_flux = flux_fjordholm_etal #flux_wintermeyer_etal
-surface_flux = flux_fjordholm_etal #flux_wintermeyer_etal  #flux_lax_friedrichs
+cells_per_dimension = (16, 16)
+n_saves = 10
+tspan = (0.0, 15.0 * SECONDS_PER_DAY)
+
+###############################################################################
+# Spatial discretization
+equations = ShallowWaterEquations3D(gravity = EARTH_GRAVITATIONAL_ACCELERATION,
+                                    rotation_rate = EARTH_ROTATION_RATE)
+
+# Create DG solver with polynomial degree = 3 and Fjordholm et al.'s flux as surface flux
+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))
 
-# 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)
-    # Gaussian density
-    rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2))
-
-    # Spherical coordinates for the point x
-    if sign(x[2]) == 0.0
-        signy = 1.0
-    else
-        signy = sign(x[2])
-    end
-    # Co-latitude
-    colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2))
-    # Latitude (auxiliary variable)
-    lat = -colat + 0.5 * pi
-    # Longitude
-    r_xy = sqrt(x[1]^2 + x[2]^2)
-    if r_xy == 0.0
-        phi = pi / 2
-    else
-        phi = signy * acos(x[1] / r_xy)
-    end
-
-    # 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
-    v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha))
-    vlat = -v0 * sin(phi) * sin(alpha)
-
-    # Transform to Cartesian coordinate system
-    v1 = -cos(colat) * cos(phi) * vlat - sin(phi) * v_long
-    v2 = -cos(colat) * sin(phi) * vlat + cos(phi) * v_long
-    v3 = sin(colat) * vlat
-
-    return prim2cons(SVector(rho, v1, v2, v3, 0), equations)
-end
-
 # Source term function to apply the entropy correction term and the Lagrange multiplier to the semi-discretization.
 # The vector contravariant_normal_vector contains the normal contravariant vectors scaled with the inverse Jacobian.
 # The Lagrange multiplier is applied with the vector x, which contains the exact normal direction to the 2D manifold.
@@ -69,27 +42,28 @@ function source_terms_entropy_correction(u, du, x, t,
 
     new_du = SVector(du[1], du[2] + s2, du[3] + s3, du[4] + s4, du[5])
 
-    # Apply Lagrange multipliers using the exact normal vector to the 2D manifold (x)
-    s = source_terms_lagrange_multiplier(u, new_du, x, t, equations, x)
+    # Apply Coriolis and Lagrange multipliers source terms
+    # using the exact normal vector to the 2D manifold (x)
+    s = source_terms_coriolis_lagrange_multiplier(u, new_du, x, t, equations, x)
 
     return SVector(s[1], s[2] + s2, s[3] + s3, s[4] + s4, s[5])
 end
 
-initial_condition = initial_condition_advection_sphere
+# Transform the initial condition to the proper set of conservative variables
+initial_condition_transformed = transform_initial_condition(initial_condition, equations)
 
-mesh = P4estMeshCubedSphere2D(5, 2.0, polydeg = polydeg,
+mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = polydeg,
                               initial_refinement_level = 0)
 
 # A semidiscretization collects data structures and functions for the spatial discretization
-semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver,
                                     metric_terms = MetricTermsInvariantCurl(),
                                     source_terms = source_terms_entropy_correction)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-# Create ODE problem with time span from 0.0 to π
-tspan = (0.0, pi)
+# Create ODE problem
 ode = semidiscretize(semi, tspan)
 
 # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
@@ -100,10 +74,11 @@ summary_callback = SummaryCallback()
 analysis_callback = AnalysisCallback(semi, interval = 10,
                                      save_analysis = true,
                                      extra_analysis_errors = (:conservation_error,),
-                                     extra_analysis_integrals = (waterheight, energy_total))
+                                     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,
+save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves,
                                      solution_variables = cons2prim)
 
 # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step