p4est solver on a spherical mesh (#25)

* Added containers and 2D p4est mesh from my spherical shell implementation in Trixi.jl (cherry-picked from f45378e)

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

* Added element container with PtrArray for performance and modified the structure

* format

* Fixed bug in the definition of init_elements and added nelements(). eltype() and ndims() functions for new struct

* Apply suggestions

* Added Cartesian weak-form kernel for 2D manifolds in 3D

* format

---------

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

Project.toml

@@ -1,6 +1,6 @@
 name = "TrixiAtmo"
 uuid = "c9ed1054-d170-44a9-8ee2-d5566f5d1389"
-authors = ["Benedict Geihe <[email protected]>", "Tristan Montoya <[email protected]>", "Hendrik Ranocha <[email protected]>", "Michael Schlottke-Lakemper <[email protected]>"]
+authors = ["Benedict Geihe <[email protected]>", "Tristan Montoya <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrés Rueda-Ramírez <[email protected]>", "Michael Schlottke-Lakemper <[email protected]>"]
 version = "0.1.0-DEV"
 
 [deps]
@@ -9,6 +9,8 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718"
+StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
 Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 
 [compat]
@@ -17,5 +19,7 @@ MuladdMacro = "0.2.2"
 Reexport = "1.0"
 Static = "0.8, 1"
 StaticArrays = "1"
+StaticArrayInterface = "1.5.1"
+StrideArrays = "0.1.28"
 Trixi = "0.7, 0.8"
 julia = "1.9"

examples/p4est_2d_dgsem/elixir_euler_spherical_advection_cartesian.jl

@@ -0,0 +1,151 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiAtmo
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+equations = CompressibleEulerEquations3D(1.4)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+polydeg = 3
+solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs)
+
+# 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::CompressibleEulerEquations3D)
+    # Gaussian density
+    rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2))
+    # Constant pressure
+    p = 1.0
+
+    # 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))
+    v_lat = -v0 * sin(phi) * sin(alpha)
+
+    # Transform to Cartesian coordinate system
+    v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long
+    v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long
+    v3 = sin(colat) * v_lat
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+
+# Source term function to transform the Euler equations into the linear advection equations with variable advection velocity
+function source_terms_convert_to_linear_advection(u, du, x, t,
+                                                  equations::CompressibleEulerEquations3D)
+    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]
+    s5 = 0.5 * (s2 * v1 + s3 * v2 + s4 * v3) - du[5]
+
+    return SVector(0.0, s2, s3, s4, s5)
+end
+
+# Custom RHS that applies a source term that depends on du (used to convert the 3D Euler equations into the 3D linear advection
+# equations with position-dependent advection speed)
+function rhs_semi_custom!(du_ode, u_ode, semi, t)
+    # Compute standard Trixi RHS
+    Trixi.rhs!(du_ode, u_ode, semi, t)
+
+    # Now apply the custom source term
+    Trixi.@trixi_timeit Trixi.timer() "custom source term" begin
+        @unpack solver, equations, cache = semi
+        @unpack node_coordinates = cache.elements
+
+        # Wrap the solution and RHS
+        u = Trixi.wrap_array(u_ode, semi)
+        du = Trixi.wrap_array(du_ode, semi)
+
+        Trixi.@threaded for element in eachelement(solver, cache)
+            for j in eachnode(solver), i in eachnode(solver)
+                u_local = Trixi.get_node_vars(u, equations, solver, i, j, element)
+                du_local = Trixi.get_node_vars(du, equations, solver, i, j, element)
+                x_local = Trixi.get_node_coords(node_coordinates, equations, solver,
+                                                i, j, element)
+                source = source_terms_convert_to_linear_advection(u_local, du_local,
+                                                                  x_local, t, equations)
+                Trixi.add_to_node_vars!(du, source, equations, solver, i, j, element)
+            end
+        end
+    end
+end
+
+initial_condition = initial_condition_advection_sphere
+
+mesh = TrixiAtmo.P4estMeshCubedSphere2D(5, 1.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)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to π
+tspan = (0.0, pi)
+ode = semidiscretize(semi, tspan)
+
+# Create custom discretization that runs with the custom RHS
+ode_semi_custom = ODEProblem(rhs_semi_custom!,
+                             ode.u0,
+                             ode.tspan,
+                             semi)
+
+# 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 = (Trixi.density,))
+
+# 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_semi_custom, 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

@@ -12,12 +12,17 @@ using Trixi
 using MuladdMacro: @muladd
 using StaticArrays: SVector
 using Static: True, False
+using StrideArrays: PtrArray
+using StaticArrayInterface: static_size
 using LinearAlgebra: norm
 using Reexport: @reexport
 @reexport using StaticArrays: SVector
 
 include("auxiliary/auxiliary.jl")
 include("equations/equations.jl")
+include("meshes/meshes.jl")
+include("solvers/solvers.jl")
+include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl")
 
 export CompressibleMoistEulerEquations2D
 

src/meshes/meshes.jl

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