Interpolate from cubedsphere mesh to structured columns

[odlomax]

Hi @wdeconinck,

I've been working on a test to demonstrate functionality that will help @MayeulDestouches with some of his SABER work. It demonstrates how to interpolate from an arbitrarily partitioned CubedSphere mesh (and possibly others) to a StructuredColumns functionspace. This is historically quite tricky, due to the constraints on placed on StructuredColumns. I basically cheat and use a couple of intermediate PointCloud functionspaces to sort out all the redistribution stuff.

One of the things I needed to write was this lambda:

    const auto makePointCloud = [](const Grid& grid, const grid::Partitioner partitioner) {

        const auto distribution = grid::Distribution(grid, partitioner);

        auto lonLats = std::vector<PointXY>{};
        auto idx = gidx_t{0};
        for (const auto& lonLat : grid.lonlat()) {
            if (distribution.partition(idx++) == mpi::rank()) {
                lonLats.emplace_back(lonLat.data());
            }
        }
        return functionspace::PointCloud(lonLats);
    };

Seems to me like there should be a constructor in PointCloud for this. Perhaps two: a Partitioner version and a Distribution version.

Do you agree? If so, it's pretty easy for me to add them.

[odlomax]

@ytremolet @fmahebert

FYI, this is how you can do a global interpolation with Atlas. It's even simpler if it's a NodeColumns to NodeColumns interpolation as you don't have to make the intermediate PointCloud functionspaces.

[odlomax]

Hi @wdeconinck,

This PR fixes a crash that @MayeulDestouches was experiencing when using the "cubedsphere" matching mesh partitioner in SABER. @MayeulDestouches, could you add a link to the issue here.

Cheers!

[wdeconinck]

@odlomax thanks for this. The only problem seems to have been the "listSize"?

I agree with a PointCloud constructor as you proposed as well, this could be in a following PR.

[MayeulDestouches]

The bug solved by this PR is documented here. I copy it here in case people don't have access to it:

As another limitation of this three-step atlas interpolation, code will fail on 1, 2 and 3 PEs for some combination of grids (eg from CS-LFR-12 to F15). For any CS-LFR grid, we can find a Gaussian grid fine enough to make the cubed-sphere matching mesh partitioning fail.

[odlomax]

@odlomax thanks for this. The only problem seems to have been the "listSize"?

I agree with a PointCloud constructor as you proposed as well, this could be in a following PR.

That's correct. Ideally, I'd like a cubed-sphere mesh generator that creates non-disjoint meshes. Then we could use the standard matching mesh tools. A dream for Cubed-Sphere 2.0 perhaps...