address comments

Author: nefrathenrici

src/CommonGrids/CommonGrids.jl

@@ -441,19 +441,6 @@ function Box3DGrid(
             periodic_x,
             periodic_y,
         ),
-        Topologies.spacefillingcurve(
-            DefaultRectangleXYMesh(
-                FT;
-                x_min,
-                x_max,
-                y_min,
-                y_max,
-                x_elem,
-                y_elem,
-                periodic_x,
-                periodic_y,
-            ),
-        ),
     ),
     z_mesh::Meshes.IntervalMesh = DefaultZMesh(
         FT;

src/Geometry/axistensors.jl

@@ -234,6 +234,8 @@ const AxisVector{T, A1, S} = AxisTensor{T, 1, Tuple{A1}, S}
 AxisVector(ax::A1, v::SVector{N, T}) where {A1 <: AbstractAxis, N, T} =
     AxisVector{T, A1, SVector{N, T}}((ax,), v)
 
+(AxisVector{T, A, SVector{0, T}} where {T})() where {A} =
+    AxisVector(A.instance, SVector{0, T}())
 (AxisVector{T, A, SVector{1, T}} where {T})(arg1::Real) where {A} =
     AxisVector(A.instance, SVector(arg1))
 (AxisVector{T, A, SVector{2, T}} where {T})(arg1::Real, arg2::Real) where {A} =
@@ -244,8 +246,6 @@ AxisVector(ax::A1, v::SVector{N, T}) where {A1 <: AbstractAxis, N, T} =
     arg3::Real,
 ) where {A} = AxisVector(A.instance, SVector(arg1, arg2, arg3))
 
-(AxisVector{T, A, SVector{0, T}})() where {A, T} =
-    AxisVector(A.instance, SVector{0, T}())
 
 const CovariantVector{T, I, S} = AxisVector{T, CovariantAxis{I}, S}
 const ContravariantVector{T, I, S} = AxisVector{T, ContravariantAxis{I}, S}
@@ -331,6 +331,10 @@ const Cartesian2Tensor{T, A, S} =
 const Local2Tensor{T, A, S} =
     Axis2Tensor{T, A, S} where {T, A <: Tuple{LocalAxis, AbstractAxis}, S}
 
+const CovariantTensor = Union{CovariantVector, Covariant2Tensor}
+const ContravariantTensor = Union{ContravariantVector, Contravariant2Tensor}
+const CartesianTensor = Union{CartesianVector, Cartesian2Tensor}
+const LocalTensor = Union{LocalVector, Local2Tensor}
 for I in [(), (1,), (2,), (3,), (1, 2), (1, 3), (2, 3), (1, 2, 3)]
 
     strI = isempty(I) ? "Null" : join(I)
@@ -355,12 +359,6 @@ for I in [(), (1,), (2,), (3,), (1, 2), (1, 3), (2, 3), (1, 2, 3)]
     end
 end
 
-const CovariantTensor = Union{CovariantVector, Covariant2Tensor}
-const ContravariantTensor =
-    Union{ContravariantVector, Contravariant2Tensor}
-const CartesianTensor = Union{CartesianVector, Cartesian2Tensor}
-const LocalTensor = Union{LocalVector, Local2Tensor}
-
 # LinearAlgebra
 
 check_axes(::A, ::A) where {A} = nothing

src/Geometry/conversions.jl

@@ -550,10 +550,6 @@ The return type when taking the divergence of a field of `V`.
 Required for statically infering the result type of the divergence operation for `AxisVector` subtypes.
 """
 @inline divergence_result_type(::Type{V}) where {V <: AxisVector} = eltype(V)
-@inline divergence_result_type(::Type{V}) where {V <: CovariantNullVector} =
-    eltype(V)
-@inline divergence_result_type(::Type{V}) where {V <: ContravariantNullVector} =
-    eltype(V)
 
 # this isn't quite right: it only is true when the Christoffel symbols are zero
 @inline divergence_result_type(
@@ -575,9 +571,6 @@ Required for statically infering the result type of the gradient operation for `
     N = length(I)
     AxisVector{V, CovariantAxis{I}, SVector{N, V}}
 end
-@inline function gradient_result_type(::Val{()}, ::Type{V}) where {V <: Number}
-    CovariantNullVector{V}
-end
 @inline function gradient_result_type(
     ::Val{I},
     ::Type{V},
@@ -625,7 +618,8 @@ Curl is only defined for `CovariantVector`` field input types.
 ) where {FT} = ContravariantNullVector{FT}
 @inline curl_result_type(::Val{()}, ::Type{<:CovariantVector{FT}}) where {FT} =
     ContravariantNullVector{FT}
-
+@inline curl_result_type(_, ::Type{CovariantNullVector{FT}}) where {FT} =
+    ContravariantNullVector{FT}
 
 @inline curl_result_type(::Val{(1,)}, ::Type{Covariant1Vector{FT}}) where {FT} =
     Contravariant1Vector{FT} # not strictly correct: should be a zero Vector

src/Operators/spectralelement.jl

@@ -188,27 +188,13 @@ Copy the slab indexed by `slabidx` from `bc` to `out`.
 """
 Base.@propagate_inbounds function copyto_slab!(out, bc, slabidx)
     space = axes(out)
-    QS = Spaces.quadrature_style(space)
-    Nq = Quadratures.degrees_of_freedom(QS)
     rbc = resolve_operator(bc, slabidx)
     @inbounds for ij in node_indices(axes(out))
         set_node!(space, out, ij, slabidx, get_node(space, rbc, ij, slabidx))
     end
     return nothing
 end
 
-Base.@propagate_inbounds function copyto_slab!(
-    out::Field{T, <:Spaces.FiniteDifferenceSpace},
-    bc,
-    slabidx,
-) where {T}
-    space = axes(out)
-    rbc = resolve_operator(bc, slabidx)
-    @inbounds for ij in node_indices(space)
-        set_node!(space, out, ij, slabidx, get_node(space, rbc, ij, slabidx))
-    end
-    return nothing
-end
 
 """
     resolve_operator(bc, slabidx)
@@ -355,7 +341,7 @@ Base.@propagate_inbounds function get_node(
     slabidx,
 )
     space = reconstruct_placeholder_space(axes(field), parent_space)
-    i, = space isa Spaces.FiniteDifferenceSpace ? (1,) : Tuple(ij)
+    i, = Tuple(ij)
     if space isa Spaces.FaceExtrudedFiniteDifferenceSpace ||
        space isa Spaces.FaceFiniteDifferenceSpace
         _v = slabidx.v + half
@@ -1635,6 +1621,7 @@ matrix_interpolate(field::Field, Nu::Integer) =
     matrix_interpolate(field, Quadratures.Uniform{Nu}())
 
 import .DataLayouts: slab_index
+import .Spaces: slab_type
 
 """
     rmatmul1(W, S, i, j)
@@ -1670,6 +1657,7 @@ function rmatmul2(W, S, i, j)
 end
 
 
+
 function apply_operator(
     op::Union{
         Divergence{()},
@@ -1686,7 +1674,18 @@ function apply_operator(
     FT = Spaces.undertype(space)
     RT = operator_return_eltype(op, eltype(arg))
     zero_value = zero(RT)
-    out = DataLayouts.IF{RT, 1}(MArray, FT)
-    out[slab_index(1)] = zero_value
+    slab_type = DataLayouts.slab_type(space)
+    QS = Spaces.quadrature_style(space)
+    Nq = Quadratures.degrees_of_freedom(QS)
+
+    if slab_type === DataLayouts.IJF
+        out = DataLayouts.IJF{RT, Nq}(MArray, FT)
+        fill!(parent(out), zero(FT))
+        out[slab_index(1, 1)] = zero_value
+    else
+        out = DataLayouts.IF{RT, 1}(MArray, FT)
+        out[slab_index(1)] = zero_value
+    end
+
     return Field(SArray(out), space)
 end

src/Spaces/Spaces.jl

@@ -212,4 +212,18 @@ set_mask!(fn, space::ExtrudedFiniteDifferenceSpace) =
 set_mask!(space::AbstractSpace, data::DataLayouts.AbstractData) =
     set_mask!(grid(space), data)
 
+"""
+    slab_type(space)
+
+Determines the appropriate slab data layout type for a given space.
+
+For spaces with 2 horizontal dimensions, returns IJF.
+For 1D spaces, returns IF.
+"""
+slab_type(space::SpectralElementSpace2D) = DataLayouts.IJF
+slab_type(space::SpectralElementSpace1D) = DataLayouts.IF
+slab_type(space::FiniteDifferenceSpace) = DataLayouts.IF
+slab_type(space::ExtrudedFiniteDifferenceSpace) =
+    slab_type(horizontal_space(space))
+
 end # module

src/Topologies/Topologies.jl

@@ -63,7 +63,6 @@ end
 
 function print_context(io::IO, context::ClimaComms.MPICommsContext)
     print(io, "MPICommsContext with ", ClimaComms.nprocs(context), " processes")
-
     print(io, " using ")
     print_device(io, context.device)
 end

src/Topologies/topology2d.jl

@@ -237,7 +237,7 @@ end
 function Topology2D(
     context::ClimaComms.AbstractCommsContext,
     mesh::Meshes.AbstractMesh{2},
-    elemorder = Meshes.elements(mesh),
+    elemorder = spacefillingcurve(mesh),
     elempid = nothing, # array of same size as elemorder, containing owning pid for each element (it should be sorted)
     orderindex = Meshes.linearindices(elemorder), # inverse mapping of elemorder (e.g. map CartesianIndex => Int)
 )

src/to_device.jl

@@ -18,10 +18,7 @@ If the input is already defined on the target device, returns a copy.
 
 This means that `out === x` will not in general be satisfied.
 """
-function to_device(
-    device::ClimaComms.AbstractDevice,
-    x,
-)
+function to_device(device::ClimaComms.AbstractDevice, x)
     return Adapt.adapt(ClimaComms.array_type(device), x)
 end