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)); 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 { 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); + } } } } 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; + } } } } 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 diff --git a/src/tests/interpolation/CMakeLists.txt b/src/tests/interpolation/CMakeLists.txt index 29d220590..296c834b9 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 4 + 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..17d794496 --- /dev/null +++ b/src/tests/interpolation/test_interpolation_structured2D_to_area.cc @@ -0,0 +1,167 @@ +/* + * (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/functionspace/PointCloud.h" +#include "atlas/grid/Partitioner.h" + +#include "atlas/util/GridPointsJSONWriter.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); +} + +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; + 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 outrank = output_rank(-1); + idx_t nb_fields = 2; + idx_t nb_levels = 19; + + // Setup Grid and functionspace + 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); + + // 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); + + 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"); + + 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); + + // setup target field_set + FieldSet fields_target = create_target_fields(outputFS, nb_fields, nb_levels); + + // execute interpolation + interpolation.execute(fields_source, fields_target); + + if (atlas::functionspace::NodeColumns(outputFS)) { + output::Gmsh gmsh{"area.msh","a"}; + gmsh.write(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); +}