Covariant form as variable-coefficient problem using auxiliary variab…

…les (#47) * linear advection on the sphere * formatter * very preliminary implementation of covariant form with p4est * removed untracked * removed allocations, now 100x faster * add DS_Store to gitignore * add DS_Store to gitignore * made equations actually 2D * integrated upstream changes * analysis callback for covariant form * cleanup * separate Andres's spherical shell containers from main Trixi * 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 * integrated Andres's new container type with covariant solver * added flux-differencing kernel for covariant form * integrated Andrés' PR with covariant solver * add tests for spherical advection in Cartesian coords * revert change to moist bubble case * revert change to moist bubble case * covariant weak form with tests * removed shallow water from this branch * reorganize a bit * hard-code dimension of node coordinates to avoid type instability * metric terms for shell of radius != 1 * test for covariant form is now on earth scale * fix name of cartesian test * moved i,j,element,cache into flux signatures for covariant form * added test for flux differencing covariant form * transformation between local coordinate systems now uses latitude-longitude as intermediate state, not Cartesian * format * make compatible with my PR trixi-framework/Trixi.jl#2068 * separate containers to store geometric information specific to covariant form * add element-local mapping from Guba et al. * add reference to Guba et al. to docstring * new mapping works with cartesian solver solver in covariant_form branch * exact metrics for element-local mapping * streamline tests and improve readability * rename default DGSEM mapping * didn't need new analysis cache * improve comments * refactor * remove DiffEqCallbacks for now * add better comments * improve comments * better documentation * add NDIMS_AMBIENT parameter to equations * fix docstring format * make link in docs work * fix link in docs * documentation for the covariant solver * added surface central flux and test * fix trixi_compat * specialize compute_coefficients! to dimension 2 * format * fix spelling * remove initial_condition parameter from rhs * no need to use Float32 for physical constant defaults * incorporate downstream changes * format * standardize test cases * move covariant form max_dt method to callbacks_step * Fix comments * generalize interface between equations and cache * relax relative tolerance on cartesian advection test) * improve docs * streamline initial condition definition and improve docs * tests for cartesian form now use standard and element-local mapping approaches * first attempt at aux vars * works with allocations * fix allocations * clean up docs * format and update comments * move auxvars into own cache field * improve comments * format comments * revert back to without basis matrix inverse * add warning for P4estMeshCubedSphere and rename volume->area element * rename function transforming initial condition * fix elixir and docs * poldeg(solver) -> Trixi.polydeg(solver) * reformat ugly line breaks * modularize covariant basis calculation for easier testing * remove duplicated cons2entropy function --------- Co-authored-by: Andrés Rueda-Ramírez <[email protected]>
trixi-framework · Nov 25, 2024 · 40800a6 · 40800a6
1 parent 49ec206
commit 40800a6
Show file tree

Hide file tree

Showing 28 changed files with 1,197 additions and 166 deletions.
diff --git a/Project.toml b/Project.toml
@@ -5,8 +5,8 @@ version = "0.1.0-DEV"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
@@ -17,13 +17,13 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 
 [compat]
 LinearAlgebra = "1"
-MuladdMacro = "0.2.2"
 LoopVectorization = "0.12.118"
+MuladdMacro = "0.2.2"
 NLsolve = "4.5.1"
 Reexport = "1.0"
 Static = "0.8, 1"
 StaticArrayInterface = "1.5.1"
 StaticArrays = "1"
 StrideArrays = "0.1.28"
-Trixi = "0.7, 0.8, 0.9"
+Trixi = "0.9"
 julia = "1.9"
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -1,7 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
-TrixiAtmo = "c9ed1054-d170-44a9-8ee2-d5566f5d1389"
 
 [compat]
 Documenter = "1.3"

diff --git a/docs/make.jl b/docs/make.jl
@@ -1,7 +1,14 @@
-using TrixiAtmo
 using Documenter
 using DocumenterInterLinks
 
+# Fix for https://github.com/trixi-framework/Trixi.jl/issues/668
+if (get(ENV, "CI", nothing) != "true") &&
+   (get(ENV, "TRIXI_DOC_DEFAULT_ENVIRONMENT", nothing) != "true")
+    push!(LOAD_PATH, dirname(@__DIR__))
+end
+
+using TrixiAtmo
+
 # Provide external links to the Trixi.jl docs (project root and inventory file)
 links = InterLinks("Trixi" => ("https://trixi-framework.github.io/Trixi.jl/stable/",
                                "https://trixi-framework.github.io/Trixi.jl/stable/objects.inv"))

diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl b/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl
@@ -45,12 +45,12 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio
     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)
+    vlat = -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
+    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

diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl b/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl
@@ -47,12 +47,12 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio
     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)
+    vlat = -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
+    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

diff --git a/...ir_euler_spherical_advection_cartesian.jl → ...llowwater_cubed_sphere_shell_advection.jl b/...ir_euler_spherical_advection_cartesian.jl → ...llowwater_cubed_sphere_shell_advection.jl
@@ -6,58 +6,19 @@ using TrixiAtmo
 ###############################################################################
 # semidiscretization of the linear advection equation
 
-# We use the Euler equations structure but modify the rhs! function to convert it to a
-# variable-coefficient advection equation
-equations = CompressibleEulerEquations3D(1.4)
+cells_per_dimension = 5
+initial_condition = initial_condition_gaussian
+
+# 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
-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
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
-# Source term function to transform the Euler equations into the linear advection equations with variable advection velocity
+# 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::CompressibleEulerEquations3D,
+                                                  equations::ShallowWaterEquations3D,
                                                   normal_direction)
     v1 = u[2] / u[1]
     v2 = u[3] / u[1]
@@ -66,29 +27,27 @@ function source_terms_convert_to_linear_advection(u, du, x, t,
     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)
+    return SVector(0.0, s2, s3, s4, 0.0)
 end
 
-initial_condition = initial_condition_advection_sphere
-
-element_local_mapping = false
-
-mesh = P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg,
+# Create a 2D cubed sphere mesh the size of the Earth
+mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = 3,
                               initial_refinement_level = 0,
-                              element_local_mapping = element_local_mapping)
+                              element_local_mapping = false)
+
+# 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, solver,
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver,
                                     source_terms = source_terms_convert_to_linear_advection)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-# Create ODE problem with time span from 0.0 to π
-tspan = (0.0, pi)
-ode = semidiscretize(semi, tspan)
+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
@@ -97,8 +56,7 @@ 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,))
+                                     extra_analysis_errors = (:conservation_error,))
 
 # The SaveSolutionCallback allows to save the solution to a file in regular intervals
 save_solution = SaveSolutionCallback(interval = 10,

diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl b/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl
@@ -41,12 +41,12 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio
     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)
+    vlat = -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
+    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

diff --git a/examples/elixir_spherical_advection_covariant.jl b/examples/elixir_spherical_advection_covariant.jl
@@ -0,0 +1,68 @@
+###############################################################################
+# DGSEM for the linear advection equation on the cubed sphere
+###############################################################################
+
+using OrdinaryDiffEq, Trixi, TrixiAtmo
+
+###############################################################################
+# Spatial discretization
+
+cells_per_dimension = 5
+initial_condition = initial_condition_gaussian
+
+equations = CovariantLinearAdvectionEquation2D()
+
+# Create DG solver with polynomial degree = p and a local Lax-Friedrichs flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs,
+               volume_integral = VolumeIntegralWeakForm())
+
+# Create a 2D cubed sphere mesh the size of the Earth. For the covariant form to work 
+# properly, we currently need polydeg to equal that of the solver, 
+# initial_refinement_level = 0, and element_local_mapping = true.
+mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS,
+                              polydeg = Trixi.polydeg(solver),
+                              initial_refinement_level = 0,
+                              element_local_mapping = true)
+
+# Convert initial condition given in terms of zonal and meridional velocity components to 
+# one given in terms of contravariant velocity components
+initial_condition_transformed = transform_to_contravariant(initial_condition, equations)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0 to T
+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 = 10,
+                                     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,
+                                     solution_variables = cons2cons)
+
+# 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, save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()
diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl
@@ -8,32 +8,40 @@ See also: [trixi-framework/TrixiAtmo.jl](https://github.com/trixi-framework/Trix
 """
 module TrixiAtmo
 
+using Reexport: @reexport
 using Trixi
 using MuladdMacro: @muladd
 using Static: True, False
 using StrideArrays: PtrArray
 using StaticArrayInterface: static_size
-using LinearAlgebra: norm
+using LinearAlgebra: norm, dot
 using Reexport: @reexport
 using LoopVectorization: @turbo
-@reexport using StaticArrays: SVector
+
+@reexport using StaticArrays: SVector, SMatrix
 
 include("auxiliary/auxiliary.jl")
 include("equations/equations.jl")
 include("meshes/meshes.jl")
 include("semidiscretization/semidiscretization.jl")
 include("solvers/solvers.jl")
 include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl")
-include("callbacks_step/stepsize_dg2d.jl")
+include("callbacks_step/callbacks_step.jl")
 
-export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D
+export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D,
+       CovariantLinearAdvectionEquation2D
 
 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 examples_dir
+export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION,
+       EARTH_ROTATION_RATE, SECONDS_PER_DAY
+export spherical2contravariant, contravariant2spherical, spherical2cartesian,
+       transform_to_cartesian, transform_to_contravariant
+export initial_condition_gaussian
 
+export examples_dir
 end # module TrixiAtmo
diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl
@@ -7,6 +7,7 @@ read-only and should *not* be modified.
 To find out which files are available, use, e.g., `readdir`:
 
 # Examples
+
 ```@example
 readdir(examples_dir())
 ```