Quad-Based Icosahedral Grid (#50)

* First implementation of an icosahedral grid with 3 quads on each triangular face

* Removed cells_per_dimension from elixir

* Changed radius in initial_condition_gaussian

* Added test for spherical advection on icosahedral grid

* Improved ascii art?

* Added information to run convergence test with P4estMeshCubedSphere2D

* Corrected advection elixir for tests

* Fixed initial condition to avoid pole singularities

* WIP: Icosahedral grid with multiple trees per quad face of icosahedron (#52)

* Added the possibility to have multiple trees per quad face of icosahedron

* Updated advection elixir (and test) that uses the quad-based icosahedral grid

* Updated spherical advection reference solutions to account for changes in the initial condition

* Improve comments

* Remove sqrt in initial condition to make it slightly cheaper

* Fix comment

* Redefined index maps as constant tuples and improved the comments

* typo

* Moved all export statements to TrixiAtmo.jl and fixed typo

examples/elixir_shallowwater_cubed_sphere_shell_advection.jl

@@ -3,11 +3,21 @@ using OrdinaryDiffEq
 using Trixi
 using TrixiAtmo
 
+# To run a convergence test, we have two options:
+# 1. Use the p4est variable initial_refinement_level to refine the grid:
+#    - To do this, line 46 ("initial_refinement_level = 0") must NOT be a comment
+#    - Call convergence_test("../examples/elixir_shallowwater_cubed_sphere_shell_advection.jl", 4, initial_refinement_level = 0)
+#    - NOT OPTIMAL: Good convergence the first iterations, but then it stagnates. Reason: The geometry does not improve with refinement.
+# 2. Use the variable trees_per_face_dimension of P4estMeshCubedSphere2D
+#    - To do this, line 46 ("initial_refinement_level = 0") MUST BE commented/removed.
+#    - Call convergence_test("../examples/elixir_shallowwater_cubed_sphere_shell_advection.jl", 4, cells_per_dimension = (3,3))
+#    - OPTIMAL convergence of polydeg + 1. Reason: The geometry improves with refinement.
+
 ###############################################################################
 # semidiscretization of the linear advection equation
 
-cells_per_dimension = 5
 initial_condition = initial_condition_gaussian
+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
@@ -32,8 +42,8 @@ function source_terms_convert_to_linear_advection(u, du, x, t,
 end
 
 # Create a 2D cubed sphere mesh the size of the Earth
-mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = 3,
-                              initial_refinement_level = 0,
+mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = 3,
+                              #initial_refinement_level = 0, # Comment to use cells_per_dimension in the convergence test
                               element_local_mapping = false)
 
 # Convert initial condition given in terms of zonal and meridional velocity components to 
@@ -54,12 +64,12 @@ ode = semidiscretize(semi, (0.0, 12 * SECONDS_PER_DAY))
 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(interval = 100,
                                      solution_variables = cons2prim)
 
 # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step

examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl

@@ -0,0 +1,91 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiAtmo
+
+# To run a convergence test, we have two options:
+# 1. Use the p4est variable initial_refinement_level to refine the grid:
+#    - To do this, line 46 ("initial_refinement_level = 0") must NOT be a comment
+#    - Call convergence_test("../examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl", 4, initial_refinement_level = 0)
+#    - NOT OPTIMAL: Good convergence the first iterations, but then it stagnates. Reason: The geometry does not improve with refinement.
+# 2. Use the variable trees_per_face_dimension of P4estMeshQuadIcosahedron2D
+#    - To do this, line 46 ("initial_refinement_level = 0") MUST BE commented/removed
+#    - Call convergence_test("../examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl", 4, cells_per_dimension = (1,1))
+#    - OPTIMAL convergence of polydeg + 1. Reason: The geometry improves with refinement.
+
+###############################################################################
+# semidiscretization of the linear advection equation
+initial_condition = initial_condition_gaussian
+polydeg = 3
+cells_per_dimension = (2, 2)
+
+# 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 = polydeg, surface_flux = flux_lax_friedrichs)
+
+# 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,
+                                                  equations::ShallowWaterEquations3D,
+                                                  normal_direction)
+    v1 = u[2] / u[1]
+    v2 = u[3] / u[1]
+    v3 = u[4] / u[1]
+
+    s2 = du[1] * v1 - du[2]
+    s3 = du[1] * v2 - du[3]
+    s4 = du[1] * v3 - du[4]
+
+    return SVector(0.0, s2, s3, s4, 0.0)
+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,
+                                  polydeg = polydeg)
+
+# Convert initial condition given in terms of zonal and meridional velocity components to 
+# one given in terms of Cartesian momentum components
+initial_condition_transformed = transform_to_cartesian(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_convert_to_linear_advection)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+ode = semidiscretize(semi, (0.0, 12 * SECONDS_PER_DAY))
+
+# 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(interval = 100,
+                                     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/TrixiAtmo.jl

@@ -37,7 +37,8 @@ export flux_chandrashekar, flux_LMARS
 export velocity, waterheight, pressure, energy_total, energy_kinetic, energy_internal,
        lake_at_rest_error, source_terms_lagrange_multiplier,
        clean_solution_lagrange_multiplier!
-export P4estCubedSphere2D, MetricTermsCrossProduct, MetricTermsInvariantCurl
+export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, MetricTermsCrossProduct,
+       MetricTermsInvariantCurl
 export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION,
        EARTH_ROTATION_RATE, SECONDS_PER_DAY
 export spherical2contravariant, contravariant2spherical, spherical2cartesian,

src/equations/reference_data.jl

@@ -52,7 +52,7 @@ for Cartesian and covariant formulations.
 @inline function initial_condition_gaussian(x, t)
     RealT = eltype(x)
 
-    a = EARTH_RADIUS  # radius of the sphere in metres
+    a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) #radius of the sphere
     V = convert(RealT, 2π) * a / (12 * SECONDS_PER_DAY)  # speed of rotation
     alpha = convert(RealT, π / 4)  # angle of rotation
     h_0 = 1000.0f0  # bump height in metres
@@ -68,7 +68,14 @@ for Cartesian and covariant formulations.
     h = h_0 * exp(-b_0 * ((x[1] - x_0[1])^2 + (x[2] - x_0[2])^2 + (x[3] - x_0[3])^2))
 
     # get zonal and meridional components of the velocity
-    lon, lat = atan(x[2], x[1]), asin(x[3] / a)
+    r_xy = x[1]^2 + x[2]^2
+    if r_xy == 0.0
+        lon = pi / 2
+    else
+        lon = atan(x[2], x[1])
+    end
+
+    lat = asin(x[3] / a)
     vlon = V * (cos(lat) * cos(alpha) + sin(lat) * cos(lon) * sin(alpha))
     vlat = -V * sin(lon) * sin(alpha)
 

src/meshes/meshes.jl

@@ -1,2 +1,2 @@
-export P4estMeshCubedSphere2D
 include("p4est_cubed_sphere.jl")
+include("p4est_icosahedron_quads.jl")