Try #1487:

Project.toml

@@ -19,6 +19,7 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Memoize = "c03570c3-d221-55d1-a50c-7939bbd78826"
 PkgVersion = "eebad327-c553-4316-9ea0-9fa01ccd7688"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
@@ -28,6 +29,7 @@ Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
+WeakValueDicts = "897b6980-f191-5a31-bcb0-bf3c4585e0c1"
 
 [compat]
 Adapt = "3"
@@ -43,13 +45,15 @@ GaussQuadrature = "0.5"
 GilbertCurves = "0.1"
 HDF5 = "0.16, 0.17"
 IntervalSets = "0.5, 0.6, 0.7"
+Memoize = "0.4"
 PkgVersion = "0.1, 0.2, 0.3"
 RecursiveArrayTools = "2"
 Requires = "1"
 RootSolvers = "0.3, 0.4"
 Static = "0.4, 0.5, 0.6, 0.7, 0.8"
 StaticArrays = "1"
 UnPack = "1"
+WeakValueDicts = "0.1"
 julia = "1.8"
 
 [extras]

docs/src/api.md

@@ -219,7 +219,6 @@ Spaces.dss_ghost!
 Spaces.create_dss_buffer
 Spaces.fill_send_buffer!
 Spaces.DSSBuffer
-Spaces.create_ghost_buffer
 Spaces.load_from_recv_buffer!
 Spaces.weighted_dss_start!
 Spaces.weighted_dss_internal!

examples/hybrid/hybrid3dcs_dss.jl

@@ -85,10 +85,8 @@ function hybrid3dcubedsphere_dss_profiler(
         c = center_initial_condition(ᶜlocal_geometry),
         f = face_initial_condition(ᶠlocal_geometry),
     )
-    ghost_buffer = (
-        c = Spaces.create_ghost_buffer(Y.c),
-        f = Spaces.create_ghost_buffer(Y.f),
-    )
+    ghost_buffer =
+        (c = Spaces.create_dss_buffer(Y.c), f = Spaces.create_dss_buffer(Y.f))
     dss_buffer_f = Spaces.create_dss_buffer(Y.f)
     dss_buffer_c = Spaces.create_dss_buffer(Y.c)
     nsamples = 10000

src/ClimaCore.jl

@@ -11,6 +11,8 @@ include("Geometry/Geometry.jl")
 include("Domains/Domains.jl")
 include("Meshes/Meshes.jl")
 include("Topologies/Topologies.jl")
+include("Quadratures/Quadratures.jl")
+include("Grids/Grids.jl")
 include("Spaces/Spaces.jl")
 include("Fields/Fields.jl")
 include("Operators/Operators.jl")

src/Fields/Fields.jl

@@ -5,6 +5,7 @@ import ..slab, ..slab_args, ..column, ..column_args, ..level
 import ..DataLayouts: DataLayouts, AbstractData, DataStyle
 import ..Domains
 import ..Topologies
+import ..Grids: ColumnIndex
 import ..Spaces: Spaces, AbstractSpace, AbstractPointSpace
 import ..Geometry: Geometry, Cartesian12Vector
 import ..Utilities: PlusHalf, half
@@ -38,13 +39,6 @@ Field(::Type{T}, space::S) where {T, S <: AbstractSpace} =
 
 ClimaComms.context(field::Field) = ClimaComms.context(axes(field))
 
-ClimaComms.context(space::Spaces.ExtrudedFiniteDifferenceSpace) =
-    ClimaComms.context(Spaces.horizontal_space(space))
-ClimaComms.context(space::Spaces.SpectralElementSpace2D) =
-    ClimaComms.context(space.topology)
-ClimaComms.context(space::S) where {S <: Spaces.AbstractSpace} =
-    ClimaComms.context(space.topology)
-
 ClimaComms.context(topology::Topologies.Topology2D) = topology.context
 ClimaComms.context(topology::T) where {T <: Topologies.AbstractTopology} =
     topology.context
@@ -93,11 +87,12 @@ const CenterExtrudedFiniteDifferenceField{V, S} = Field{
 } where {V <: AbstractData, S <: Spaces.CenterExtrudedFiniteDifferenceSpace}
 
 # Cubed Sphere Fields
+#=
 const CubedSphereSpectralElementField2D{V, S} = Field{
     V,
     S,
 } where {V <: AbstractData, S <: Spaces.CubedSphereSpectralElementSpace2D}
-
+=#
 
 Base.propertynames(field::Field) = propertynames(getfield(field, :values))
 @inline field_values(field::Field) = getfield(field, :values)
@@ -333,22 +328,12 @@ function interpcoord(elemrange, x::Real)
     return z, ξ
 end
 
-"""
-    Spaces.variational_solve!(field)
-
-Divide `field` by the mass matrix.
-"""
-function Spaces.variational_solve!(field::Field)
-    Spaces.variational_solve!(field_values(field), axes(field))
-    return field
-end
-
 """
     Spaces.weighted_dss!(f::Field[, ghost_buffer = Spaces.create_dss_buffer(field)])
 
 Apply weighted direct stiffness summation (DSS) to `f`. This operates in-place
 (i.e. it modifies the `f`). `ghost_buffer` contains the necessary information
-for communication in a distributed setting, see [`Spaces.create_ghost_buffer`](@ref).
+for communication in a distributed setting, see [`Spaces.create_dss_buffer`](@ref).
 
 This is a projection operation from the piecewise polynomial space
 ``\\mathcal{V}_0`` to the continuous space ``\\mathcal{V}_1 = \\mathcal{V}_0
@@ -389,13 +374,6 @@ Spaces.weighted_dss_internal!(field::Field, dss_buffer) =
 Spaces.weighted_dss_ghost!(field::Field, dss_buffer) =
     Spaces.weighted_dss_ghost!(field_values(field), axes(field), dss_buffer)
 
-"""
-    Spaces.create_ghost_buffer(field::Field)
-
-Create a buffer for communicating neighbour information of `field`.
-"""
-Spaces.create_ghost_buffer(field::Field) = Spaces.create_dss_buffer(field)
-
 """
     Spaces.create_dss_buffer(field::Field)
 
@@ -406,20 +384,6 @@ function Spaces.create_dss_buffer(field::Field)
     hspace = Spaces.horizontal_space(space)
     Spaces.create_dss_buffer(field_values(field), hspace)
 end
-# Add definitions for backward compatibility
-Spaces.weighted_dss2!(
-    field::Field,
-    dss_buffer = Spaces.create_dss_buffer(field),
-) = Spaces.weighted_dss!(field, dss_buffer)
-
-Spaces.weighted_dss_start2!(field::Field, ghost_buffer) =
-    Spaces.weighted_dss_start!(field, ghost_buffer)
-
-Spaces.weighted_dss_internal2!(field::Field, ghost_buffer) =
-    Spaces.weighted_dss_internal!(field, ghost_buffer)
-
-Spaces.weighted_dss_ghost2!(field, ghost_buffer) =
-    Spaces.weighted_dss_ghost!(field, ghost_buffer)
 
 Base.@propagate_inbounds function level(
     field::Union{

src/Fields/indices.jl

@@ -1,19 +1,3 @@
-"""
-    ColumnIndex(ij,h)
-
-An index into a column of a field. This can be used as an argument to `getindex`
-of a `Field`, to return a field on that column.
-
-# Example
-```julia
-colidx = ColumnIndex((1,1),1)
-field[colidx]
-```
-"""
-struct ColumnIndex{N}
-    ij::NTuple{N, Int}
-    h::Int
-end
 
 Base.@propagate_inbounds Base.getindex(field::Field, colidx::ColumnIndex) =
     column(field, colidx)
@@ -145,7 +129,7 @@ bycolumn(
     fn,
     space::Spaces.ExtrudedFiniteDifferenceSpace,
     device::ClimaComms.AbstractCPUDevice,
-) = bycolumn(fn, space.horizontal_space, device)
+) = bycolumn(fn, Spaces.horizontal_space(space), device)
 
 
 function bycolumn(fn, space::AbstractSpace, ::ClimaComms.CUDADevice)
@@ -164,7 +148,7 @@ Number of columns in a given space.
 ncolumns(field::Field) = ncolumns(axes(field))
 
 ncolumns(space::Spaces.ExtrudedFiniteDifferenceSpace) =
-    ncolumns(space.horizontal_space)
+    ncolumns(Spaces.horizontal_space(space))
 
 function ncolumns(space::Spaces.SpectralElementSpace1D)
     Nh = Topologies.nlocalelems(space)
@@ -219,7 +203,7 @@ function byslab(
     ::ClimaComms.CPUSingleThreaded,
     space::Spaces.AbstractSpectralElementSpace,
 )
-    Nh = Topologies.nlocalelems(space.topology)::Int
+    Nh = Topologies.nlocalelems(Spaces.topology(space))::Int
     @inbounds for h in 1:Nh
         fn(SlabIndex(nothing, h))
     end

src/Geometry/localgeometry.jl

@@ -51,3 +51,89 @@ end
 
 undertype(::Type{LocalGeometry{I, C, FT, S}}) where {I, C, FT, S} = FT
 undertype(::Type{SurfaceGeometry{FT, N}}) where {FT, N} = FT
+
+
+function blockmat(
+    a::Geometry.Axis2Tensor{
+        FT,
+        Tuple{Geometry.UAxis, Geometry.Covariant1Axis},
+        SMatrix{1, 1, FT, 1},
+    },
+    b::Geometry.Axis2Tensor{
+        FT,
+        Tuple{Geometry.WAxis, Geometry.Covariant3Axis},
+        SMatrix{1, 1, FT, 1},
+    },
+) where {FT}
+    A = Geometry.components(a)
+    B = Geometry.components(b)
+    Geometry.AxisTensor(
+        (Geometry.UWAxis(), Geometry.Covariant13Axis()),
+        SMatrix{2, 2}(A[1, 1], zero(FT), zero(FT), B[1, 1]),
+    )
+end
+
+function blockmat(
+    a::Geometry.Axis2Tensor{
+        FT,
+        Tuple{Geometry.VAxis, Geometry.Covariant2Axis},
+        SMatrix{1, 1, FT, 1},
+    },
+    b::Geometry.Axis2Tensor{
+        FT,
+        Tuple{Geometry.WAxis, Geometry.Covariant3Axis},
+        SMatrix{1, 1, FT, 1},
+    },
+) where {FT}
+    A = Geometry.components(a)
+    B = Geometry.components(b)
+    Geometry.AxisTensor(
+        (Geometry.VWAxis(), Geometry.Covariant23Axis()),
+        SMatrix{2, 2}(A[1, 1], zero(FT), zero(FT), B[1, 1]),
+    )
+end
+
+function blockmat(
+    a::Geometry.Axis2Tensor{
+        FT,
+        Tuple{Geometry.UVAxis, Geometry.Covariant12Axis},
+        SMatrix{2, 2, FT, 4},
+    },
+    b::Geometry.Axis2Tensor{
+        FT,
+        Tuple{Geometry.WAxis, Geometry.Covariant3Axis},
+        SMatrix{1, 1, FT, 1},
+    },
+) where {FT}
+    A = Geometry.components(a)
+    B = Geometry.components(b)
+    Geometry.AxisTensor(
+        (Geometry.UVWAxis(), Geometry.Covariant123Axis()),
+        SMatrix{3, 3}(
+            A[1, 1],
+            A[2, 1],
+            zero(FT),
+            A[1, 2],
+            A[2, 2],
+            zero(FT),
+            zero(FT),
+            zero(FT),
+            B[1, 1],
+        ),
+    )
+end
+
+function product_geometry(
+    horizontal_local_geometry::Geometry.LocalGeometry,
+    vertical_local_geometry::Geometry.LocalGeometry,
+)
+    coordinates = Geometry.product_coordinates(
+        horizontal_local_geometry.coordinates,
+        vertical_local_geometry.coordinates,
+    )
+    J = horizontal_local_geometry.J * vertical_local_geometry.J
+    WJ = horizontal_local_geometry.WJ * vertical_local_geometry.WJ
+    ∂x∂ξ =
+        blockmat(horizontal_local_geometry.∂x∂ξ, vertical_local_geometry.∂x∂ξ)
+    return Geometry.LocalGeometry(coordinates, J, WJ, ∂x∂ξ)
+end

src/Grids/Grids.jl

@@ -0,0 +1,73 @@
+module Grids
+
+import ClimaComms, Adapt, ForwardDiff, LinearAlgebra
+import LinearAlgebra: det, norm
+import Memoize: @memoize
+import WeakValueDicts: WeakValueDict
+import ..DataLayouts,
+    ..Domains, ..Meshes, ..Topologies, ..Geometry, ..Quadratures
+import ..Utilities: PlusHalf, half
+import ..slab, ..column, ..level
+
+using StaticArrays
+
+"""
+    Grids.AbstractGrid
+
+Grids should define the following
+
+
+- [`topology`](@ref): the topology of the grid
+- [`mesh`](@ref): the mesh of the grid
+- [`domain`](@ref): the domain of the grid
+- `ClimaComms.context`
+- `ClimaComms.device`
+
+- [`local_geometry_data`](@ref): the `DataLayout` object containing the local geometry data information
+
+"""
+abstract type AbstractGrid end
+
+
+"""
+    Grids.topology(grid::AbstractGrid)
+
+Get the topology of a grid.
+"""
+function topology end
+
+function vertical_topology end
+
+"""
+    Grids.local_geometry_data(
+        grid       :: AbstractGrid,
+        staggering :: Union{Staggering, Nothing},
+    )
+
+Get the `DataLayout` object containing the local geometry data information of
+the `grid` with staggering `staggering`.
+
+If the grid is not staggered, `staggering` should be `nothing`.
+"""
+function local_geometry_data end
+
+function local_dss_weights end
+function quadrature_style end
+function vertical_topology end
+
+
+
+ClimaComms.context(grid::AbstractGrid) = ClimaComms.context(topology(grid))
+ClimaComms.device(grid::AbstractGrid) = ClimaComms.device(topology(grid))
+
+Meshes.domain(grid::AbstractGrid) = Meshes.domain(Topologies.topology(grid))
+
+include("finitedifference.jl")
+include("spectralelement.jl")
+include("extruded.jl")
+include("column.jl")
+include("level.jl")
+
+
+
+end # module