From 5a83e47257b7f6c2a6850cbc98922416594e2f28 Mon Sep 17 00:00:00 2001 From: Antonino Bonanni Date: Thu, 7 Mar 2024 16:31:22 +0000 Subject: [PATCH 1/9] area cutout unit-test --- src/tests/interpolation/CMakeLists.txt | 8 ++ ...test_interpolation_structured2D_to_area.cc | 113 ++++++++++++++++++ 2 files changed, 121 insertions(+) create mode 100644 src/tests/interpolation/test_interpolation_structured2D_to_area.cc diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index 29d220590..1ec6d1f63 100644 --- a/src/tests/interpolation/CMakeLists.txt +++ b/src/tests/interpolation/CMakeLists.txt @@ -113,6 +113,14 @@ ecbuild_add_test( TARGET atlas_test_interpolation_structured2D_to_points ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) +ecbuild_add_test( TARGET atlas_test_interpolation_structured2D_to_area + SOURCES test_interpolation_structured2D_to_area.cc + LIBS atlas + MPI 2 + CONDITION eckit_HAVE_MPI + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) + ecbuild_add_test( TARGET atlas_test_interpolation_cubedsphere SOURCES test_interpolation_cubedsphere.cc diff --git a/src/tests/interpolation/test_interpolation_structured2D_to_area.cc b/src/tests/interpolation/test_interpolation_structured2D_to_area.cc new file mode 100644 index 000000000..f3ee35818 --- /dev/null +++ b/src/tests/interpolation/test_interpolation_structured2D_to_area.cc @@ -0,0 +1,113 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * 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. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include "atlas/array.h" +#include "atlas/field/Field.h" +#include "atlas/field/FieldSet.h" + +#include "atlas/functionspace/StructuredColumns.h" +#include "atlas/grid/Grid.h" +#include "atlas/grid/Iterator.h" +#include "atlas/interpolation.h" +#include "atlas/mesh/Mesh.h" +#include "atlas/meshgenerator.h" +#include "atlas/output/Gmsh.h" +#include "atlas/util/CoordinateEnums.h" +#include "atlas/util/function/VortexRollup.h" +#include "atlas/util/PolygonXY.h" + +#include "atlas/domain/Domain.h" +#include "atlas/parallel/mpi/mpi.h" +#include "atlas/functionspace/NodeColumns.h" +#include "atlas/grid/Partitioner.h" + +#include "tests/AtlasTestEnvironment.h" + +using atlas::util::Config; + +namespace atlas { +namespace test { + +//----------------------------------------------------------------------------- + +std::string input_gridname(const std::string& default_grid) { + return eckit::Resource("--input-grid", default_grid); +} + +FieldSet create_source_fields(functionspace::StructuredColumns& fs, idx_t nb_fields, idx_t nb_levels) { + using Value = double; + FieldSet fields_source; + auto lonlat = array::make_view(fs.xy()); + for (idx_t f = 0; f < nb_fields; ++f) { + auto field_source = fields_source.add(fs.createField()); + auto source = array::make_view(field_source); + for (idx_t n = 0; n < fs.size(); ++n) { + for (idx_t k = 0; k < nb_levels; ++k) { + source(n, k) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); + } + }; + } + return fields_source; +} +FieldSet create_target_fields(FunctionSpace& fs, idx_t nb_fields, idx_t nb_levels) { + using Value = double; + FieldSet fields_target; + for (idx_t f = 0; f < nb_fields; ++f) { + fields_target.add(fs.createField(option::levels(nb_levels))); + } + return fields_target; +} + +void do_test( std::string type, int input_halo, bool expect_failure ) { + idx_t nb_fields = 2; + idx_t nb_levels = 137; + + std::cout << "======>> mpi::comm().rank(): " << mpi::comm().rank() << std::endl; + + // Setup Grid and functionspace + Grid inputGrid(input_gridname("O400")); + functionspace::StructuredColumns inputFS(inputGrid, option::levels(nb_levels) | option::halo(input_halo)); + + // Setup source field_set + FieldSet fields_source = create_source_fields(inputFS, nb_fields, nb_levels); + + // Define cutout area and grid + double boundingBoxNorth = 45; + double boundingBoxSouth = -45; + double boundingBoxEast = 140; + double boundingBoxWest = 50; + + RectangularLonLatDomain rd({boundingBoxWest, boundingBoxEast}, {boundingBoxSouth, boundingBoxNorth}); + Grid areaGrid(inputGrid, rd); + + Mesh areaMesh = MeshGenerator("structured").generate( areaGrid, grid::MatchingPartitioner(inputFS) ); + atlas::FunctionSpace outputFS = atlas::functionspace::NodeColumns{areaMesh}; + + // setup interpolation + Interpolation interpolation(option::type(type), inputFS, outputFS); + + // setup target field_set + FieldSet fields_target = create_target_fields(outputFS, nb_fields, nb_levels); + + // execute interpolation + interpolation.execute(fields_source, fields_target); + +} + +CASE("test structured-linear2D, halo 3") { + EXPECT_NO_THROW( do_test("structured-linear2D", 3, false) ); +} + +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { + return atlas::test::run(argc, argv); +} From af247645fd255b0f25725733eb9026b929637a73 Mon Sep 17 00:00:00 2001 From: Antonino Bonanni Date: Fri, 8 Mar 2024 08:32:13 +0000 Subject: [PATCH 2/9] area cutout unit test: set -np=4 and nlev=19 --- src/tests/interpolation/CMakeLists.txt | 2 +- .../interpolation/test_interpolation_structured2D_to_area.cc | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index 1ec6d1f63..296c834b9 100644 --- a/src/tests/interpolation/CMakeLists.txt +++ b/src/tests/interpolation/CMakeLists.txt @@ -116,7 +116,7 @@ ecbuild_add_test( TARGET atlas_test_interpolation_structured2D_to_points ecbuild_add_test( TARGET atlas_test_interpolation_structured2D_to_area SOURCES test_interpolation_structured2D_to_area.cc LIBS atlas - MPI 2 + MPI 4 CONDITION eckit_HAVE_MPI ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) diff --git a/src/tests/interpolation/test_interpolation_structured2D_to_area.cc b/src/tests/interpolation/test_interpolation_structured2D_to_area.cc index f3ee35818..27682f67e 100644 --- a/src/tests/interpolation/test_interpolation_structured2D_to_area.cc +++ b/src/tests/interpolation/test_interpolation_structured2D_to_area.cc @@ -67,9 +67,7 @@ FieldSet create_target_fields(FunctionSpace& fs, idx_t nb_fields, idx_t nb_level void do_test( std::string type, int input_halo, bool expect_failure ) { idx_t nb_fields = 2; - idx_t nb_levels = 137; - - std::cout << "======>> mpi::comm().rank(): " << mpi::comm().rank() << std::endl; + idx_t nb_levels = 19; // Setup Grid and functionspace Grid inputGrid(input_gridname("O400")); From b23036043cebf52f86e6b69b272a24f3c6276620 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 8 Mar 2024 10:58:16 +0100 Subject: [PATCH 3/9] Fix StructuredMeshGeneration situation with multiple OpenMP threads --- src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc index 96da93938..c89ed6bf6 100644 --- a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc @@ -360,6 +360,12 @@ Find min and max latitudes used by this part. } lat_north = *std::min_element(thread_reduce_lat_north.begin(), thread_reduce_lat_north.end()); lat_south = *std::max_element(thread_reduce_lat_south.begin(), thread_reduce_lat_south.end()); + if (lat_north == std::numeric_limits::max()) { + lat_north = -1; + } + if (lat_south == std::numeric_limits::min()) { + lat_south = -1; + } } } } From 763fe37d1e50c8459dc6443a00a2aa5611a11452 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 8 Mar 2024 11:00:31 +0100 Subject: [PATCH 4/9] Fix StructuredInterpolation2D using matrix to partitions without points --- .../method/structured/StructuredInterpolation2D.tcc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc index b40436f97..7e559647b 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc @@ -423,8 +423,10 @@ void StructuredInterpolation2D::setup( const FunctionSpace& source ) { // fill sparse matrix if( failed_points.empty() ) { idx_t inp_npts = source.size(); - Matrix A( out_npts_, inp_npts, triplets ); - setMatrix(A); + if (out_npts_ > 0) { + Matrix A( out_npts_, inp_npts, triplets ); + setMatrix(A); + } } } } From 83b3bf066790f7015b6ebb6a5ce34fc8f417fb68 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 8 Mar 2024 11:01:37 +0100 Subject: [PATCH 5/9] Reduce requirements on area cutout unit-test and output results --- .../test_interpolation_structured2D_to_area.cc | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/tests/interpolation/test_interpolation_structured2D_to_area.cc b/src/tests/interpolation/test_interpolation_structured2D_to_area.cc index 27682f67e..c7f1b3207 100644 --- a/src/tests/interpolation/test_interpolation_structured2D_to_area.cc +++ b/src/tests/interpolation/test_interpolation_structured2D_to_area.cc @@ -70,7 +70,7 @@ void do_test( std::string type, int input_halo, bool expect_failure ) { idx_t nb_levels = 19; // Setup Grid and functionspace - Grid inputGrid(input_gridname("O400")); + Grid inputGrid(input_gridname("O40")); functionspace::StructuredColumns inputFS(inputGrid, option::levels(nb_levels) | option::halo(input_halo)); // Setup source field_set @@ -88,8 +88,11 @@ void do_test( std::string type, int input_halo, bool expect_failure ) { Mesh areaMesh = MeshGenerator("structured").generate( areaGrid, grid::MatchingPartitioner(inputFS) ); atlas::FunctionSpace outputFS = atlas::functionspace::NodeColumns{areaMesh}; + output::Gmsh gmsh{"area.msh"}; + gmsh.write(areaMesh); + // setup interpolation - Interpolation interpolation(option::type(type), inputFS, outputFS); + Interpolation interpolation(option::type(type)|Config("matrix_free",false), inputFS, outputFS); // setup target field_set FieldSet fields_target = create_target_fields(outputFS, nb_fields, nb_levels); @@ -97,6 +100,7 @@ void do_test( std::string type, int input_halo, bool expect_failure ) { // execute interpolation interpolation.execute(fields_source, fields_target); + gmsh.write(fields_target); } CASE("test structured-linear2D, halo 3") { From 102a734458d7433c07050e4b3a859b82a270ae8c Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 11 Mar 2024 10:29:36 +0000 Subject: [PATCH 6/9] Remove assertion in PointCloud to cater for empty partitions --- src/atlas/functionspace/PointCloud.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/src/atlas/functionspace/PointCloud.cc b/src/atlas/functionspace/PointCloud.cc index 61388a2f5..4f2dd4c71 100644 --- a/src/atlas/functionspace/PointCloud.cc +++ b/src/atlas/functionspace/PointCloud.cc @@ -170,7 +170,6 @@ PointCloud::PointCloud(const Grid& grid, const grid::Partitioner& _partitioner, if (halo_radius == 0. || nb_partitions_ == 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)); From 3c6fb6ab5507bda66f3a597a49645e147377e2d1 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 11 Mar 2024 14:40:40 +0000 Subject: [PATCH 7/9] Robustness --- src/atlas/grid/StructuredPartitionPolygon.cc | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/atlas/grid/StructuredPartitionPolygon.cc b/src/atlas/grid/StructuredPartitionPolygon.cc index 4a52d5948..ab7929146 100644 --- a/src/atlas/grid/StructuredPartitionPolygon.cc +++ b/src/atlas/grid/StructuredPartitionPolygon.cc @@ -35,6 +35,11 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect const auto grid = fs.grid(); const auto dom = RectangularDomain(grid.domain()); + if ( fs.size() == 0 ) { + bb.clear(); + return; + } + if (_halo > 0 && _halo != fs.halo()) { throw_Exception("halo must zero or matching the one of the StructuredColumns", Here()); } @@ -324,13 +329,15 @@ StructuredPartitionPolygon::StructuredPartitionPolygon(const functionspace::Func fs_(fs), halo_(halo) { ATLAS_TRACE("StructuredPartitionPolygon"); compute(fs, halo, points_, inner_bounding_box_); - auto min = Point2(std::numeric_limits::max(), std::numeric_limits::max()); - auto max = Point2(std::numeric_limits::lowest(), std::numeric_limits::lowest()); - for (size_t i = 0; i < inner_bounding_box_.size() - 1; ++i) { - min = Point2::componentsMin(min, inner_bounding_box_[i]); - max = Point2::componentsMax(max, inner_bounding_box_[i]); + if (inner_bounding_box_.size()) { + auto min = Point2(std::numeric_limits::max(), std::numeric_limits::max()); + auto max = Point2(std::numeric_limits::lowest(), std::numeric_limits::lowest()); + for (size_t i = 0; i < inner_bounding_box_.size() - 1; ++i) { + min = Point2::componentsMin(min, inner_bounding_box_[i]); + max = Point2::componentsMax(max, inner_bounding_box_[i]); + } + inscribed_domain_ = {{min[XX], max[XX]}, {min[YY], max[YY]}}; } - inscribed_domain_ = {{min[XX], max[XX]}, {min[YY], max[YY]}}; } size_t StructuredPartitionPolygon::footprint() const { From dc35fed8c7783edb01e3c5f5ee7a67e81523a312 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 11 Mar 2024 14:43:02 +0000 Subject: [PATCH 8/9] Update GridPointsJSONWriter to accept also an already created Partitioner --- src/atlas/util/GridPointsJSONWriter.cc | 26 +++++++++++++++++++++++++- src/atlas/util/GridPointsJSONWriter.h | 6 ++++-- 2 files changed, 29 insertions(+), 3 deletions(-) diff --git a/src/atlas/util/GridPointsJSONWriter.cc b/src/atlas/util/GridPointsJSONWriter.cc index 7888a0c9a..042985bca 100644 --- a/src/atlas/util/GridPointsJSONWriter.cc +++ b/src/atlas/util/GridPointsJSONWriter.cc @@ -60,6 +60,30 @@ std::vector points_from_list(const std::string& list, long base) { //------------------------------------------------------------------------------ +GridPointsJSONWriter::GridPointsJSONWriter(Grid grid, grid::Partitioner partitioner, const eckit::Parametrisation& args) : grid_{grid} { + args.get("json.precision",precision_=-1); + args.get("verbose",verbose_=0); + if (not args.get("partition",partition_=-1)) { + args.get("partition",partition_); + } + args.get("json.pretty", pretty_=false); + args.get("field",field_="lonlat"); + args.get("field_base",field_base_=0); + std::string points_list; + if (args.get("index",points_list)) { + args.get("index_base",points_base_ = 0); + points_ = points_from_list(points_list,points_base_); + } + + nb_partitions_ = partitioner.nb_partitions(); + ATLAS_DEBUG_VAR(nb_partitions_); + if( nb_partitions_ > 0 ) { + distribution_ = grid::Distribution{grid_, partitioner}; + } +} + +//------------------------------------------------------------------------------ + GridPointsJSONWriter::GridPointsJSONWriter(Grid grid, const eckit::Parametrisation& args) : grid_{grid} { args.get("json.precision",precision_=-1); args.get("verbose",verbose_=0); @@ -85,7 +109,7 @@ GridPointsJSONWriter::GridPointsJSONWriter(Grid grid, const eckit::Parametrisati partitioner_config.set("type", partitioner); } partitioner_config.set("partitions", nb_partitions_); - distribution_ = grid::Distribution{grid_,partitioner_config}; + distribution_ = grid::Distribution{grid_, partitioner_config}; } } diff --git a/src/atlas/util/GridPointsJSONWriter.h b/src/atlas/util/GridPointsJSONWriter.h index 0164c6169..4901d226b 100644 --- a/src/atlas/util/GridPointsJSONWriter.h +++ b/src/atlas/util/GridPointsJSONWriter.h @@ -22,7 +22,9 @@ namespace util { class GridPointsJSONWriter { public: - GridPointsJSONWriter(Grid grid, const eckit::Parametrisation& args); + GridPointsJSONWriter(Grid, const eckit::Parametrisation&); + + GridPointsJSONWriter(Grid, grid::Partitioner, const eckit::Parametrisation&); void write(std::ostream& out, eckit::Channel& info) const; @@ -45,4 +47,4 @@ class GridPointsJSONWriter { //------------------------------------------------------------------------------ } // namespace util -} // namespace atlas \ No newline at end of file +} // namespace atlas From c1dd11bd5c1a73e09244984a3ce81e53dd7ca5c1 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 12 Mar 2024 16:23:09 +0000 Subject: [PATCH 9/9] Add option to change output functionspace in atlas_test_interpolation_structured2D_to_area --- ...test_interpolation_structured2D_to_area.cc | 62 +++++++++++++++++-- 1 file changed, 57 insertions(+), 5 deletions(-) diff --git a/src/tests/interpolation/test_interpolation_structured2D_to_area.cc b/src/tests/interpolation/test_interpolation_structured2D_to_area.cc index c7f1b3207..17d794496 100644 --- a/src/tests/interpolation/test_interpolation_structured2D_to_area.cc +++ b/src/tests/interpolation/test_interpolation_structured2D_to_area.cc @@ -26,8 +26,11 @@ #include "atlas/domain/Domain.h" #include "atlas/parallel/mpi/mpi.h" #include "atlas/functionspace/NodeColumns.h" +#include "atlas/functionspace/PointCloud.h" #include "atlas/grid/Partitioner.h" +#include "atlas/util/GridPointsJSONWriter.h" + #include "tests/AtlasTestEnvironment.h" using atlas::util::Config; @@ -41,6 +44,15 @@ std::string input_gridname(const std::string& default_grid) { return eckit::Resource("--input-grid", default_grid); } +int output_rank(const int& default_rank) { + return eckit::Resource("--output-rank", default_rank); +} + +std::string output_functionspace(const std::string& default_fs) { + return eckit::Resource("--output-functionspace", default_fs); +} + + FieldSet create_source_fields(functionspace::StructuredColumns& fs, idx_t nb_fields, idx_t nb_levels) { using Value = double; FieldSet fields_source; @@ -66,6 +78,7 @@ FieldSet create_target_fields(FunctionSpace& fs, idx_t nb_fields, idx_t nb_level } void do_test( std::string type, int input_halo, bool expect_failure ) { + idx_t outrank = output_rank(-1); idx_t nb_fields = 2; idx_t nb_levels = 19; @@ -73,6 +86,8 @@ void do_test( std::string type, int input_halo, bool expect_failure ) { Grid inputGrid(input_gridname("O40")); functionspace::StructuredColumns inputFS(inputGrid, option::levels(nb_levels) | option::halo(input_halo)); + inputFS.polygon(0).outputPythonScript("input.py"); + // Setup source field_set FieldSet fields_source = create_source_fields(inputFS, nb_fields, nb_levels); @@ -85,11 +100,45 @@ void do_test( std::string type, int input_halo, bool expect_failure ) { RectangularLonLatDomain rd({boundingBoxWest, boundingBoxEast}, {boundingBoxSouth, boundingBoxNorth}); Grid areaGrid(inputGrid, rd); - Mesh areaMesh = MeshGenerator("structured").generate( areaGrid, grid::MatchingPartitioner(inputFS) ); - atlas::FunctionSpace outputFS = atlas::functionspace::NodeColumns{areaMesh}; + util::GridPointsJSONWriter input_writer(inputGrid, + util::Config + ("partition",mpi::rank()) + ("partitioner.partitions",mpi::size()) + ("field","lonlat") + ); + if( mpi::rank() == outrank ) { + Log::info() << "input grid coordinates of rank " << outrank << ": \n"; + input_writer.write(Log::info()); + } + + util::GridPointsJSONWriter output_writer(areaGrid, grid::MatchingPartitioner(inputFS), + util::Config + ("partition",mpi::rank()) + ("field","lonlat") + ); + if( mpi::rank() == outrank ) { + Log::info() << "output grid coordinates of rank " << outrank << ": \n"; + output_writer.write(Log::info()); + } + + FunctionSpace outputFS; + std::string output_functionspace_type = output_functionspace("NodeColumns"); - output::Gmsh gmsh{"area.msh"}; - gmsh.write(areaMesh); + if (output_functionspace_type=="NodeColumns") { + Mesh areaMesh = MeshGenerator("structured").generate( areaGrid, grid::MatchingPartitioner(inputFS) ); + output::Gmsh gmsh{"area.msh"}; + gmsh.write(areaMesh); + outputFS = atlas::functionspace::NodeColumns{areaMesh}; + } + else if (output_functionspace_type=="StructuredColumns") { + outputFS = functionspace::StructuredColumns(areaGrid, grid::MatchingPartitioner(inputFS) ); + } + else if (output_functionspace_type=="PointCloud") { + outputFS = atlas::functionspace::PointCloud{areaGrid, grid::MatchingPartitioner(inputFS)}; + } + else { + ATLAS_NOTIMPLEMENTED; + } // setup interpolation Interpolation interpolation(option::type(type)|Config("matrix_free",false), inputFS, outputFS); @@ -100,7 +149,10 @@ void do_test( std::string type, int input_halo, bool expect_failure ) { // execute interpolation interpolation.execute(fields_source, fields_target); - gmsh.write(fields_target); + if (atlas::functionspace::NodeColumns(outputFS)) { + output::Gmsh gmsh{"area.msh","a"}; + gmsh.write(fields_target); + } } CASE("test structured-linear2D, halo 3") {