merge changes from PR #36 (work in progress)

.github/workflows/FormatCheck.yml

@@ -29,7 +29,7 @@ jobs:
         # TODO: Change the call below to
         #       format(".")
         run: |
-          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.60"))'
           julia  -e 'using JuliaFormatter; format(".")'
       - name: Format check
         run: |

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.0-DEV"
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
@@ -17,6 +18,7 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 [compat]
 LinearAlgebra = "1"
 MuladdMacro = "0.2.2"
+LoopVectorization = "0.12.118"
 NLsolve = "4.5.1"
 Reexport = "1.0"
 Static = "0.8, 1"

examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl

@@ -0,0 +1,125 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiAtmo
+###############################################################################
+# Entropy conservation for the spherical shallow water equations in Cartesian
+# form obtained through an entropy correction term
+
+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_fjordholm_etal #flux_wintermeyer_etal
+surface_flux = flux_fjordholm_etal #flux_wintermeyer_etal  #flux_lax_friedrichs
+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))
+    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, 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.
+function source_terms_entropy_correction(u, du, x, t,
+                                         equations::ShallowWaterEquations3D,
+                                         contravariant_normal_vector)
+
+    # Entropy correction term
+    s2 = -contravariant_normal_vector[1] * equations.gravity * u[1]^2
+    s3 = -contravariant_normal_vector[2] * equations.gravity * u[1]^2
+    s4 = -contravariant_normal_vector[3] * equations.gravity * u[1]^2
+
+    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)
+
+    return SVector(s[1], s[2] + s2, s[3] + s3, s[4] + s4, s[5])
+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_entropy_correction)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to π
+tspan = (0.0, pi)
+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 = 10,
+                                     save_analysis = true,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight, energy_total))
+
+# 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()

examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl

@@ -0,0 +1,120 @@
+
+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_fjordholm_etal 
+surface_flux = flux_wintermeyer_etal # flux_fjordholm_etal #flux_lax_friedrichs
+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))
+    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, 0), 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))
+
+# 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()

examples/elixir_shallowwater_cubed_sphere_shell_standard.jl

@@ -0,0 +1,109 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiAtmo
+###############################################################################
+# Entropy consistency test for the spherical shallow water equations in Cartesian
+# form using the standard DGSEM with LLF dissipation
+
+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)
+
+# 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))
+    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, 0), equations)
+end
+
+# Source term function to apply the "standard" Lagrange multiplier to the semi-discretization
+# We call source_terms_lagrange_multiplier, but pass x as the normal direction, as this is the
+# standard way to compute the Lagrange multipliers.
+function source_terms_lagrange_multiplier_standard(u, du, x, t,
+                                                   equations::ShallowWaterEquations3D,
+                                                   contravariant_normal_vector)
+    return source_terms_lagrange_multiplier(u, du, x, t, equations, x)
+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,
+                                    source_terms = source_terms_lagrange_multiplier_standard)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to π
+tspan = (0.0, pi)
+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 = 10,
+                                     save_analysis = true,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight, energy_total))
+
+# 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()