wip

Author: nefrathenrici

src/prognostic_equations/implicit/manual_sparse_jacobian.jl

@@ -71,12 +71,11 @@ function jacobian_cache(alg::ManualSparseJacobian, Y, atmos)
         approximate_solve_iters,
     ) = alg
     FT = Spaces.undertype(axes(Y.c))
-    CTh = CTh_vector_type(axes(Y.c))
 
     DiagonalRow = DiagonalMatrixRow{FT}
     TridiagonalRow = TridiagonalMatrixRow{FT}
     BidiagonalRow_C3 = BidiagonalMatrixRow{C3{FT}}
-    TridiagonalRow_ACTh = TridiagonalMatrixRow{Adjoint{FT, CTh{FT}}}
+    TridiagonalRow_ACT12 = TridiagonalMatrixRow{Adjoint{FT, CT12{FT}}}
     BidiagonalRow_ACT3 = BidiagonalMatrixRow{Adjoint{FT, CT3{FT}}}
     BidiagonalRow_C3xACT12 =
         BidiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CT12{FT})')}
@@ -150,7 +149,7 @@ function jacobian_cache(alg::ManualSparseJacobian, Y, atmos)
             MatrixFields.unrolled_map(
                 name ->
                     (name, @name(c.uₕ)) =>
-                        similar(Y.c, TridiagonalRow_ACTh),
+                        similar(Y.c, TridiagonalRow_ACT12),
                 active_scalar_names,
             ) : ()
         )...,

src/utils/utilities.jl

@@ -222,43 +222,23 @@ Extracts the `g³ʰ` sub-tensor from the `gⁱʲ` tensor.
 """
 function g³ʰ(gⁱʲ)
     full_CT_axis = axes(gⁱʲ)[1]
-    CTh_axis = if full_CT_axis == Geometry.Contravariant123Axis()
-        Geometry.Contravariant12Axis()
+    N = length(full_CT_axis)
+    g³ʰ_components_view = view(Geometry.components(gⁱʲ), N:N, 1:(N - 1))
+    g³ʰ_components = if full_CT_axis == Geometry.Contravariant123Axis()
+        g³ʰ_components_view
     elseif full_CT_axis == Geometry.Contravariant13Axis()
-        Geometry.Contravariant1Axis()
+        vcat(g³ʰ_components_view, 0)
     elseif full_CT_axis == Geometry.Contravariant23Axis()
-        Geometry.Contravariant2Axis()
+        vcat(0, g³ʰ_components_view)
     else
         error("$full_CT_axis is missing either vertical or horizontal sub-axes")
     end
-    N = length(full_CT_axis)
     return Geometry.AxisTensor(
-        (Geometry.Contravariant3Axis(), CTh_axis),
-        view(Geometry.components(gⁱʲ), N:N, 1:(N - 1)),
+        (Geometry.Contravariant3Axis(), Geometry.Contravariant12Axis()),
+        g³ʰ_components,
     )
 end
 
-"""
-    CTh_vector_type(space)
-
-Extracts the (abstract) horizontal contravariant vector type from the given
-`AbstractSpace`.
-"""
-function CTh_vector_type(space)
-    full_CT_axis = axes(eltype(Fields.local_geometry_field(space).gⁱʲ))[1]
-    return if full_CT_axis == Geometry.Contravariant123Axis()
-        Geometry.Contravariant12Vector
-    elseif full_CT_axis == Geometry.Contravariant13Axis()
-        Geometry.Contravariant1Vector
-    elseif full_CT_axis == Geometry.Contravariant23Axis()
-        Geometry.Contravariant2Vector
-    elseif full_CT_axis == Geometry.Contravariant3Axis()
-        Geometry.ContravariantNullVector
-    else
-        error("$full_CT_axis is missing either vertical or horizontal sub-axes")
-    end
-end
-
 has_topography(space::Spaces.FiniteDifferenceSpace) = false
 has_topography(space) = Spaces.grid(space).hypsography != Spaces.Grids.Flat()