Merge pull request #2375 from CliMA/ne/fd_grid

Add support for FiniteDifferenceSpaces

Author: nefrathenrici

NEWS.md

@@ -4,6 +4,12 @@ ClimaCore.jl Release Notes
 main
 -------
 
+
+v0.14.43
+-------
+
+- Add support for horizontal operators on column spaces [2375](https://github.com/Clima/ClimaCore.jl/pull/2375)
+
 v0.14.42
 -------
 

Project.toml

@@ -1,7 +1,7 @@
 name = "ClimaCore"
 uuid = "d414da3d-4745-48bb-8d80-42e94e092884"
 authors = ["CliMA Contributors <[email protected]>"]
-version = "0.14.42"
+version = "0.14.43"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

src/CommonGrids/CommonGrids.jl

@@ -40,7 +40,7 @@ z_domain = Domains.IntervalDomain(
 )
 z_mesh = Meshes.IntervalMesh(z_domain; nelems = z_elem)
 h_grid = Grids.SpectralElementGrid2D(h_topology, quad)
-z_topology = Topologies.IntervalTopology(context, z_mesh)
+z_topology = Topologies.IntervalTopology(ClimaComms.SingletonCommsContext(device), z_mesh)
 z_grid = Grids.FiniteDifferenceGrid(z_topology)
 grid = Grids.ExtrudedFiniteDifferenceGrid(
     h_grid,
@@ -441,6 +441,19 @@ 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/CommonGrids/Helpers.jl

@@ -27,13 +27,12 @@ function DefaultSliceXMesh(
     x_elem::Integer,
 ) where {FT}
 
-    x1boundary = (:east, :west)
-    z_boundary_names = (:bottom, :top)
+    x1boundary = periodic_x ? nothing : (:east, :west)
     h_domain = Domains.IntervalDomain(
         Geometry.XPoint{FT}(x_min),
         Geometry.XPoint{FT}(x_max);
         periodic = periodic_x,
-        boundary_names = (:east, :west),
+        boundary_names = x1boundary,
     )
     return Meshes.IntervalMesh(h_domain; nelems = x_elem)
 end
@@ -93,8 +92,9 @@ function DefaultRectangleXYMesh(
     periodic_x::Bool,
     periodic_y::Bool,
 ) where {FT <: AbstractFloat}
-    x1boundary = (:east, :west)
-    x2boundary = (:south, :north)
+    x1boundary = periodic_x ? nothing : (:east, :west)
+    x2boundary = periodic_y ? nothing : (:south, :north)
+
     domain = Domains.RectangleDomain(
         Domains.IntervalDomain(
             Geometry.XPoint{FT}(x_min),

src/Domains/Domains.jl

@@ -56,7 +56,7 @@ function IntervalDomain(
     if !periodic && isnothing(boundary_names)
         throw(
             ArgumentError(
-                "if `periodic=false` then an `boundary_names::Tuple{Symbol,Symbol}` keyword argument is required.",
+                "if `periodic=false` then a `boundary_names::Tuple{Symbol,Symbol}` keyword argument is required.",
             ),
         )
     end

src/Fields/indices.jl

@@ -245,6 +245,32 @@ function byslab(
     end
 end
 
+function byslab(
+    fn,
+    ::ClimaComms.AbstractCPUDevice,
+    space::Spaces.CenterFiniteDifferenceSpace,
+)
+    Nv = Spaces.nlevels(space)
+    @inbounds begin
+        for v in 1:Nv
+            fn(SlabIndex(v, 1))
+        end
+    end
+end
+
+function byslab(
+    fn,
+    ::ClimaComms.AbstractCPUDevice,
+    space::Spaces.FaceFiniteDifferenceSpace,
+)
+    Nv = Spaces.nlevels(space)
+    @inbounds begin
+        for v in 1:Nv
+            fn(SlabIndex(v - half, 1))
+        end
+    end
+end
+
 universal_index(colidx::Fields.ColumnIndex{2}) =
     CartesianIndex(colidx.ij[1], colidx.ij[2], 1, 1, colidx.h)
 

src/Geometry/axistensors.jl

@@ -219,6 +219,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} =
@@ -306,11 +308,12 @@ 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)]
 
-for I in [(1,), (2,), (3,), (1, 2), (1, 3), (2, 3), (1, 2, 3)]
-    strI = join(I)
+    strI = isempty(I) ? "Null" : join(I)
     N = length(I)
-    strUVW = join(map(i -> [:U, :V, :W][i], I))
+
+    strUVW = isempty(I) ? "Null" : join(map(i -> [:U, :V, :W][i], I))
     @eval begin
         const $(Symbol(:Covariant, strI, :Axis)) = CovariantAxis{$I}
         const $(Symbol(:Covariant, strI, :Vector)){T} =
@@ -467,7 +470,7 @@ end
         end
         push!(vals, val)
     end
-    return :(@inbounds AxisVector(ato, SVector($(vals...))))
+    return :(@inbounds AxisVector(ato, SVector{$(length(Ito)), $T}($(vals...))))
 end
 
 function _transform(

src/Geometry/conversions.jl

@@ -501,13 +501,16 @@ Curl is only defined for `CovariantVector`` field input types.
 |  Covariant12Vector | (1,2) | Contravariant3Vector |
 |  Covariant3Vector | (1,2) | Contravariant12Vector |
 |  Covariant123Vector | (1,2) | Contravariant123Vector |
-|  Covariant1Vector | (1,) | Contravariant1Vector |
+|  Covariant1Vector | (1,) | ContravariantNullVector |
 |  Covariant2Vector | (1,) | Contravariant3Vector |
 |  Covariant3Vector | (1,) | Contravariant2Vector |
+|  Covariant13Vector | (1,) | Contravariant2Vector |
+|  Covariant123Vector | (1,) | Contravariant23Vector |
 |  Covariant12Vector | (3,) | Contravariant12Vector |
 |  Covariant1Vector | (3,) | Contravariant2Vector |
 |  Covariant2Vector | (3,) | Contravariant1Vector |
 |  Covariant3Vector | (3,) | Contravariant3Vector |
+|  Any CovariantVector | () | ContravariantNullVector |
 """
 @inline curl_result_type(::Val{(1, 2)}, ::Type{Covariant3Vector{FT}}) where {FT} =
     Contravariant12Vector{FT}
@@ -516,9 +519,8 @@ Curl is only defined for `CovariantVector`` field input types.
 @inline curl_result_type(::Val{(1, 2)}, ::Type{Covariant123Vector{FT}}) where {FT} =
     Contravariant123Vector{FT}
 
-
 @inline curl_result_type(::Val{(1,)}, ::Type{Covariant1Vector{FT}}) where {FT} =
-    Contravariant1Vector{FT} # not strictly correct: should be a zero Vector
+    ContravariantNullVector{FT}
 @inline curl_result_type(::Val{(1,)}, ::Type{Covariant2Vector{FT}}) where {FT} =
     Contravariant3Vector{FT}
 @inline curl_result_type(::Val{(1,)}, ::Type{Covariant3Vector{FT}}) where {FT} =
@@ -537,6 +539,9 @@ Curl is only defined for `CovariantVector`` field input types.
 @inline curl_result_type(::Val{(3,)}, ::Type{Covariant3Vector{FT}}) where {FT} =
     Contravariant3Vector{FT}
 
+@inline curl_result_type(_, ::Type{<:CovariantVector{FT}}) where {FT} =
+    ContravariantNullVector{FT}
+
 _norm_sqr(x, local_geometry::LocalGeometry) = sum(x -> _norm_sqr(x, local_geometry), x)
 _norm_sqr(x::Number, ::LocalGeometry) = LinearAlgebra.norm_sqr(x)
 _norm_sqr(x::AbstractArray, ::LocalGeometry) = LinearAlgebra.norm_sqr(x)