Skip to content

Add space with multiple columns#2525

Open
ph-kev wants to merge 1 commit into
mainfrom
kp/multiple-columns
Open

Add space with multiple columns#2525
ph-kev wants to merge 1 commit into
mainfrom
kp/multiple-columns

Conversation

@ph-kev

@ph-kev ph-kev commented Jun 10, 2026

Copy link
Copy Markdown
Member

This PR adds a space that can supports multiple columns. This is needed to complete CliMA/ClimaLand.jl#1723.

This PR have been tested with the FluxnetSimulations in ClimaLand and it produces identical results with the current single column used for the FluxnetSimulation.

The remapping feature will not be done in this PR. For the purpose of calibration, this is not needed because the data is extracted from the ClimaCore fields rather than from the NetCDF files produced from the diagnostics.

Note that most of the code in this PR have been written with Claude.

@ph-kev ph-kev force-pushed the kp/multiple-columns branch 6 times, most recently from 7567810 to e78f8a4 Compare June 10, 2026 23:20

@nefrathenrici nefrathenrici left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks good, I only saw a few immediate correctness issues.

Comment thread test/TestUtilities/TestUtilities.jl Outdated
Comment thread src/Grids/pointcloud.jl Outdated
Comment thread src/Spaces/extruded.jl Outdated
@ph-kev ph-kev marked this pull request as ready for review June 12, 2026 17:53
@ph-kev ph-kev force-pushed the kp/multiple-columns branch from e78f8a4 to 02c2eb2 Compare June 12, 2026 20:28
Comment thread ext/cuda/adapt.jl

# PointCloudGrid has no quadrature_style field; pass nothing so that FD
# operators on PointColumnEnsembleSpace work on CUDA.
Adapt.adapt_structure(

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it be better to combine this with the above method

Adapt.adapt_structure(
    to::CUDA.KernelAdaptor,
    grid::Grids.ExtrudedFiniteDifferenceGrid,
) = Grids.DeviceExtrudedFiniteDifferenceGrid(
    Adapt.adapt(to, Grids.vertical_topology(grid)),
    hasproperty(grid, :horizontal_grid.quadrature_style), Adapt.adapt(to, grid.horizontal_grid.quadrature_style),
    Adapt.adapt(to, grid.global_geometry),
    Adapt.adapt(to, grid.center_local_geometry),
    Adapt.adapt(to, grid.face_local_geometry),
)

Comment thread src/Grids/extruded.jl

function ExtrudedFiniteDifferenceGrid(
horizontal_grid::Union{SpectralElementGrid1D, SpectralElementGrid2D},
horizontal_grid::Union{SpectralElementGrid1D, SpectralElementGrid2D, PointCloudGrid},

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This union appears a lot. Could it just be a const or abstract type?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah aren't these all subtypes of AbstractSpectralElementGrid? Could just use that instead.

Comment thread src/Grids/pointcloud.jl
diagonal matrix `diag(R·π/180, R·cosd(lat)·π/180)` reflecting the sphere metric at
each point, and `J = R²·cosd(lat)·(π/180)²`.
"""
struct PointCloudGrid{

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it might be better to add a field for topology, and make it always nothing.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This could simplify all the topology(grid) methods, but subsequent operations on the result of topology would still be invalid. We could make some sort of NoTopology struct, or add methods to topology-related functions so they can handle Nothing, but that feels roughly equivalent to what this PR is doing anyway. So I don't think adding a field for the topology would be very useful.

@dennisYatunin dennisYatunin left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks great! It has a fair bit of code duplication, both in terms of method definitions and const type declarations, but that's a higher-level design issue that can be fixed separately. I left a few nitpicks about readability and error handling, but overall it looks ready to merge in.

Comment on lines +776 to +777
@assert context isa ClimaComms.SingletonCommsContext "PointColumnEnsembleGrid only supports SingletonCommsContext."
@assert ClimaComms.device(context) == device "The given device and context device do not match."

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
@assert context isa ClimaComms.SingletonCommsContext "PointColumnEnsembleGrid only supports SingletonCommsContext."
@assert ClimaComms.device(context) == device "The given device and context device do not match."
context = ClimaComms.SingletonCommsContext(device)

No need to provide a keyword argument for context if it can only have one possible value.

Comment thread src/Grids/extruded.jl

function ExtrudedFiniteDifferenceGrid(
horizontal_grid::Union{SpectralElementGrid1D, SpectralElementGrid2D},
horizontal_grid::Union{SpectralElementGrid1D, SpectralElementGrid2D, PointCloudGrid},

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah aren't these all subtypes of AbstractSpectralElementGrid? Could just use that instead.

Comment thread src/Grids/extruded.jl
Comment on lines +129 to +134
# For PointCloudGrid there is no horizontal topology; route context/device
# through the horizontal grid directly (which stores its own context).
ClimaComms.context(grid::ExtrudedFiniteDifferenceGrid{<:PointCloudGrid}) =
ClimaComms.context(grid.horizontal_grid)
ClimaComms.device(grid::ExtrudedFiniteDifferenceGrid{<:PointCloudGrid}) =
ClimaComms.device(grid.horizontal_grid)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# For PointCloudGrid there is no horizontal topology; route context/device
# through the horizontal grid directly (which stores its own context).
ClimaComms.context(grid::ExtrudedFiniteDifferenceGrid{<:PointCloudGrid}) =
ClimaComms.context(grid.horizontal_grid)
ClimaComms.device(grid::ExtrudedFiniteDifferenceGrid{<:PointCloudGrid}) =
ClimaComms.device(grid.horizontal_grid)
ClimaComms.context(grid::ExtrudedFiniteDifferenceGrid) =
ClimaComms.context(grid.horizontal_grid)
ClimaComms.device(grid::ExtrudedFiniteDifferenceGrid) =
ClimaComms.device(grid.horizontal_grid)

This could also be the default for every ExtrudedFiniteDifferenceGrid. Special cases for such specific subtypes are generally harder to maintain.

Comment thread src/Grids/pointcloud.jl
Comment on lines +15 to +18
The `local_geometry` is stored as an `IFH{LG, 1, N}` data layout — one node per
"element" (`Ni = 1`), N "elements" (`Nh = N`). The horizontal Jacobian `∂x∂ξ` is the
diagonal matrix `diag(R·π/180, R·cosd(lat)·π/180)` reflecting the sphere metric at
each point, and `J = R²·cosd(lat)·(π/180)²`.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
The `local_geometry` is stored as an `IFH{LG, 1, N}` data layout — one node per
"element" (`Ni = 1`), N "elements" (`Nh = N`). The horizontal Jacobian `∂x∂ξ` is the
diagonal matrix `diag(R·π/180, R·cosd(lat)·π/180)` reflecting the sphere metric at
each point, and `J = R²·cosd(lat)·(π/180)²`.
The `local_geometry` is stored as an `IFH{LG, 1, N}` data layout, with each of
the `N` locations represented by an element with one nodal point. Based on the
[metric tensor](https://en.wikipedia.org/wiki/Metric_tensor#The_round_metric_on_a_sphere)
of a sphere, the horizontal Jacobian `∂x∂ξ` is given by the diagonal matrix
`diag(R·π/180, R·cosd(lat)·π/180)`, with the determinant `J = R²·cosd(lat)·(π/180)²`.

Wording here was a bit unclear, and I thought it would be good to add a link to some reference.

Comment thread src/Grids/pointcloud.jl
diagonal matrix `diag(R·π/180, R·cosd(lat)·π/180)` reflecting the sphere metric at
each point, and `J = R²·cosd(lat)·(π/180)²`.
"""
struct PointCloudGrid{

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This could simplify all the topology(grid) methods, but subsequent operations on the result of topology would still be invalid. We could make some sort of NoTopology struct, or add methods to topology-related functions so they can handle Nothing, but that feels roughly equivalent to what this PR is doing anyway. So I don't think adding a field for the topology would be very useful.

Comment thread src/Grids/pointcloud.jl
Comment on lines +38 to +40
# No topology — return nothing. Callers that need a topology (e.g. DSS) should
# not be called on a PointCloudGrid.
topology(::PointCloudGrid) = nothing

@dennisYatunin dennisYatunin Jun 26, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could this throw an error instead? Otherwise, users will see a MethodError about how something does not work when the topology is Nothing, rather than a more informative message about how they should not access the topology of a PointCloudGrid.

Comment thread src/Grids/pointcloud.jl
device::ClimaComms.AbstractDevice = ClimaComms.device(),
context::ClimaComms.AbstractCommsContext = ClimaComms.SingletonCommsContext(device),
) where {FT}
@assert context isa ClimaComms.SingletonCommsContext "PointCloudGrid only supports SingletonCommsContext."

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
@assert context isa ClimaComms.SingletonCommsContext "PointCloudGrid only supports SingletonCommsContext."
context = ClimaComms.SingletonCommsContext(device)

As in CommonGrids, you can avoid the keyword argument if it only accepts one possible value.

Comment thread src/Spaces/multicolumn.jl

`local_geometry` holds an `IFH{LG, 1, N}` data layout.
"""
struct PointCloudLevelSpace{

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can this be called a PointCloudSpace instead? The name of every other ___Space type matches that of its corresponding ___Grid type.

Comment thread src/Spaces/multicolumn.jl
Comment on lines +164 to +165
set_mask!(fn, space::MultiColumnFiniteDifferenceSpace) =
set_mask!(fn, grid(space).horizontal_grid)

@dennisYatunin dennisYatunin Jun 26, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't this do nothing if get_mask always returns NoMask?

Comment thread src/Spaces/multicolumn.jl
Comment on lines +178 to +179
return size(level_space.local_geometry, 5) == size(
Grids.local_geometry_data(hgrid, nothing), 5)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could this check the local_geometry itself, instead of just its size? Maybe something like level_space.local_geometry === Grids.local_geometry_data(hgrid, nothing), if that's possible.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants