output relative vorticity, add barotropic instability case

Project.toml

@@ -9,6 +9,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718"

examples/elixir_spherical_shallow_water_covariant.jl

@@ -4,46 +4,23 @@
 
 using OrdinaryDiffEq, Trixi, TrixiAtmo
 
-###############################################################################
-# Problem definition
-
-function initial_condition_geostrophic_balance(x, t,
-                                               equations::CovariantShallowWaterEquations2D)
-    (; gravitational_acceleration, rotation_rate) = equations
-
-    radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2)
-    lat = asin(x[3] / radius)
-
-    # compute zonal and meridional components of the velocity
-    V = convert(eltype(x), 2π) * radius / (12 * SECONDS_PER_DAY)
-    v_lon, v_lat = V * cos(lat), zero(eltype(x))
-
-    # compute geopotential height 
-    h = 1 / gravitational_acceleration *
-        (2.94f4 - (radius * rotation_rate * V + 0.5f0 * V^2) * (sin(lat))^2)
-
-    # convert to conservative variables
-    return SVector(h, h * v_lon, h * v_lat)
-end
-
 ###############################################################################
 # Spatial discretization
 
 polydeg = 3
-cells_per_dimension = 2
-
+cells_per_dimension = 5
 element_local_mapping = true
 
-mesh = P4estMeshCubedSphere2D(5, EARTH_RADIUS, polydeg = polydeg,
+mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = polydeg,
                               initial_refinement_level = 0,
                               element_local_mapping = element_local_mapping)
 
 equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION,
                                              EARTH_ROTATION_RATE)
 
-initial_condition = initial_condition_geostrophic_balance
+initial_condition = initial_condition_convergence_test
 source_terms = source_terms_convergence_test
-tspan = (0.0, 12 * SECONDS_PER_DAY)
+tspan = (0.0, 7 * SECONDS_PER_DAY)
 
 # Standard weak-form volume integral
 volume_integral = VolumeIntegralWeakForm()
@@ -67,12 +44,12 @@ ode = semidiscretize(semi, tspan)
 summary_callback = SummaryCallback()
 
 # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
-analysis_callback = AnalysisCallback(semi, interval = 10,
+analysis_callback = AnalysisCallback(semi, interval = 100,
                                      save_analysis = true,
                                      extra_analysis_errors = (:conservation_error,))
 
 # 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]) / 50,
                                      solution_variables = cons2cons)
 
 # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step

examples/elixir_spherical_shallow_water_covariant_barotropic_instability.jl

@@ -0,0 +1,121 @@
+###############################################################################
+# DGSEM for the shallow water equations on the cubed sphere
+###############################################################################
+
+using OrdinaryDiffEq, Trixi, TrixiAtmo, QuadGK
+
+###############################################################################
+# Problem definition
+@inline function galewsky_velocity(θ, u_0, θ_0, θ_1)
+    if (θ_0 < θ) && (θ < θ_1)
+        u = u_0 / exp(-4 / (θ_1 - θ_0)^2) * exp(1 / (θ - θ_0) * 1 / (θ - θ_1))
+    else
+        u = zero(θ)
+    end
+    return u
+end
+
+@inline function galewsky_integrand(θ, u_0, θ_0, θ_1, a,
+                                    equations::CovariantShallowWaterEquations2D)
+    (; rotation_rate) = equations
+    u = galewsky_velocity(θ, u_0, θ_0, θ_1)
+    return u * (2 * rotation_rate * sin(θ) + u * tan(θ) / a)
+end
+
+function initial_condition_barotropic_instability(x, t,
+                                                  equations::CovariantShallowWaterEquations2D)
+    (; gravitational_acceleration, rotation_rate) = equations
+    realT = eltype(x)
+    radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2)
+    lon = atan(x[2], x[1])
+    lat = asin(x[3] / radius)
+
+    # compute zonal and meridional velocity components
+    u_0 = 80.0f0
+    lat_0 = convert(realT, π / 7)
+    lat_1 = convert(realT, π / 2) - lat_0
+    v_lon = galewsky_velocity(lat, u_0, lat_0, lat_1)
+    v_lat = zero(eltype(x))
+
+    # numerically integrate (here we use the QuadGK package) to get height
+    galewsky_integral, _ = quadgk(latp -> galewsky_integrand(latp, u_0, lat_0, lat_1,
+                                                             radius,
+                                                             equations), -π / 2, lat)
+    h = 10158.0f0 - radius / gravitational_acceleration * galewsky_integral
+
+    # add perturbation to initiate instability
+    α = convert(realT, 1 / 3)
+    β = convert(realT, 1 / 15)
+    lat_2 = convert(realT, π / 4)
+    if (-π < lon) && (lon < π)
+        h = h + 120.0f0 * cos(lat) * exp(-((lon / α)^2)) * exp(-((lat_2 - lat) / β)^2)
+    end
+    # convert to conservative variables
+    return SVector(h, h * v_lon, h * v_lat)
+end
+
+###############################################################################
+# Spatial discretization
+
+polydeg = 5
+cells_per_dimension = 5
+element_local_mapping = true
+
+mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = polydeg,
+                              initial_refinement_level = 0,
+                              element_local_mapping = element_local_mapping)
+
+equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION,
+                                             EARTH_ROTATION_RATE)
+
+initial_condition = initial_condition_barotropic_instability
+source_terms = source_terms_convergence_test
+tspan = (0.0, 6 * SECONDS_PER_DAY)
+
+# Standard weak-form volume integral
+volume_integral = VolumeIntegralWeakForm()
+
+# Create DG solver with polynomial degree = p and a local Lax-Friedrichs flux
+solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs,
+               volume_integral = volume_integral)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms)
+
+###############################################################################
+# 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
+analysis_callback = AnalysisCallback(semi, interval = 100,
+                                     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]) / 50,
+                                     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
+tsave = LinRange(tspan..., 100)
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 100.0, save_everystep = false, saveat = tsave, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

examples/elixir_spherical_shallow_water_covariant_rossby_haurwitz.jl

@@ -0,0 +1,107 @@
+###############################################################################
+# DGSEM for the shallow water equations on the cubed sphere
+###############################################################################
+
+using OrdinaryDiffEq, Trixi, TrixiAtmo
+
+###############################################################################
+# Problem definition
+
+function initial_condition_rossby_haurwitz(x, t,
+                                           equations::CovariantShallowWaterEquations2D)
+    (; gravitational_acceleration, rotation_rate) = equations
+
+    radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2)
+    lon = atan(x[2], x[1])
+    lat = asin(x[3] / radius)
+
+    h_0 = 8.0f3
+    K = 7.848f-6
+    R = 4.0f0
+
+    A = 0.5f0 * K * (2 * rotation_rate + K) * (cos(lat))^2 +
+        0.25f0 * K^2 * (cos(lat))^(2 * R) *
+        ((R + 1) * (cos(lat))^2 +
+         (2 * R^2 - R - 2) - 2 * R^2 / ((cos(lat))^2))
+    B = 2 * (rotation_rate + K) * K / ((R + 1) * (R + 2)) * (cos(lat))^R *
+        ((R^2 + 2R + 2) - (R + 1)^2 * (cos(lat))^2)
+    C = 0.25f0 * K^2 * (cos(lat))^(2 * R) * ((R + 1) * (cos(lat))^2 - (R + 2))
+
+    h = h_0 +
+        (1 / gravitational_acceleration) *
+        (radius^2 * A + radius^2 * B * cos(R * lon) + radius^2 * C * cos(2 * R * lon))
+
+    v_lon = radius * K * cos(lat) +
+            radius * K * (cos(lat))^(R - 1) * (R * (sin(lat))^2 - (cos(lat))^2) *
+            cos(R * lon)
+    v_lat = -radius * K * R * (cos(lat))^(R - 1) * sin(lat) * sin(R * lon)
+
+    # convert to conservative variables
+    return SVector(h, h * v_lon, h * v_lat)
+end
+
+###############################################################################
+# Spatial discretization
+
+polydeg = 5
+cells_per_dimension = 5
+element_local_mapping = true
+
+mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = polydeg,
+                              initial_refinement_level = 0,
+                              element_local_mapping = element_local_mapping)
+
+equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION,
+                                             EARTH_ROTATION_RATE)
+
+initial_condition = initial_condition_rossby_haurwitz
+source_terms = source_terms_convergence_test
+tspan = (0.0, 7 * SECONDS_PER_DAY)
+
+# Standard weak-form volume integral
+volume_integral = VolumeIntegralWeakForm()
+
+# Create DG solver with polynomial degree = p and a local Lax-Friedrichs flux
+solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs,
+               volume_integral = volume_integral)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms)
+
+###############################################################################
+# 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
+analysis_callback = AnalysisCallback(semi, interval = 100,
+                                     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]) / 50,
+                                     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
+tsave = LinRange(tspan..., 100)
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 100.0, save_everystep = false, saveat = tsave, callback = callbacks);
+
+# Print the timer summary
+summary_callback()