Merge remote-tracking branch 'origin/tm/test_spherical_advection' int…

…o tm/covariant_form
trixi-framework · Aug 20, 2024 · 2fd7610 · 2fd7610
2 parents 73dda27 + 317033a
commit 2fd7610
Showing 1 changed file with 13 additions and 6 deletions.
diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl
@@ -151,7 +151,7 @@ function calc_node_coordinates_2d_shell!(node_coordinates,
     # Hanging nodes will cause holes in the mesh if its polydeg is higher
     # than the polydeg of the solver.
     @assert length(basis.nodes)>=length(mesh.nodes) "The solver can't have a lower polydeg than the mesh"
-
+    @assert size(mesh.tree_node_coordinates, 1)==3 "Shell must use 3D node coordinates"
     calc_node_coordinates_2d_shell!(node_coordinates, mesh, basis.nodes)
 end
 
@@ -163,7 +163,7 @@ function calc_node_coordinates_2d_shell!(node_coordinates,
     # places and the additional information passed to the compiler makes them faster
     # than native `Array`s.
     tmp1 = Trixi.StrideArray(undef, real(mesh),
-                             Trixi.StaticInt(size(mesh.tree_node_coordinates, 1)),
+                             Trixi.StaticInt(3),
                              Trixi.static_length(nodes),
                              Trixi.static_length(mesh.nodes))
     matrix1 = Trixi.StrideArray(undef, real(mesh),
@@ -174,7 +174,6 @@ function calc_node_coordinates_2d_shell!(node_coordinates,
 
     # Macros from `p4est`
     p4est_root_len = 1 << Trixi.P4EST_MAXLEVEL
-    p4est_root_len = 1 << Trixi.P4EST_MAXLEVEL
     p4est_quadrant_len(l) = 1 << (Trixi.P4EST_MAXLEVEL - l)
 
     trees = Trixi.unsafe_wrap_sc(Trixi.p4est_tree_t, mesh.p4est.trees)
@@ -209,6 +208,7 @@ function calc_node_coordinates_2d_shell!(node_coordinates,
     return node_coordinates
 end
 
+# This only works for a sphere
 function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArray{<:Any,
                                                                                    5},
                                               element,
@@ -217,6 +217,10 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
     @unpack derivative_matrix = basis
 
     for j in eachnode(basis), i in eachnode(basis)
+        radius = sqrt(node_coordinates[1, i, j, element]^2 +
+                      node_coordinates[2, i, j, element]^2 +
+                      node_coordinates[3, i, j, element]^2)
+
         for n in 1:3
             # (n, m, l) cyclic
             m = (n % 3) + 1
@@ -230,7 +234,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
                                                           jacobian_matrix[l, 2, i, j,
                                                                           element] *
                                                           node_coordinates[m, i, j,
-                                                                           element])
+                                                                           element]) /
+                                                         radius
 
             contravariant_vectors[n, 2, i, j, element] = (jacobian_matrix[l, 1, i, j,
                                                                           element] *
@@ -240,7 +245,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
                                                           jacobian_matrix[m, 1, i, j,
                                                                           element] *
                                                           node_coordinates[l, i, j,
-                                                                           element])
+                                                                           element]) /
+                                                         radius
 
             contravariant_vectors[n, 3, i, j, element] = (jacobian_matrix[m, 1, i, j,
                                                                           element] *
@@ -250,7 +256,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
                                                           jacobian_matrix[m, 2, i, j,
                                                                           element] *
                                                           jacobian_matrix[l, 1, i, j,
-                                                                          element])
+                                                                          element]) /
+                                                         radius
         end
     end