diff --git a/src/atlas/functionspace/PointCloud.cc b/src/atlas/functionspace/PointCloud.cc index 2ce1d64c7..25594c800 100644 --- a/src/atlas/functionspace/PointCloud.cc +++ b/src/atlas/functionspace/PointCloud.cc @@ -17,18 +17,22 @@ #include "atlas/array.h" #include "atlas/field/Field.h" #include "atlas/field/FieldSet.h" +#include "atlas/grid/Distribution.h" #include "atlas/grid/Grid.h" #include "atlas/grid/Iterator.h" #include "atlas/option/Options.h" #include "atlas/parallel/HaloExchange.h" #include "atlas/parallel/mpi/mpi.h" #include "atlas/runtime/Exception.h" +#include "atlas/util/KDTree.h" #include "atlas/util/CoordinateEnums.h" #include "atlas/util/Metadata.h" #include "atlas/util/Point.h" #include "atlas/util/Unique.h" + #include "eckit/mpi/Comm.h" +#include "eckit/log/Bytes.h" #if ATLAS_HAVE_FORTRAN #define REMOTE_IDX_BASE 1 @@ -88,18 +92,141 @@ PointCloud::PointCloud(const FieldSet & flds): lonlat_(flds["lonlat"]) { } } -PointCloud::PointCloud(const Grid& grid) { - lonlat_ = Field("lonlat", array::make_datatype(), array::make_shape(grid.size(), 2)); - auto lonlat = array::make_view(lonlat_); +grid::Partitioner make_partitioner(const Grid& grid, const eckit::Configuration& config) { + auto partitioner = grid.partitioner(); + if( config.has("partitioner") ) { + partitioner.set("type",config.getString("partitioner")); + } + if( not partitioner.has("type") ) { + partitioner.set("type","equal_regions"); + } + return grid::Partitioner(partitioner); +} + +PointCloud::PointCloud(const Grid& grid, const eckit::Configuration& config) : + PointCloud(grid, make_partitioner(grid,config), config) { +} + +PointCloud::PointCloud(const Grid& grid, const grid::Partitioner& _partitioner, const eckit::Configuration& config) { + ATLAS_TRACE("PointCloud(grid,partitioner,config)"); + int part = mpi::rank(); + double halo_radius; + config.get("halo_radius", halo_radius = 0.); + + grid::Partitioner partitioner(_partitioner); + if ( not partitioner ) { + partitioner = grid::Partitioner("equal_regions"); + } + + auto distribution = partitioner.partition(grid); + auto size_owned = distribution.nb_pts()[part]; + size_owned_ = size_owned; + + if (halo_radius == 0. || mpi::size() == 1) { + idx_t size_halo = size_owned; + ATLAS_ASSERT(size_owned > 0); + lonlat_ = Field("lonlat", array::make_datatype(), array::make_shape(size_halo, 2)); + ghost_ = Field("ghost", array::make_datatype(), array::make_shape(size_halo)); + global_index_ = Field("global_index", array::make_datatype(), array::make_shape(size_halo)); + partition_ = Field("partition", array::make_datatype(), array::make_shape(size_halo)); + remote_index_ = Field("remote_index", array::make_datatype(), array::make_shape(size_halo)); + + auto lonlat = array::make_view(lonlat_); + auto ridx = array::make_indexview(remote_index_); + auto gidx = array::make_view(global_index_); + array::make_view(ghost_).assign(0); + array::make_view(partition_).assign(part); + + idx_t j{0}; + gidx_t g{0}; + for (auto p : grid.lonlat()) { + if( distribution.partition(g++) == part ) { + gidx(j) = g+1; + ridx(j) = j; + lonlat(j, 0) = p.lon(); + lonlat(j, 1) = p.lat(); + ++j; + } + } + } + else { + + std::vector keep; + + { + ATLAS_TRACE("Build list of points to keep"); + std::vector owned_lonlat; + std::vector owned_grid_idx; + owned_lonlat.reserve(size_owned); + owned_grid_idx.reserve(size_owned); + + auto kdtree = util::IndexKDTree(); + { + ATLAS_TRACE("build kdtree"); + kdtree.reserve(grid.size()); + idx_t j{0}; + for (auto p : grid.lonlat()) { + if( distribution.partition(j) == part ) { + owned_lonlat.emplace_back(p); + owned_grid_idx.emplace_back(j); + } + kdtree.insert(p,j); + ++j; + } + kdtree.build(); + } + + keep.resize(grid.size()); + { + ATLAS_TRACE("search kdtree"); + for (idx_t j=0; j(), array::make_shape(size_halo, 2)); + partition_ = Field("partition", array::make_datatype(), array::make_shape(size_halo)); + ghost_ = Field("ghost", array::make_datatype(), array::make_shape(size_halo)); + global_index_ = Field("global_index", array::make_datatype(), array::make_shape(size_halo)); + max_glb_idx_ = grid.size(); + auto lonlat = array::make_view(lonlat_); + auto partition = array::make_view(partition_); + auto ghost = array::make_view(ghost_); + auto glb_idx = array::make_view(global_index_); + + gidx_t g=0; + idx_t j = 0; + for (auto p : grid.lonlat()) { + if (keep[g]) { + lonlat(j, 0) = p.lon(); + lonlat(j, 1) = p.lat(); + partition(j) = distribution.partition(g); + ghost(j) = partition(j) != part; + glb_idx(j) = g+1; + ++j; + } + ++g; + } + + } - idx_t j{0}; - for (auto p : grid.lonlat()) { - lonlat(j, 0) = p.lon(); - lonlat(j, 1) = p.lat(); - ++j; } + + setupHaloExchange(); } + Field PointCloud::ghost() const { if (not ghost_) { ghost_ = Field("ghost", array::make_datatype(), array::make_shape(size())); @@ -316,124 +443,384 @@ void PointCloud::haloExchange(const Field& field, bool on_device) const { haloExchange(fieldset, on_device); } -void PointCloud::setupHaloExchange(){ +void PointCloud::create_remote_index() const { + ATLAS_TRACE(); + const eckit::mpi::Comm& comm = atlas::mpi::comm(); + const int mpi_rank = comm.rank(); + const int mpi_size = comm.size(); + auto size_halo = lonlat_.shape(0); + + remote_index_ = Field("remote_idx", array::make_datatype(), array::make_shape(size_halo)); + auto remote_idx = array::make_indexview(remote_index_); + auto ghost = array::make_view(ghost_); + auto part = array::make_view(partition_); + + atlas_omp_parallel_for(idx_t n = 0; n < size_halo; ++n) { + if(not ghost(n)) { + remote_idx(n) = n; + } + else { + remote_idx(n) = -1; + } + } + + auto build_partition_graph = [this,mpi_size,mpi_rank,size_halo,&comm]() -> std::unique_ptr { + + auto part = array::make_view(this->partition_); + auto ghost = array::make_view(this->ghost_); + + std::vector others_set(mpi_size, 0); + others_set[mpi_rank] = 1; + for( idx_t i = 0; i others; + others.reserve(mpi_size); + for (idx_t p = 0; p < mpi_size; ++p) { + if (others_set[p]) { + others.emplace_back(p); + } + } + + eckit::mpi::Buffer recv_others(mpi_size); + + ATLAS_TRACE_MPI(ALLGATHER) { comm.allGatherv(others.begin(), others.end(), recv_others); } + + std::vector counts(recv_others.counts.begin(), recv_others.counts.end()); + std::vector displs(recv_others.displs.begin(), recv_others.displs.end()); + std::vector values(recv_others.buffer.begin(), recv_others.buffer.end()); + return std::unique_ptr( + new Mesh::PartitionGraph(values.data(), mpi_size, displs.data(), counts.data())); + }; + + std::unique_ptr graph_ptr; + ATLAS_TRACE_SCOPE("Building partition graph...") { graph_ptr = build_partition_graph(); } + const Mesh::PartitionGraph& graph = *graph_ptr; + + ATLAS_TRACE_SCOPE("Setup remote_index fields...") { + auto p = array::make_view(partition_); + auto g = array::make_view(global_index_); + + auto neighbours = graph.nearestNeighbours(mpi_rank); + const idx_t nb_neighbours = static_cast(neighbours.size()); + std::unordered_map part_to_neighbour; + ATLAS_TRACE_SCOPE("part_to_neighbour") { + part_to_neighbour.reserve(nb_neighbours); + for (idx_t j = 0; j < nb_neighbours; ++j) { + part_to_neighbour[neighbours[j]] = j; + } + } + std::vector halo_per_neighbour(neighbours.size(), 0); + ATLAS_TRACE_SCOPE("set halo_per_neighbour") + for( idx_t i = 0; i> g_per_neighbour(neighbours.size()); + ATLAS_TRACE_SCOPE("assemble g_per_neighbour") { + for (idx_t j = 0; j < nb_neighbours; ++j) { + g_per_neighbour[j].reserve(halo_per_neighbour[j]); + } + for( idx_t j = 0; j send_requests(neighbours.size()); + std::vector recv_requests(neighbours.size()); + + std::vector recv_size(neighbours.size()); + std::vector send_size(neighbours.size()); + + int tag = 0; + ATLAS_TRACE_SCOPE("send-receive g_per_neighbour size") { + for (idx_t j = 0; j < nb_neighbours; ++j) { + send_size[j] = static_cast(g_per_neighbour[j].size()); + send_requests[j] = comm.iSend(send_size[j], neighbours[j], tag); + recv_requests[j] = comm.iReceive(recv_size[j], neighbours[j], tag); + } + + ATLAS_TRACE_MPI(WAIT) { + for (idx_t j = 0; j < nb_neighbours; ++j) { + comm.wait(send_requests[j]); + } + + for (idx_t j = 0; j < nb_neighbours; ++j) { + comm.wait(recv_requests[j]); + } + } + } + + std::vector> recv_g_per_neighbour(neighbours.size()); + ATLAS_TRACE_SCOPE("send-receive g_per_neighbour") + for (idx_t j = 0; j < nb_neighbours; ++j) { + recv_g_per_neighbour[j].resize(recv_size[j]); + + send_requests[j] = comm.iSend(g_per_neighbour[j].data(), g_per_neighbour[j].size(), neighbours[j], tag); + recv_requests[j] = + comm.iReceive(recv_g_per_neighbour[j].data(), recv_g_per_neighbour[j].size(), neighbours[j], tag); + } + + std::vector> send_r_per_neighbour(neighbours.size()); + + // Assemble "send_r_per_neighbour" + // Depending if we can afford memory for a globally sized array, we can have a + // much faster version of g_to_r map using std::vector. + // TODO: using c++14 we can create a polymorphic lambda to avoid duplicated + // code in the two branches of following if. + + // size_t max_glb_idx = grid().size(); + ATLAS_ASSERT(max_glb_idx_); + + atlas::vector g_to_r_vector; + bool use_unordered_map_fallback = false; + try { + g_to_r_vector.resize(max_glb_idx_ + 1); + } + catch (std::bad_alloc& e) { + if (comm.size() > 1) { + Log::warning() << "Could not allocate " << eckit::Bytes((max_glb_idx_ + 1) * sizeof(idx_t)) << Here() + << "\n" + << "Using slower unordered_map fallback to map global to remote indices" + << std::endl; + use_unordered_map_fallback = true; + } + else { + throw_Exception( + "Could not allocate " + std::string(eckit::Bytes((max_glb_idx_ + 1) * sizeof(idx_t))), Here()); + } + } + if (not use_unordered_map_fallback) { + auto& g_to_r = g_to_r_vector; + ATLAS_TRACE_SCOPE("g_to_r (using vector)") { + atlas_omp_parallel_for(idx_t j = 0; j < size_halo; ++j) { + if (not ghost(j)) { + // parallel omp possible for ``` g_to_r[g(j)] = j ``` as we only loop over size_owned, + // where global_index is such that race-conditions cannot occur +#if ATLAS_ARRAYVIEW_BOUNDS_CHECKING + if (g(j) >= max_glb_idx_ + 1) { + ATLAS_DEBUG_VAR(g(j)); + throw_OutOfRange("g_to_r", g(j), max_glb_idx_ + 1, Here()); + } +#endif + g_to_r[g(j)] = j; + } + } + } + for (idx_t j = 0; j < nb_neighbours; ++j) { + send_r_per_neighbour[j].reserve(recv_size[j]); + + comm.wait(recv_requests[j]); // wait for recv_g_per_neighbour[j] + for (gidx_t gidx : recv_g_per_neighbour[j]) { + send_r_per_neighbour[j].emplace_back(g_to_r[gidx]); + } + } + } + else { + std::unordered_map g_to_r; + g_to_r.reserve(size_halo); + + ATLAS_TRACE_SCOPE("g_to_r (using unordered set)") { + for (idx_t j = 0; j < size_halo; ++j) { + if (not ghost(j)) { + g_to_r[g(j)] = j; + } + } + } + ATLAS_TRACE_SCOPE("assemble r_per_neighbour") + for (idx_t j = 0; j < nb_neighbours; ++j) { + send_r_per_neighbour[j].reserve(recv_size[j]); + + comm.wait(recv_requests[j]); // wait for recv_g_per_neighbour[j] + for (gidx_t gidx : recv_g_per_neighbour[j]) { + send_r_per_neighbour[j].emplace_back(g_to_r[gidx]); + } + } + } + + std::vector> r_per_neighbour(neighbours.size()); + + ATLAS_TRACE_SCOPE("send-receive r_per_neighbour") { + for (idx_t j = 0; j < nb_neighbours; ++j) { + r_per_neighbour[j].resize(halo_per_neighbour[j]); + } + + ATLAS_TRACE_MPI(ALLTOALL) { + for (idx_t j = 0; j < nb_neighbours; ++j) { + comm.wait(send_requests[j]); + send_requests[j] = comm.iSend(send_r_per_neighbour[j].data(), send_r_per_neighbour[j].size(), + neighbours[j], tag); + recv_requests[j] = + comm.iReceive(r_per_neighbour[j].data(), r_per_neighbour[j].size(), neighbours[j], tag); + } + } + ATLAS_TRACE_MPI(WAIT) { + for (idx_t j = 0; j < nb_neighbours; ++j) { + comm.wait(recv_requests[j]); + } + } + } + + + std::vector counters(neighbours.size(), 0); + ATLAS_TRACE_SCOPE("set remote_idx") + for (idx_t j = 0; j < size_halo; ++j) { + if (ghost(j)) { + idx_t neighbour = part_to_neighbour[p(j)]; + remote_idx(j) = r_per_neighbour[neighbour][counters[neighbour]++]; + } + ATLAS_ASSERT(remote_idx(j) >= 0); + } + + ATLAS_TRACE_MPI(WAIT) { + for (idx_t j = 0; j < nb_neighbours; ++j) { + comm.wait(send_requests[j]); + } + } + } +} + +void PointCloud::setupHaloExchange() { + ATLAS_TRACE(); const eckit::mpi::Comm& comm = atlas::mpi::comm(); - std::size_t mpi_rank = comm.rank(); - const std::size_t mpi_size = comm.size(); - - if (not partition_ and not remote_index_) { - - auto lonlat_v = array::make_view(lonlat_); - // data structure containing a flag to identify the 'ghost points'; - // 0={is not a ghost point}, 1={is a ghost point} - auto is_ghost = array::make_view(ghost_); - - std::vector opoints_local; - std::vector gpoints_local; - std::vector lonlat_u; - std::vector opoints_local_u; - - for (idx_t i = 0; i < lonlat_v.shape(0); ++i){ - lonlat_u.emplace_back(util::unique_lonlat(lonlat_v(i, XX), lonlat_v(i, YY))); - } - - idx_t j {0}; - for (idx_t i = 0; i < is_ghost.shape(0); ++i) { - PointXY loc(lonlat_v(j, XX), lonlat_v(j, YY)); - if (is_ghost(i)) { - gpoints_local.emplace_back(loc); - } else { - opoints_local.emplace_back(loc); - opoints_local_u.emplace_back(util::unique_lonlat(loc.x(), loc.y())); - } - ++j; - } - - std::vector coords_gp_local; - coords_gp_local.reserve(gpoints_local.size() * 2); - - for (auto& gp : gpoints_local) { - coords_gp_local.push_back(gp[XX]); - coords_gp_local.push_back(gp[YY]); - } - - eckit::mpi::Buffer buffers_rec(mpi_size); - - comm.allGatherv(coords_gp_local.begin(), coords_gp_local.end(), buffers_rec); - - std::vector gpoints_global; - - for (std::size_t pe = 0; pe < mpi_size; ++pe) { - for (std::size_t j = 0; j < buffers_rec.counts[pe]/2; ++j) { - PointXY loc_gp(*(buffers_rec.begin() + buffers_rec.displs[pe] + 2 * j + XX), - *(buffers_rec.begin() + buffers_rec.displs[pe] + 2 * j + YY)); - gpoints_global.emplace_back(loc_gp); - } - } - - std::vector gpoints_global_u; - for (atlas::PointXY& loc : gpoints_global) { - gpoints_global_u.emplace_back(util::unique_lonlat(loc.x(), loc.y())); - } - - std::vector partition_ids_gp_global(gpoints_global.size(), -1); - std::vector remote_index_gp_global(gpoints_global.size(), -1); - - std::vector::iterator iter_xy_gp_01; - - for (std::size_t idx = 0; idx < gpoints_global_u.size(); ++idx) { - iter_xy_gp_01 = std::find(opoints_local_u.begin(), - opoints_local_u.end(), gpoints_global_u.at(idx)); - if (iter_xy_gp_01 != opoints_local_u.end()) { - std::size_t ridx = std::distance(opoints_local_u.begin(), iter_xy_gp_01); - partition_ids_gp_global.at(idx) = mpi_rank; - remote_index_gp_global.at(idx) = ridx; - } - } - - comm.allReduceInPlace(partition_ids_gp_global.begin(), - partition_ids_gp_global.end(), eckit::mpi::max()); - comm.allReduceInPlace(remote_index_gp_global.begin(), - remote_index_gp_global.end(), eckit::mpi::max()); - - std::vector partition_ids_local(lonlat_v.shape(0), -1); - std::vector remote_index_local(lonlat_v.shape(0), -1); - - idx_t idx_loc {0}; - std::vector::iterator iter_xy_gp_02; - - for (idx_t i = 0; i < lonlat_v.shape(0); ++i){ - iter_xy_gp_02 = std::find(gpoints_global_u.begin(), - gpoints_global_u.end(), lonlat_u.at(i)); - if (iter_xy_gp_02 != gpoints_global_u.end()) { - std::size_t idx_gp = std::distance(gpoints_global_u.begin(), iter_xy_gp_02); - partition_ids_local[idx_loc] = partition_ids_gp_global[idx_gp]; - remote_index_local[idx_loc] = remote_index_gp_global[idx_gp]; - } else { - partition_ids_local[idx_loc] = mpi_rank; - remote_index_local[idx_loc] = idx_loc; - } - ++idx_loc; - } - - partition_ = Field("partition", array::make_datatype(), - array::make_shape(partition_ids_local.size())); - - auto partitionv = array::make_view(partition_); - for (idx_t i = 0; i < partitionv.shape(0); ++i) { - partitionv(i) = partition_ids_local.at(i); - } - - remote_index_ = Field("remote_index", array::make_datatype(), - array::make_shape(remote_index_local.size())); - - auto remote_indexv = array::make_indexview(remote_index_); - for (idx_t i = 0; i < remote_indexv.shape(0); ++i) { - remote_indexv(i) = remote_index_local.at(i); - } + const int mpi_rank = comm.rank(); + const int mpi_size = comm.size(); + + + if (ghost_ and partition_ and global_index_ and not remote_index_) { + create_remote_index(); + } + else if (not partition_ or not remote_index_) { + ATLAS_TRACE("do setup"); + + auto lonlat_v = array::make_view(lonlat_); + // data structure containing a flag to identify the 'ghost points'; + // 0={is not a ghost point}, 1={is a ghost point} + auto is_ghost = array::make_view(ghost_); + + std::vector opoints_local; + std::vector gpoints_local; + std::vector lonlat_u; + std::vector opoints_local_u; + lonlat_u.reserve(lonlat_v.shape(0)); + opoints_local_u.reserve(is_ghost.shape(0)); + gpoints_local.reserve(is_ghost.shape(0)); + opoints_local.reserve(is_ghost.shape(0)); + + for (idx_t i = 0; i < lonlat_v.shape(0); ++i) { + lonlat_u.emplace_back(util::unique_lonlat(lonlat_v(i, XX), lonlat_v(i, YY))); + } + + idx_t j{0}; + for (idx_t i = 0; i < is_ghost.shape(0); ++i) { + PointXY loc(lonlat_v(j, XX), lonlat_v(j, YY)); + if (is_ghost(i)) { + gpoints_local.emplace_back(loc); + } + else { + opoints_local.emplace_back(loc); + opoints_local_u.emplace_back(util::unique_lonlat(loc.x(), loc.y())); + } + ++j; + } + + std::vector coords_gp_local; + coords_gp_local.reserve(gpoints_local.size() * 2); + + for (auto& gp : gpoints_local) { + coords_gp_local.push_back(gp[XX]); + coords_gp_local.push_back(gp[YY]); + } + + eckit::mpi::Buffer buffers_rec(mpi_size); + + ATLAS_TRACE_SCOPE("mpi all gather") { + comm.allGatherv(coords_gp_local.begin(), coords_gp_local.end(), buffers_rec); + } + + std::vector gpoints_global; + + for (std::size_t pe = 0; pe < mpi_size; ++pe) { + for (std::size_t j = 0; j < buffers_rec.counts[pe] / 2; ++j) { + PointXY loc_gp(*(buffers_rec.begin() + buffers_rec.displs[pe] + 2 * j + XX), + *(buffers_rec.begin() + buffers_rec.displs[pe] + 2 * j + YY)); + gpoints_global.emplace_back(loc_gp); + } + } + + std::vector gpoints_global_u; + gpoints_global_u.reserve(gpoints_global.size()); + for (atlas::PointXY& loc : gpoints_global) { + gpoints_global_u.emplace_back(util::unique_lonlat(loc.x(), loc.y())); + } + + std::vector partition_ids_gp_global(gpoints_global.size(), -1); + std::vector remote_index_gp_global(gpoints_global.size(), -1); + + std::vector::iterator iter_xy_gp_01; + ATLAS_TRACE_SCOPE("find global") + for (std::size_t idx = 0; idx < gpoints_global_u.size(); ++idx) { + iter_xy_gp_01 = std::find(opoints_local_u.begin(), opoints_local_u.end(), gpoints_global_u.at(idx)); + if (iter_xy_gp_01 != opoints_local_u.end()) { + std::size_t ridx = std::distance(opoints_local_u.begin(), iter_xy_gp_01); + partition_ids_gp_global.at(idx) = mpi_rank; + remote_index_gp_global.at(idx) = ridx; + } + } + + ATLAS_TRACE_SCOPE("mpi all reduce") { + comm.allReduceInPlace(partition_ids_gp_global.begin(), partition_ids_gp_global.end(), eckit::mpi::max()); + comm.allReduceInPlace(remote_index_gp_global.begin(), remote_index_gp_global.end(), eckit::mpi::max()); + } + + std::vector partition_ids_local(lonlat_v.shape(0), -1); + std::vector remote_index_local(lonlat_v.shape(0), -1); + + idx_t idx_loc{0}; + std::vector::iterator iter_xy_gp_02; + + ATLAS_TRACE_SCOPE("find local") + for (idx_t i = 0; i < lonlat_v.shape(0); ++i) { + iter_xy_gp_02 = std::find(gpoints_global_u.begin(), gpoints_global_u.end(), lonlat_u.at(i)); + if (iter_xy_gp_02 != gpoints_global_u.end()) { + std::size_t idx_gp = std::distance(gpoints_global_u.begin(), iter_xy_gp_02); + partition_ids_local[idx_loc] = partition_ids_gp_global[idx_gp]; + remote_index_local[idx_loc] = remote_index_gp_global[idx_gp]; + } + else { + partition_ids_local[idx_loc] = mpi_rank; + remote_index_local[idx_loc] = idx_loc; + } + ++idx_loc; + } + + partition_ = Field("partition", array::make_datatype(), array::make_shape(partition_ids_local.size())); + + auto partitionv = array::make_view(partition_); + for (idx_t i = 0; i < partitionv.shape(0); ++i) { + partitionv(i) = partition_ids_local.at(i); + } + + remote_index_ = + Field("remote_index", array::make_datatype(), array::make_shape(remote_index_local.size())); + + auto remote_indexv = array::make_indexview(remote_index_); + for (idx_t i = 0; i < remote_indexv.shape(0); ++i) { + remote_indexv(i) = remote_index_local.at(i); + } } + ATLAS_ASSERT(partition_); + ATLAS_ASSERT(ghost_); + ATLAS_ASSERT(remote_index_); ATLAS_ASSERT(ghost_.size() == remote_index_.size()); ATLAS_ASSERT(ghost_.size() == partition_.size()); @@ -502,8 +889,11 @@ PointCloud::PointCloud(const std::initializer_list : new detail::PointCloud{std::vector(points.begin(), points.end())})), functionspace_(dynamic_cast(get())) {} -PointCloud::PointCloud(const Grid& grid): - FunctionSpace(new detail::PointCloud(grid)), functionspace_(dynamic_cast(get())) {} +PointCloud::PointCloud(const Grid& grid, const eckit::Configuration& config): + FunctionSpace(new detail::PointCloud(grid, config)), functionspace_(dynamic_cast(get())) {} + +PointCloud::PointCloud(const Grid& grid, const grid::Partitioner& partitioner, const eckit::Configuration& config): + FunctionSpace(new detail::PointCloud(grid, partitioner,config)), functionspace_(dynamic_cast(get())) {} } // namespace functionspace diff --git a/src/atlas/functionspace/PointCloud.h b/src/atlas/functionspace/PointCloud.h index 7c0e14ccd..a6affc4e3 100644 --- a/src/atlas/functionspace/PointCloud.h +++ b/src/atlas/functionspace/PointCloud.h @@ -23,6 +23,7 @@ #include "atlas/runtime/Exception.h" #include "atlas/util/Config.h" #include "atlas/util/Point.h" +#include "atlas/grid/Partitioner.h" namespace atlas { namespace parallel { @@ -46,7 +47,8 @@ class PointCloud : public functionspace::FunctionSpaceImpl { PointCloud(const Field& lonlat); PointCloud(const Field& lonlat, const Field& ghost); PointCloud(const FieldSet&); // assuming lonlat ghost ridx and partition present - PointCloud(const Grid&); + PointCloud(const Grid&, const eckit::Configuration& = util::NoConfig()); + PointCloud(const Grid&, const grid::Partitioner&, const eckit::Configuration& = util::NoConfig()); ~PointCloud() override {} std::string type() const override { return "PointCloud"; } operator bool() const override { return true; } @@ -149,14 +151,17 @@ class PointCloud : public functionspace::FunctionSpaceImpl { void set_field_metadata(const eckit::Configuration& config, Field& field) const; + void create_remote_index() const; private: Field lonlat_; Field vertical_; mutable Field ghost_; - Field remote_index_; + mutable Field remote_index_; Field global_index_; Field partition_; + idx_t size_owned_; + idx_t max_glb_idx_{0}; std::unique_ptr halo_exchange_; idx_t levels_{0}; @@ -179,7 +184,8 @@ class PointCloud : public FunctionSpace { PointCloud(const std::vector&); PointCloud(const std::vector&); PointCloud(const std::initializer_list>&); - PointCloud(const Grid&); + PointCloud(const Grid&, const eckit::Configuration& = util::NoConfig()); + PointCloud(const Grid&, const grid::Partitioner&, const eckit::Configuration& = util::NoConfig()); operator bool() const { return valid(); } bool valid() const { return functionspace_; } diff --git a/src/atlas/grid/detail/grid/Unstructured.cc b/src/atlas/grid/detail/grid/Unstructured.cc index f5047a756..258436b6e 100644 --- a/src/atlas/grid/detail/grid/Unstructured.cc +++ b/src/atlas/grid/detail/grid/Unstructured.cc @@ -221,7 +221,7 @@ Grid::Config Unstructured::meshgenerator() const { return Config("type", "delaunay"); } Grid::Config Unstructured::partitioner() const { - return Config("type", "serial"); + return Config("type", "equal_regions"); } void Unstructured::print(std::ostream& os) const { diff --git a/src/tests/functionspace/CMakeLists.txt b/src/tests/functionspace/CMakeLists.txt index 0859497ee..137478faf 100644 --- a/src/tests/functionspace/CMakeLists.txt +++ b/src/tests/functionspace/CMakeLists.txt @@ -92,6 +92,15 @@ ecbuild_add_test( TARGET atlas_test_pointcloud CONDITION eckit_HAVE_MPI AND NOT HAVE_GRIDTOOLS_STORAGE ) + +ecbuild_add_test( TARGET atlas_test_pointcloud_halo_creation + SOURCES test_pointcloud_halo_creation.cc + LIBS atlas + MPI 4 + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} + CONDITION eckit_HAVE_MPI AND NOT HAVE_GRIDTOOLS_STORAGE AND MPI_SLOTS GREATER_EQUAL 4 +) + ecbuild_add_test( TARGET atlas_test_pointcloud_he_2PE SOURCES test_pointcloud_haloexchange_2PE.cc LIBS atlas diff --git a/src/tests/functionspace/fctest_functionspace.F90 b/src/tests/functionspace/fctest_functionspace.F90 index a705790a6..fd2011d11 100644 --- a/src/tests/functionspace/fctest_functionspace.F90 +++ b/src/tests/functionspace/fctest_functionspace.F90 @@ -628,8 +628,8 @@ module fcta_FunctionSpace_fxt field2 = fs%create_field(name="field2",kind=atlas_real(8),levels=3,variables=2) -FCTEST_CHECK_EQUAL( field%shape(), ([5,int(grid%size())]) ) -FCTEST_CHECK_EQUAL( field2%shape(), ([2,3,int(grid%size())]) ) +FCTEST_CHECK_EQUAL( field%shape(), ([5,int(fs%size())]) ) +FCTEST_CHECK_EQUAL( field2%shape(), ([2,3,int(fs%size())]) ) #ifndef _CRAYFTN FCTEST_CHECK_EQUAL( field%owners(), 1 ) diff --git a/src/tests/functionspace/test_pointcloud_halo_creation.cc b/src/tests/functionspace/test_pointcloud_halo_creation.cc new file mode 100644 index 000000000..6228c2a34 --- /dev/null +++ b/src/tests/functionspace/test_pointcloud_halo_creation.cc @@ -0,0 +1,269 @@ +/* + * (C) Copyright 2013 ECMWF + * (C) Crown Copyright 2023 Met Office + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * + */ + + +#include + +#include "atlas/array.h" +#include "atlas/field.h" +#include "atlas/functionspace/PointCloud.h" +#include "atlas/option.h" +#include "atlas/parallel/mpi/mpi.h" + +#include "tests/AtlasTestEnvironment.h" + +#include "atlas/grid.h" +#include "atlas/util/KDTree.h" + +using namespace eckit; +using namespace atlas::functionspace; +using namespace atlas::util; + +namespace atlas { +namespace test { + +//----------------------------------------------------------------------------- + +CASE("Distributed creation from unstructured grid with halo") { +std::vector points = { + {180,0}, + {90,0}, + {-90,0}, + {0,90}, + {0,-90}, + {0,0}, + {18,0}, + {36,0}, + {54,0}, + {72,0}, + {108,0}, + {126,0}, + {144,0}, + {162,0}, + {-162,0}, + {-144,0}, + {-126,0}, + {-108,0}, + {-72,0}, + {-54,0}, + {-36,0}, + {-18,0}, + {0,18}, + {0,36}, + {0,54}, + {0,72}, + {180,72}, + {180,54}, + {180,36}, + {180,18}, + {180,-18}, + {180,-36}, + {180,-54}, + {180,-72}, + {0,-72}, + {0,-54}, + {0,-36}, + {0,-18}, + {90,18}, + {90,36}, + {90,54}, + {90,72}, + {-90,72}, + {-90,54}, + {-90,36}, + {-90,18}, + {-90,-18}, + {-90,-36}, + {-90,-54}, + {-90,-72}, + {90,-72}, + {90,-54}, + {90,-36}, + {90,-18}, + {123.974,-58.6741}, + {154.087,-16.9547}, + {154.212,-58.8675}, + {114.377,-41.9617}, + {125.567,-23.5133}, + {137.627,-40.8524}, + {106.162,-24.5874}, + {158.508,-38.55}, + {137.826,-72.8109}, + {142.103,-26.799}, + {138.256,-13.8871}, + {168.39,-24.3266}, + {168.954,-12.0094}, + {117.333,-12.35}, + {102.254,-11.1537}, + {120.307,59.7167}, + {107.196,26.0167}, + {144.768,28.3721}, + {150.891,60.0343}, + {164.566,25.5053}, + {116.851,14.0295}, + {124.84,28.3978}, + {157.901,42.042}, + {111.41,43.1056}, + {134.333,44.6677}, + {103.277,11.707}, + {135.358,73.2119}, + {135.349,14.2311}, + {153.48,13.386}, + {168.071,11.5344}, + {-162.99,26.3775}, + {-147.519,56.1313}, + {-122.579,27.4824}, + {-117.909,59.2376}, + {-104.052,27.3616}, + {-153.107,14.9717}, + {-110.833,41.7436}, + {-144.847,32.8534}, + {-161.546,42.1031}, + {-129.866,44.5201}, + {-133.883,72.4163}, + {-166.729,11.8907}, + {-135.755,15.2529}, + {-106.063,14.4869}, + {-119.452,11.7037}, + {-146.026,-58.6741}, + {-115.913,-16.9547}, + {-115.788,-58.8675}, + {-155.623,-41.9617}, + {-144.433,-23.5133}, + {-132.373,-40.8524}, + {-163.838,-24.5874}, + {-111.492,-38.55}, + {-132.174,-72.8109}, + {-127.897,-26.799}, + {-131.744,-13.8871}, + {-101.61,-24.3266}, + {-101.046,-12.0094}, + {-152.667,-12.35}, + {-167.746,-11.1537}, + {-14.0127,-27.2963}, + {-59.193,-57.0815}, + {-56.465,-19.5751}, + {-27.056,-59.3077}, + {-57.124,-35.9752}, + {-33.4636,-28.3914}, + {-74.8037,-46.8602}, + {-40.089,-45.1376}, + {-74.8149,-28.3136}, + {-21.3072,-42.2177}, + {-44.0778,-72.6353}, + {-19.6969,-12.8527}, + {-40.1318,-12.1601}, + {-72.691,-11.4129}, + {-56.0261,58.6741}, + {-25.9127,16.9547}, + {-25.7876,58.8675}, + {-65.6229,41.9617}, + {-54.4335,23.5133}, + {-42.373,40.8524}, + {-73.838,24.5874}, + {-21.4917,38.55}, + {-42.1744,72.8109}, + {-37.8974,26.799}, + {-41.7437,13.8871}, + {-11.6095,24.3266}, + {-11.0459,12.0094}, + {-62.667,12.35}, + {-77.7456,11.1537}, + {30.3071,59.7167}, + {17.1956,26.0167}, + {54.7676,28.3721}, + {60.8915,60.0343}, + {74.5657,25.5053}, + {26.8506,14.0295}, + {34.8398,28.3978}, + {67.9014,42.042}, + {21.41,43.1056}, + {44.3335,44.6677}, + {13.2772,11.707}, + {45.3579,73.2119}, + {45.3492,14.2311}, + {63.4799,13.386}, + {78.0714,11.5344}, + {17.01,-26.3775}, + {32.4806,-56.1313}, + {57.4213,-27.4824}, + {62.0912,-59.2376}, + {75.9483,-27.3616}, + {26.893,-14.9717}, + {69.1672,-41.7436}, + {35.1527,-32.8534}, + {18.4543,-42.1031}, + {50.1344,-44.5201}, + {46.1172,-72.4163}, + {13.2711,-11.8907}, + {44.2448,-15.2529}, + {73.9368,-14.4869}, + {60.5478,-11.7037} +}; +auto grid = UnstructuredGrid(points); + +auto pointcloud = functionspace::PointCloud(grid,util::Config("halo_radius",2000000)); + +auto field = pointcloud.createField(option::name("field")); + +auto lonlat = array::make_view(pointcloud.lonlat()); +auto ghost = array::make_view(pointcloud.ghost()); + +auto view = array::make_view(field); + +size_t count_ghost{0}; +for (idx_t i=0; i