diff --git a/src/atlas/mesh/Mesh.h b/src/atlas/mesh/Mesh.h index 3c4e09af7..cfe643162 100644 --- a/src/atlas/mesh/Mesh.h +++ b/src/atlas/mesh/Mesh.h @@ -39,6 +39,7 @@ class PartitionPolygons; namespace atlas { namespace mesh { +class MeshBuilder; class Nodes; class HybridElements; typedef HybridElements Edges; @@ -140,6 +141,7 @@ class Mesh : DOXYGEN_HIDE(public util::ObjectHandle) { return s; } + friend class mesh::MeshBuilder; friend class meshgenerator::MeshGeneratorImpl; void setProjection(const Projection& p) { get()->setProjection(p); } void setGrid(const Grid& p) { get()->setGrid(p); } diff --git a/src/atlas/mesh/MeshBuilder.cc b/src/atlas/mesh/MeshBuilder.cc index 454e53d3e..7d1cf2f15 100644 --- a/src/atlas/mesh/MeshBuilder.cc +++ b/src/atlas/mesh/MeshBuilder.cc @@ -7,7 +7,13 @@ #include "atlas/mesh/MeshBuilder.h" +#include + #include "atlas/array/MakeView.h" +#include "atlas/grid/Grid.h" +#include "atlas/grid/Iterator.h" +#include "atlas/grid/StructuredGrid.h" +#include "atlas/grid/UnstructuredGrid.h" #include "atlas/mesh/ElementType.h" #include "atlas/mesh/HybridElements.h" #include "atlas/mesh/Mesh.h" @@ -16,6 +22,106 @@ #include "atlas/util/CoordinateEnums.h" namespace atlas { + +namespace detail { +atlas::UnstructuredGrid assemble_unstructured_grid(size_t nb_nodes, const double lons[], const double lats[], + const int ghosts[], const eckit::mpi::Comm& comm) { + + // First serialize owned lons and lats into single vector + const size_t nb_owned_nodes = std::count(ghosts, ghosts + nb_nodes, 0); + std::vector owned_lonlats(2 * nb_owned_nodes); + int counter = 0; + for (size_t i = 0; i < nb_nodes; ++i) { + if (ghosts[i] == 0) { + owned_lonlats[counter] = lons[i]; + counter++; + owned_lonlats[counter] = lats[i]; + counter++; + } + } + ATLAS_ASSERT(counter == 2 * nb_owned_nodes); + + // Gather points across MPI ranks + size_t nb_nodes_global = 0; + comm.allReduce(nb_owned_nodes, nb_nodes_global, eckit::mpi::sum()); + + std::vector global_lonlats(2 * nb_nodes_global); + eckit::mpi::Buffer buffer(comm.size()); + comm.allGatherv(owned_lonlats.begin(), owned_lonlats.end(), buffer); + global_lonlats = std::move(buffer.buffer); + + std::vector points(nb_nodes_global); + for (size_t i = 0; i < nb_nodes_global; ++i) { + points[i] = atlas::PointXY({global_lonlats[2*i], global_lonlats[2*i + 1]}); + } + + return atlas::UnstructuredGrid(new std::vector(points.begin(), points.end())); +} + +// Check if grid points match mesh points -- not obviously true as the MeshBuilder sets the Grid +// independently from the Mesh. This check can give confidence that the grid and mesh are specified +// consistently with each other. +// +// This check makes a few assumptions +// - the global_indices passed to the MeshBuilder are a 1-based contiguous index +// - the Grid is one of StructuredGrid or UnstructedGrid. This restriction is because a +// CubedSphereGrid doesn't present a simple interface for the coordinates at the N'th point. +// - the Mesh grid points are the grid nodes, and not the grid cell-centers (e.g., HEALPix or CubedSphere meshes) +void validate_grid_vs_mesh(const atlas::Grid& grid, size_t nb_nodes, const double lons[], const double lats[], + const int ghosts[], const gidx_t global_indices[], const eckit::mpi::Comm& comm) { + // Check assumption that global_indices look like a 1-based contiguous index + const size_t nb_owned_nodes = std::count(ghosts, ghosts + nb_nodes, 0); + size_t nb_nodes_global = 0; + comm.allReduce(nb_owned_nodes, nb_nodes_global, eckit::mpi::sum()); + for (size_t i = 0; i < nb_nodes; ++i) { + if (ghosts[i] == 0) { + // Check global_indices is consistent with a 1-based contiguous index over nodes + ATLAS_ASSERT(global_indices[i] >= 1); + ATLAS_ASSERT(global_indices[i] <= nb_nodes_global); + } + } + + double lonlat[2]; + const auto equal_within_roundoff = [](const double a, const double b) -> bool { + return fabs(a - b) <= 360.0 * 1.0e-16; + }; + + // Check lonlats for each supported grid type + if (grid.type() == "unstructured") { + const atlas::UnstructuredGrid ugrid(grid); + for (size_t i = 0; i < nb_nodes; ++i) { + if (ghosts[i] == 0) { + ugrid.lonlat(global_indices[i] - 1, lonlat); + if (!equal_within_roundoff(lonlat[0], lons[i]) || !equal_within_roundoff(lonlat[1], lats[i])) { + throw_Exception("In MeshBuilder: UnstructuredGrid from config does not match mesh coordinates", Here()); + } + } + } + } else if (grid.type() == "structured") { + const atlas::StructuredGrid sgrid(grid); + for (size_t i = 0; i < nb_nodes; ++i) { + if (ghosts[i] == 0) { + idx_t i, j; + sgrid.index2ij(global_indices[i] - 1, i, j); + sgrid.lonlat(i, j, lonlat); + if (!equal_within_roundoff(lonlat[0], lons[i]) || !equal_within_roundoff(lonlat[1], lats[i])) { + throw_Exception("In MeshBuilder: StructuredGrid from config does not match mesh coordinates", Here()); + } + } + } + } else { + for (size_t i = 0; i < nb_nodes; ++i) { + if (ghosts[i] == 0) { + auto point = *(grid.lonlat().begin() + static_cast(global_indices[i] - 1)); + if (!equal_within_roundoff(point.lon(), lons[i]) || !equal_within_roundoff(point.lat(), lats[i])) { + throw_Exception("In MeshBuilder: Grid from config does not match mesh coordinates", Here()); + } + } + } + } +} +} // namespace detail + namespace mesh { //---------------------------------------------------------------------------------------------------------------------- @@ -27,7 +133,8 @@ Mesh MeshBuilder::operator()(const std::vector& lons, const std::vector< const std::vector>& tri_boundary_nodes, const std::vector& tri_global_indices, const std::vector>& quad_boundary_nodes, - const std::vector& quad_global_indices) const { + const std::vector& quad_global_indices, + const eckit::Configuration& config) const { const size_t nb_nodes = global_indices.size(); const size_t nb_tris = tri_global_indices.size(); const size_t nb_quads = quad_global_indices.size(); @@ -43,16 +150,42 @@ Mesh MeshBuilder::operator()(const std::vector& lons, const std::vector< return operator()(nb_nodes, lons.data(), lats.data(), ghosts.data(), global_indices.data(), remote_indices.data(), remote_index_base, partitions.data(), nb_tris, reinterpret_cast(tri_boundary_nodes.data()), tri_global_indices.data(), nb_quads, - reinterpret_cast(quad_boundary_nodes.data()), quad_global_indices.data()); + reinterpret_cast(quad_boundary_nodes.data()), quad_global_indices.data(), + config); } Mesh MeshBuilder::operator()(size_t nb_nodes, const double lons[], const double lats[], const int ghosts[], const gidx_t global_indices[], const idx_t remote_indices[], const idx_t remote_index_base, const int partitions[], size_t nb_tris, const gidx_t tri_boundary_nodes[], const gidx_t tri_global_indices[], size_t nb_quads, const gidx_t quad_boundary_nodes[], - const gidx_t quad_global_indices[]) const { + const gidx_t quad_global_indices[], + const eckit::Configuration& config) const { + // Get MPI comm from config name or fall back to atlas default comm + auto mpi_comm_name = [](const auto& config) { + return config.getString("mpi_comm", atlas::mpi::comm().name()); + }; + const eckit::mpi::Comm& comm = eckit::mpi::comm(mpi_comm_name(config).c_str()); + Mesh mesh{}; + // Setup a grid, if requested via config argument + if (config.has("grid")) { + atlas::Grid grid; + if (config.has("grid.type") && (config.getString("grid.type") == "unstructured") + && !config.has("grid.xy")) { + // Assemble the unstructured grid by gathering input lons,lats across ranks + grid = ::atlas::detail::assemble_unstructured_grid(nb_nodes, lons, lats, ghosts, comm); + } else { + // Build grid directly from config + grid = atlas::Grid(config.getSubConfiguration("grid")); + const bool validate = config.getBool("validate", false); + if (validate) { + ::atlas::detail::validate_grid_vs_mesh(grid, nb_nodes, lons, lats, ghosts, global_indices, comm); + } + } + mesh.setGrid(grid); + } + // Populate node data mesh.nodes().resize(nb_nodes); @@ -82,11 +215,11 @@ Mesh MeshBuilder::operator()(size_t nb_nodes, const double lons[], const double // First, count how many cells of each type are on this processor // Then optimize away the element type if globally nb_tris or nb_quads is zero size_t sum_nb_tris = 0; - atlas::mpi::comm().allReduce(nb_tris, sum_nb_tris, eckit::mpi::sum()); + comm.allReduce(nb_tris, sum_nb_tris, eckit::mpi::sum()); const bool add_tris = (sum_nb_tris > 0); size_t sum_nb_quads = 0; - atlas::mpi::comm().allReduce(nb_quads, sum_nb_quads, eckit::mpi::sum()); + comm.allReduce(nb_quads, sum_nb_quads, eckit::mpi::sum()); const bool add_quads = (sum_nb_quads > 0); if (add_tris) { @@ -133,7 +266,7 @@ Mesh MeshBuilder::operator()(size_t nb_nodes, const double lons[], const double ATLAS_ASSERT(idx == nb_tris + nb_quads); - cells_part.assign(atlas::mpi::comm().rank()); + cells_part.assign(comm.rank()); return mesh; } diff --git a/src/atlas/mesh/MeshBuilder.h b/src/atlas/mesh/MeshBuilder.h index 5f2f42127..b7ab42547 100644 --- a/src/atlas/mesh/MeshBuilder.h +++ b/src/atlas/mesh/MeshBuilder.h @@ -8,8 +8,11 @@ #pragma once #include "atlas/mesh/Mesh.h" +#include "atlas/parallel/mpi/mpi.h" #include "atlas/util/Config.h" +#include "eckit/config/Configuration.h" + #include #include @@ -29,6 +32,8 @@ namespace mesh { * - cannot import halos, i.e., cells owned by other MPI tasks; halos can still be subsequently * computed by calling the BuildMesh action. * - cannot import node-to-cell connectivity information. + * + * The Mesh's Grid can be initialized via the call operator's optional config argument. */ class MeshBuilder { public: @@ -48,12 +53,21 @@ class MeshBuilder { * MPI task), in other words, must be an element of the node global_indices. The boundary nodes * are ordered node-varies-fastest, element-varies-slowest order. The cell global index is, * here also, a uniform labeling over the of the cells across all MPI tasks. + * + * The config argument can be used to + * - Request the Mesh's Grid to be constructed, usually from the config. If the Grid is either + * an UnstructuredGrid or a Structured grid, the `validate` bool option can be used to trigger + * a simple check that the grid is consistent with the lons/lats passed in to the MeshBuilder. + * In the special case where `grid.type` is unstructured and the `grid.xy` coordinates are + * _not_ given, then the grid is constructed from the lons/lats passed to the MeshBuilder. + * - Select which MPI communicator to use. */ Mesh operator()(size_t nb_nodes, const double lons[], const double lats[], const int ghosts[], const gidx_t global_indices[], const idx_t remote_indices[], const idx_t remote_index_base, const int partitions[], size_t nb_tris, const gidx_t tri_boundary_nodes[], const gidx_t tri_global_indices[], size_t nb_quads, const gidx_t quad_boundary_nodes[], - const gidx_t quad_global_indices[]) const; + const gidx_t quad_global_indices[], + const eckit::Configuration& config = util::NoConfig()) const; /** * \brief C++-interface to construct a Mesh from external connectivity data @@ -66,7 +80,8 @@ class MeshBuilder { const std::vector>& tri_boundary_nodes, const std::vector& tri_global_indices, const std::vector>& quad_boundary_nodes, - const std::vector& quad_global_indices) const; + const std::vector& quad_global_indices, + const eckit::Configuration& config = util::NoConfig()) const; }; //----------------------------------------------------------------------------- diff --git a/src/tests/mesh/test_mesh_builder.cc b/src/tests/mesh/test_mesh_builder.cc index 91f236000..67dcadecb 100644 --- a/src/tests/mesh/test_mesh_builder.cc +++ b/src/tests/mesh/test_mesh_builder.cc @@ -10,11 +10,13 @@ #include #include "atlas/array/MakeView.h" +#include "atlas/grid.h" #include "atlas/mesh/Elements.h" #include "atlas/mesh/HybridElements.h" #include "atlas/mesh/Mesh.h" #include "atlas/mesh/MeshBuilder.h" #include "atlas/mesh/Nodes.h" +#include "atlas/util/Config.h" #include "tests/AtlasTestEnvironment.h" #include "tests/mesh/helper_mesh_builder.h" @@ -129,15 +131,78 @@ CASE("test_cs_c2_mesh_serial") { std::iota(quad_global_indices.begin(), quad_global_indices.end(), 9); // nb_tris + 1 const MeshBuilder mesh_builder{}; - const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions, - tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices); - - helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base, - partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, - quad_global_indices); - //Gmsh gmsh("out.msh", util::Config("coordinates", "xyz")); - //gmsh.write(mesh); + SECTION("Build Mesh without a Grid") { + const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions, + tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices); + + helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base, + partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, + quad_global_indices); + EXPECT(!mesh.grid()); + + //Gmsh gmsh("out.msh", util::Config("coordinates", "xyz")); + //gmsh.write(mesh); + } + + SECTION("Build Mesh with an UnstructuredGrid from config") { + std::vector lonlats(2 * lons.size()); + int counter = 0; + for (size_t i = 0; i < lons.size(); ++i) { + lonlats[counter] = lons[i]; + counter++; + lonlats[counter] = lats[i]; + counter++; + } + + util::Config config{}; + config.set("grid.type", "unstructured"); + config.set("grid.xy", lonlats); + config.set("validate", true); + + const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions, + tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices, + config); + + helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base, + partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, + quad_global_indices); + + EXPECT(mesh.grid()); + EXPECT(mesh.grid().type() == "unstructured"); + } + + SECTION("Build Mesh with an UnstructuredGrid, with Grid assembled from MeshBuilder arguments") { + util::Config config{}; + config.set("grid.type", "unstructured"); + + const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions, + tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices, + config); + + helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base, + partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, + quad_global_indices); + + EXPECT(mesh.grid()); + EXPECT(mesh.grid().type() == "unstructured"); + } + + SECTION("Build Mesh with a CubedSphereGrid from config") { // but can't validate this + util::Config config{}; + config.set("grid.name", "CS-LFR-2"); + + const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions, + tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices, + config); + + helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base, + partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, + quad_global_indices); + + EXPECT(mesh.grid()); + EXPECT(mesh.grid().type() == "cubedsphere"); + } } //----------------------------------------------------------------------------- diff --git a/src/tests/mesh/test_mesh_builder_parallel.cc b/src/tests/mesh/test_mesh_builder_parallel.cc index fd4e0c7b3..6d824a284 100644 --- a/src/tests/mesh/test_mesh_builder_parallel.cc +++ b/src/tests/mesh/test_mesh_builder_parallel.cc @@ -152,15 +152,69 @@ CASE("test_cs_c2_mesh_parallel") { [&global_partitions](const gidx_t idx) { return global_partitions[idx - 1]; }); const MeshBuilder mesh_builder{}; - const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions, - tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices); - helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base, - partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, - quad_global_indices); + SECTION("Build Mesh without a Grid") { + const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions, + tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices); - //Gmsh gmsh("out.msh", util::Config("coordinates", "xyz")); - //gmsh.write(mesh); + helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base, + partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, + quad_global_indices); + EXPECT(!mesh.grid()); + + //Gmsh gmsh("out.msh", util::Config("coordinates", "xyz")); + //gmsh.write(mesh); + } + + SECTION("Build Mesh with an UnstructuredGrid from config") { + const size_t nb_owned_nodes = std::count(ghosts.begin(), ghosts.end(), 0); + std::vector owned_lonlats(2 * nb_owned_nodes); + int counter = 0; + for (size_t i = 0; i < lons.size(); ++i) { + if (ghosts[i] == 0) { + owned_lonlats[counter] = lons[i]; + counter++; + owned_lonlats[counter] = lats[i]; + counter++; + } + } + size_t nb_nodes_global = 0; + mpi::comm().allReduce(nb_owned_nodes, nb_nodes_global, eckit::mpi::sum()); + std::vector global_lonlats(2 * nb_nodes_global); + eckit::mpi::Buffer buffer(mpi::comm().size()); + mpi::comm().allGatherv(owned_lonlats.begin(), owned_lonlats.end(), buffer); + global_lonlats = std::move(buffer.buffer); + + util::Config config{}; + config.set("grid.type", "unstructured"); + config.set("grid.xy", global_lonlats); + config.set("validate", true); + + const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions, + tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices, + config); + + helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base, + partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, + quad_global_indices); + EXPECT(mesh.grid()); + EXPECT(mesh.grid().type() == "unstructured"); + } + + SECTION("Build Mesh with an UnstructuredGrid, with Grid assembled from MeshBuilder arguments") { + util::Config config{}; + config.set("grid.type", "unstructured"); + + const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions, + tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices, + config); + + helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base, + partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, + quad_global_indices); + EXPECT(mesh.grid()); + EXPECT(mesh.grid().type() == "unstructured"); + } } //-----------------------------------------------------------------------------