add NDIMS_AMBIENT parameter to equations

src/equations/covariant_advection.jl

@@ -1,16 +1,20 @@
 @muladd begin
 #! format: noindent
 
-# Variable-coefficient linear advection equation in covariant form
-struct CovariantLinearAdvectionEquation2D <: AbstractCovariantEquations{2, 3} end
+# Variable-coefficient linear advection equation in covariant form on a two-dimensional
+# surface in three-dimensional space
+struct CovariantLinearAdvectionEquation2D <: AbstractCovariantEquations{2, 3, 3} end
 
 # The first variable is the scalar conserved quantity, and the second two are the 
 # contravariant velocity components, which are spatially varying but constant in time.
 function Trixi.varnames(::typeof(cons2cons), ::CovariantLinearAdvectionEquation2D)
     return ("scalar", "v_con_1", "v_con_2")
 end
 
-Trixi.cons2entropy(u, ::CovariantLinearAdvectionEquation2D) = u
+# We will measure the entropy/energy simply as half of the square of the scalar variable
+Trixi.cons2entropy(u, ::CovariantLinearAdvectionEquation2D) = SVector{3}(u[1],
+                                                                         zero(eltype(u)),
+                                                                         zero(eltype(u)))
 
 # The flux for the covariant form takes in the element container and node/element indices
 # in order to give the flux access to the geometric information

src/equations/equations.jl

@@ -7,8 +7,9 @@ const EARTH_ROTATION_RATE = 7.292f-5  # rad/s
 const SECONDS_PER_DAY = 8.64f4
 
 # Abstract type used to dispatch specialized solvers for the covariant form
-abstract type AbstractCovariantEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, 
-                                                                            NVARS} end
+abstract type AbstractCovariantEquations{NDIMS,
+                                         NDIMS_AMBIENT,
+                                         NVARS} <: AbstractEquations{NDIMS, NVARS} end
 
 # Numerical flux plus dissipation which passes node/element indices and cache. 
 # We assume that u_ll and u_rr have been transformed into the same local coordinate system.

src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl

@@ -31,11 +31,16 @@ end
     return uEltype
 end
 
-@inline function Trixi.get_node_coords(x, ::AbstractCovariantEquations{2}, ::DG,
-                                       indices...)
-    return SVector(ntuple(@inline(idx->x[idx, indices...]), 3))
+# Get Cartesian node positions, dispatching on the dimension of the equation system
+@inline function Trixi.get_node_coords(x,
+                                       ::AbstractCovariantEquations{NDIMS,
+                                                                    NDIMS_AMBIENT},
+                                       ::DG,
+                                       indices...) where {NDIMS, NDIMS_AMBIENT}
+    return SVector(ntuple(@inline(idx->x[idx, indices...]), NDIMS_AMBIENT))
 end
 
+# Volume element J = sqrt(det(G)), where G is the matrix of covariant metric components
 @inline function volume_element(elements, i, j, element)
     return 1 / elements.inverse_jacobian[i, j, element]
 end
@@ -58,14 +63,17 @@ end
            Jv_con_2 * elements.inverse_jacobian[i, j, element]
 end
 
-# Create element container and initialize element data.
-# This function dispatches on the dimensions of the mesh and the equation type 
-function Trixi.init_elements(mesh::P4estMesh{NDIMS},
-                             equations::AbstractCovariantEquations{NDIMS},
-                             basis, ::Type{uEltype}) where {NDIMS, uEltype <: Real}
-    RealT = real(mesh)
-    NDIMS_AMBIENT = size(mesh.tree_node_coordinates, 1)
-
+# Create element container and initialize element data for a mesh of dimension NDIMS in 
+# an ambient space of dimension NDIMS_AMBIENT, with the equations expressed in covariant 
+# form
+function Trixi.init_elements(mesh::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT},
+                             equations::AbstractCovariantEquations{NDIMS,
+                                                                   NDIMS_AMBIENT},
+                             basis,
+                             ::Type{uEltype}) where {NDIMS,
+                                                     NDIMS_AMBIENT,
+                                                     RealT <: Real,
+                                                     uEltype <: Real}
     nelements = Trixi.ncells(mesh)
 
     _node_coordinates = Vector{RealT}(undef,
@@ -116,7 +124,7 @@ end
 # Compute the node positions and metric terms for the covariant form, assuming that the
 # domain is a spherical shell. We do not make any assumptions on the mesh topology.
 function Trixi.init_elements!(elements, mesh::P4estMesh{2, 3},
-                              equations::AbstractCovariantEquations{2},
+                              equations::AbstractCovariantEquations{2, 3},
                               basis::LobattoLegendreBasis)
     (; node_coordinates, covariant_basis, inverse_jacobian) = elements
 

src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl

@@ -3,8 +3,8 @@
 
 # Custom rhs! for the covariant form. Note that the mortar flux is not included, as we do
 # not yet support non-conforming meshes for the covariant solver.
-function Trixi.rhs!(du, u, t,
-                    mesh::P4estMesh{2}, equations::AbstractCovariantEquations{2},
+function Trixi.rhs!(du, u, t, mesh::P4estMesh{2},
+                    equations::AbstractCovariantEquations{2},
                     initial_condition, boundary_conditions, source_terms::Source,
                     dg::DG, cache) where {Source}
     # Reset du
@@ -65,12 +65,13 @@ end
 
 # Evaluate the initial condition in spherical vector components, then 
 # transform to contravariant components for use as prognostic variables
-function Trixi.compute_coefficients!(u, func, t, mesh::Trixi.AbstractMesh{2},
-                                     equations::AbstractCovariantEquations{2}, dg::DG,
+function Trixi.compute_coefficients!(u, func, t, mesh::P4estMesh,
+                                     equations::AbstractCovariantEquations, dg::DG,
                                      cache)
     for element in eachelement(dg, cache)
         for j in eachnode(dg), i in eachnode(dg)
-            x_node = Trixi.get_node_coords(node_coordinates, equations, dg, i, j, element)
+            x_node = Trixi.get_node_coords(cache.elements.node_coordinates,
+                                           equations, dg, i, j, element)
             u_spherical = func(x_node, t, equations)
             u_con = spherical2contravariant(u_spherical, equations,
                                             cache.elements, i, j, element)
@@ -149,11 +150,23 @@ end
     return nothing
 end
 
+# Outward unit normal vector in reference coordinates for a quadrilateral element
+@inline function reference_normal_vector(direction)
+    orientation = (direction + 1) >> 1
+    sign = isodd(direction) ? -1 : 1
+
+    if orientation == 1
+        return SVector(sign, 0.0f0, 0.0f0)
+    else
+        return SVector(0.0f0, sign, 0.0f0)
+    end
+end
+
 # Interface flux which transforms the contravariant prognostic variables into the same 
 # reference coordinate system on each side of the interface, then applies the numerical 
 # flux in reference space, taking in the reference normal vector
 function Trixi.calc_interface_flux!(surface_flux_values,
-                                    mesh::Union{P4estMesh{2}, T8codeMesh{2}},
+                                    mesh::P4estMesh{2},
                                     nonconservative_terms,
                                     equations::AbstractCovariantEquations{2},
                                     surface_integral, dg::DG, cache)