diff --git a/src/Domains/Domains.jl b/src/Domains/Domains.jl index 5336f30a2c..4838c7ca95 100644 --- a/src/Domains/Domains.jl +++ b/src/Domains/Domains.jl @@ -35,6 +35,10 @@ struct IntervalDomain{CT, B} <: AbstractDomain where { boundary_names::B end +const XIntervalDomain{FT, B} = IntervalDomain{Geometry.XPoint{FT}, B} +const YIntervalDomain{FT, B} = IntervalDomain{Geometry.YPoint{FT}, B} +const ZIntervalDomain{FT, B} = IntervalDomain{Geometry.ZPoint{FT}, B} + isperiodic(domain::IntervalDomain) = isnothing(domain.boundary_names) boundary_names(domain::IntervalDomain) = isperiodic(domain) ? () : unique(domain.boundary_names) @@ -79,6 +83,49 @@ The domain minimum along the z-direction. """ z_min(domain::IntervalDomain) = domain.coord_min.z +function XIntervalDomain(; + x_min::Real, + x_max::Real, + x_periodic::Bool = false, + x_boundary_names = (:west, :east), +) + x_min, x_max = promote(x_min, x_max) + IntervalDomain( + Geometry.XPoint(x_min), + Geometry.XPoint(x_max); + periodic = x_periodic, + boundary_names = x_boundary_names, + ) +end +function YIntervalDomain(; + y_min::Real, + y_max::Real, + y_periodic::Bool = false, + y_boundary_names = (:south, :north), +) + y_min, y_max = promote(y_min, y_max) + IntervalDomain( + Geometry.YPoint(y_min), + Geometry.YPoint(y_max); + periodic = y_periodic, + boundary_names = y_boundary_names, + ) +end +function ZIntervalDomain(; + z_min::Real, + z_max::Real, + z_periodic::Bool = false, + z_boundary_names = (:bottom, :top), +) + z_min, z_max = promote(z_min, z_max) + IntervalDomain( + Geometry.ZPoint(z_min), + Geometry.ZPoint(z_max); + periodic = z_periodic, + boundary_names = z_boundary_names, + ) +end + coordinate_type(::IntervalDomain{CT}) where {CT} = CT Base.eltype(domain::IntervalDomain) = coordinate_type(domain) @@ -108,6 +155,9 @@ struct RectangleDomain{I1 <: IntervalDomain, I2 <: IntervalDomain} <: interval1::I1 interval2::I2 end + +const XYRectangleDomain = RectangleDomain{<:XIntervalDomain, <:YIntervalDomain} + Base.:*(interval1::IntervalDomain, interval2::IntervalDomain) = RectangleDomain(interval1, interval2) @@ -144,6 +194,20 @@ function RectangleDomain( return interval1 * interval2 end +function XYRectangleDomain(; + x_min::Real, + x_max::Real, + x_periodic::Bool = false, + x_boundary_names = (:west, :east), + y_min::Real, + y_max::Real, + y_periodic::Bool = false, + y_boundary_names = (:south, :north), +) + x_domain = XIntervalDomain(; x_min, x_max, x_periodic, x_boundary_names) + y_domain = YIntervalDomain(; y_min, y_max, y_periodic, y_boundary_names) + RectangleDomain(x_domain, y_domain) +end function Base.show(io::IO, domain::RectangleDomain) print(io, nameof(typeof(domain)), ": ") diff --git a/src/Grids/column.jl b/src/Grids/column.jl index d31825f6cb..5c30d9a9c0 100644 --- a/src/Grids/column.jl +++ b/src/Grids/column.jl @@ -21,14 +21,14 @@ end """ - ColumnGrid( + ColumnViewGrid( full_grid :: ExtrudedFiniteDifferenceGrid, colidx :: ColumnIndex, ) A view into a column of a `ExtrudedFiniteDifferenceGrid`. This can be used as an """ -struct ColumnGrid{ +struct ColumnViewGrid{ G <: AbstractExtrudedFiniteDifferenceGrid, C <: ColumnIndex, } <: AbstractFiniteDifferenceGrid @@ -40,14 +40,15 @@ local_geometry_type(::Type{ColumnGrid{G, C}}) where {G, C} = local_geometry_type(G) column(grid::AbstractExtrudedFiniteDifferenceGrid, colidx::ColumnIndex) = - ColumnGrid(grid, colidx) + ColumnViewGrid(grid, colidx) -topology(colgrid::ColumnGrid) = vertical_topology(colgrid.full_grid) -vertical_topology(colgrid::ColumnGrid) = vertical_topology(colgrid.full_grid) +topology(colgrid::ColumnViewGrid) = vertical_topology(colgrid.full_grid) +vertical_topology(colgrid::ColumnViewGrid) = + vertical_topology(colgrid.full_grid) -local_geometry_data(colgrid::ColumnGrid, staggering::Staggering) = column( +local_geometry_data(colgrid::ColumnViewGrid, staggering::Staggering) = column( local_geometry_data(colgrid.full_grid, staggering::Staggering), colgrid.colidx.ij..., colgrid.colidx.h, ) -global_geometry(colgrid::ColumnGrid) = global_geometry(colgrid.full_grid) +global_geometry(colgrid::ColumnViewGrid) = global_geometry(colgrid.full_grid) diff --git a/src/Grids/extruded.jl b/src/Grids/extruded.jl index 15b68c9423..f8c90d41b4 100644 --- a/src/Grids/extruded.jl +++ b/src/Grids/extruded.jl @@ -175,7 +175,134 @@ const ExtrudedSpectralElementGrid2D = ExtrudedFiniteDifferenceGrid{<:SpectralElementGrid1D} const ExtrudedSpectralElementGrid3D = ExtrudedFiniteDifferenceGrid{<:SpectralElementGrid2D} -const ExtrudedRectilinearSpectralElementGrid3D = - ExtrudedFiniteDifferenceGrid{<:RectilinearSpectralElementGrid2D} -const ExtrudedCubedSphereSpectralElementGrid3D = + + + +const SliceGrid = ExtrudedFiniteDifferenceGrid{<:LineSpectralElementGrid} +function SliceGrid(; + x_min::Real, + x_max::Real, + x_elem, + x_periodic::Bool = false, + x_boundary_names = (:west, :east), + poly_degree = 3, + z_min::Real, + z_max::Real, + z_periodic::Bool = false, + z_boundary_names = (:bottom, :top), + z_elem::Integer, + z_stretch = Meshes.Uniform(), + context = ClimaComms.context(), +) + h_grid = LineSpectralElementGrid(; + x_min, + x_max, + x_elem, + x_periodic, + x_boundary_names, + poly_degree, + context, + ) + v_grid = ColumnGrid(; + z_min, + z_max, + z_periodic, + z_boundary_names, + z_elem, + z_stretch, + context = ClimaComms.SingletonCommsContext(ClimaComms.device(context)), + ) + ExtrudedFiniteDifferenceGrid(h_grid, v_grid) +end + + + +const BoxGrid = ExtrudedFiniteDifferenceGrid{<:RectilinearSpectralElementGrid} + +function BoxGrid(; + x_min::Real, + x_max::Real, + x_elem, + x_periodic::Bool = false, + x_boundary_names = (:west, :east), + y_min::Real, + y_max::Real, + y_elem, + y_periodic::Bool = false, + y_boundary_names = (:south, :north), + poly_degree = 3, + z_min::Real, + z_max::Real, + z_periodic::Bool = false, + z_boundary_names = (:bottom, :top), + z_elem::Integer, + z_stretch = Meshes.Uniform(), + context = ClimaComms.context(), +) + h_grid = RectilinearSpectralElementGrid2D(; + x_min, + x_max, + x_elem, + x_periodic, + x_boundary_names, + y_min, + y_max, + y_elem, + y_periodic, + y_boundary_names, + poly_degree, + context, + ) + v_grid = ColumnGrid(; + z_min, + z_max, + z_periodic, + z_boundary_names, + z_elem, + z_stretch, + context = ClimaComms.SingletonCommsContext(ClimaComms.device(context)), + ) + ExtrudedFiniteDifferenceGrid(h_grid, v_grid) +end + + +const ExtrudedCubedSphereGrid = ExtrudedFiniteDifferenceGrid{<:CubedSphereSpectralElementGrid2D} + +function ExtrudedCubedSphereGrid(; + radius::Real, + panel_elem::Integer, + cubed_sphere_type = Meshes.EquiangularCubedSphere, + poly_degree = 3, + bubble = true, + z_min::Real, + z_max::Real, + z_periodic::Bool = false, + z_boundary_names = (:bottom, :top), + z_elem::Integer, + z_stretch = Meshes.Uniform(), + context = ClimaComms.context(), +) + h_grid = CubedSphereGrid(; + radius, + panel_elem, + cubed_sphere_type, + poly_degree, + bubble, + context, + ) + v_grid = ColumnGrid(; + z_min, + z_max, + z_periodic, + z_boundary_names, + z_elem, + z_stretch, + context = ClimaComms.SingletonCommsContext(ClimaComms.device(context)), + ) + ExtrudedFiniteDifferenceGrid(h_grid, v_grid) +end + +## to be deprecated +const ExtrudedCubedSphereSpectralElementGrid3D = ExtrudedCubedSphereGrid +const ExtrudedRectilinearSpectralElementGrid3D = BoxGrid diff --git a/src/Grids/finitedifference.jl b/src/Grids/finitedifference.jl index 8f594ec417..a39189412b 100644 --- a/src/Grids/finitedifference.jl +++ b/src/Grids/finitedifference.jl @@ -206,3 +206,34 @@ local_geometry_data(grid::DeviceFiniteDifferenceGrid, ::CellCenter) = local_geometry_data(grid::DeviceFiniteDifferenceGrid, ::CellFace) = grid.face_local_geometry global_geometry(grid::DeviceFiniteDifferenceGrid) = grid.global_geometry + + +## aliases +const ColumnGrid = FiniteDifferenceGrid{<:Topologies.ColumnTopology1D} + +""" + Grids.ColumnGrid(; + z_min, z_max, + context = ClimaComms.SingletonCommsContext() + ) +""" +function ColumnGrid(; + z_min::Real, + z_max::Real, + z_periodic::Bool = false, + z_boundary_names = (:bottom, :top), + z_elem::Integer, + z_stretch = Meshes.Uniform(), + context = ClimaComms.SingletonCommsContext(), +) + mesh = Meshes.ZIntervalMesh(; + z_min, + z_max, + z_periodic, + z_boundary_names, + z_elem, + z_stretch, + ) + topology = Topologies.IntervalTopology(mesh) + FiniteDifferenceGrid(topology) +end diff --git a/src/Grids/spectralelement.jl b/src/Grids/spectralelement.jl index e49fe84b9b..fbd454ccc8 100644 --- a/src/Grids/spectralelement.jl +++ b/src/Grids/spectralelement.jl @@ -148,7 +148,7 @@ local_geometry_type( SpectralElementSpace2D(topology, quadrature_style; enable_bubble, horizontal_layout_type = DataLayouts.IJFH) Construct a `SpectralElementSpace2D` instance given a `topology` and `quadrature`. The -flag `enable_bubble` enables the `bubble correction` for more accurate element areas. +flag `bubble` enables the "bubble correction" for more accurate element areas. # Input arguments: - topology: Topology2D @@ -156,7 +156,7 @@ flag `enable_bubble` enables the `bubble correction` for more accurate element a - enable_bubble: Bool - horizontal_layout_type: Type{<:AbstractData} -The idea behind the so-called `bubble_correction` is that the numerical area +The idea behind the so-called "bubble correction" is that the numerical area of the domain (e.g., the sphere) is given by the sum of nodal integration weights times their corresponding Jacobians. However, this discrete sum is not exactly equal to the exact geometric area (4pi*radius^2 for the sphere). To make these equal, @@ -322,7 +322,7 @@ function _SpectralElementGrid2D( end # If enabled, apply bubble correction - if enable_bubble + if bubble if abs(elem_area - high_order_elem_area) ≤ eps(FT) for i in 1:Nq, j in 1:Nq ξ = SVector(quad_points[i], quad_points[j]) @@ -605,7 +605,95 @@ Adapt.adapt_structure(to, grid::SpectralElementGrid2D) = ) ## aliases -const RectilinearSpectralElementGrid2D = + +const LineSpectralElementGrid = + SpectralElementGrid1D{<:Topologies.LineTopology1D} + +function LineSpectralElementGrid(; + x_min::Real, + x_max::Real, + x_periodic::Bool = false, + x_boundary_names = (:west, :east), + x_elem::Integer, + x_stretch = Uniform(), + poly_degree = 3, + context = ClimaComms.SingletonCommsContext(), +) + mesh = XIntervalMesh(; + x_min, + x_max, + x_periodic, + x_boundary_names, + x_elem, + x_stretch, + ) + topology = IntervalTopology(context, mesh) + quadrature_style = Quadratures.GLL{poly_degree + 1}() + SpectralElementGrid1D(topology, quadrature_style) +end +const RectilinearSpectralElementGrid = SpectralElementGrid2D{<:Topologies.RectilinearTopology2D} -const CubedSphereSpectralElementGrid2D = + +""" + Grids.RectilinearSpectralElementGrid(; + x_min, x_max, x_elem, x_periodic=false, x_boundary_names = (:west, :east), + y_min, y_max, y_elem, y_periodic=false, y_boundary_names = (:south, :north), + context = ClimaComms.context(), + poly_degree = 3, + ) + +Constructor for a `SpectralElementGrid2D` on a `RectangleDomain`, with `XYPoint` coordinates. +""" +function RectilinearSpectralElementGrid(; + x_min::Real, + x_max::Real, + x_elem, + x_periodic::Bool = false, + x_boundary_names = (:west, :east), + y_min::Real, + y_max::Real, + y_elem, + y_periodic::Bool = false, + y_boundary_names = (:south, :north), + context = ClimaComms.context(), + poly_degree = 3, +) + mesh = Meshes.RectilinearMesh(; + x_min, + x_max, + x_periodic, + x_boundary_names, + x_elem, + y_min, + y_max, + y_periodic, + y_boundary_names, + y_elem, + ) + topology = Topologies.Topology2D(context, mesh) + quadrature_style = Quadratures.GLL{poly_degree + 1}() + SpectralElementGrid2D(topology, quadrature_style) +end + +const CubedSphereGrid = SpectralElementGrid2D{<:Topologies.CubedSphereTopology2D} + +function CubedSphereGrid(; + radius::Real, + panel_elem::Integer, + cubed_sphere_type = Meshes.EquiangularCubedSphere, + context = ClimaComms.context(), + poly_degree = 3, + bubble = true, +) + domain = Domains.SphereDomain(radius) + mesh = cubed_sphere_type(domain, panel_elem) + topology = Topologies.Topology2D(context, mesh) + quadrature_style = Quadratures.GLL{poly_degree + 1}() + SpectralElementGrid2D(topology, quadrature_style; bubble) +end + + +## to be deprecated +const RectilinearSpectralElementGrid2D = RectilinearSpectralElementGrid +const CubedSphereSpectralElementGrid2D = CubedSphereGrid diff --git a/src/Meshes/interval.jl b/src/Meshes/interval.jl index 7a2ea8d755..4e02807eb6 100644 --- a/src/Meshes/interval.jl +++ b/src/Meshes/interval.jl @@ -454,3 +454,33 @@ function truncate_mesh( ) return IntervalMesh(new_domain, new_stretch; nelems = new_nelems) end + +## aliases + +const XIntervalMesh = IntervalMesh{<:Domains.XIntervalDomain} +function XIntervalMesh(; + x_min::Real, + x_max::Real, + x_periodic::Bool = false, + x_boundary_names = (:west, :east), + x_elem::Integer, + x_stretch = Uniform(), +) + domain = + Domains.XIntervalDomain(; x_min, x_max, x_periodic, x_boundary_names) + mesh = IntervalMesh(domain, x_stretch; nelem = x_elem) +end + +const ZIntervalMesh = IntervalMesh{<:Domains.ZIntervalDomain} +function ZIntervalMesh(; + z_min::Real, + z_max::Real, + z_periodic::Bool = false, + z_boundary_names = (:bottom, :top), + z_elem::Integer, + z_stretch = Uniform(), +) + domain = + Domains.ZIntervalDomain(; z_min, z_max, z_periodic, z_boundary_names) + mesh = IntervalMesh(domain, z_stretch; nelem = z_elem) +end diff --git a/src/Meshes/rectangle.jl b/src/Meshes/rectangle.jl index a9b0393682..9f050c156e 100644 --- a/src/Meshes/rectangle.jl +++ b/src/Meshes/rectangle.jl @@ -32,6 +32,27 @@ function Base.hash(mesh::RectilinearMesh, h::UInt) end +function RectilinearMesh(; + x_min::Real, + x_max::Real, + x_elem::Integer, + x_periodic::Bool = false, + x_boundary_names = (:west, :east), + y_min::Real, + y_max::Real, + y_elem::Integer, + y_periodic::Bool = false, + y_boundary_names = (:south, :north), +) + x_domain = XIntervalDomain(; x_min, x_max, x_periodic, x_boundary_names) + y_domain = YIntervalDomain(; y_min, y_max, y_periodic, y_boundary_names) + + RectilinearMesh( + IntervalMesh(x_domain, x_elem), + IntervalMesh(y_domain, y_elem), + ) +end + function Base.summary(io::IO, mesh::RectilinearMesh) n1 = nelements(mesh.intervalmesh1) n2 = nelements(mesh.intervalmesh2) diff --git a/src/Topologies/interval.jl b/src/Topologies/interval.jl index f205a79e74..86f07b112d 100644 --- a/src/Topologies/interval.jl +++ b/src/Topologies/interval.jl @@ -165,3 +165,11 @@ function local_neighboring_elements(topology::AbstractIntervalTopology, elem) end end ghost_neighboring_elements(topology::AbstractIntervalTopology, elem) = () + + +## aliases +const LineTopology1D = + IntervalTopology{<:ClimaComms.AbstractCommsContext, <:Meshes.XIntervalMesh} + +const ColumnTopology1D = + IntervalTopology{<:ClimaComms.AbstractCommsContext, <:Meshes.ZIntervalMesh}