Feature/area extraction test

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<double>(), array::make_shape(size_halo, 2));
         ghost_        = Field("ghost", array::make_datatype<int>(), array::make_shape(size_halo));
         global_index_ = Field("global_index", array::make_datatype<gidx_t>(), array::make_shape(size_halo));

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<double>::max(), std::numeric_limits<double>::max());
-    auto max = Point2(std::numeric_limits<double>::lowest(), std::numeric_limits<double>::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<double>::max(), std::numeric_limits<double>::max());
+        auto max = Point2(std::numeric_limits<double>::lowest(), std::numeric_limits<double>::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 {

src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc

@@ -423,8 +423,10 @@ void StructuredInterpolation2D<Kernel>::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);
+            }
         }
     }
 }

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<idx_t>::max()) {
+                    lat_north = -1;
+                }
+                if (lat_south == std::numeric_limits<idx_t>::min()) {
+                    lat_south = -1;
+                }
             }
         }
     }

src/atlas/util/GridPointsJSONWriter.cc

@@ -60,6 +60,30 @@ std::vector<long> 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};
     }
 }
 

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
+} // namespace atlas

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

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<std::string>("--input-grid", default_grid);
+}
+
+int output_rank(const int& default_rank) {
+    return eckit::Resource<int>("--output-rank", default_rank);
+}
+
+std::string output_functionspace(const std::string& default_fs) {
+    return eckit::Resource<std::string>("--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<double, 2>(fs.xy());
+    for (idx_t f = 0; f < nb_fields; ++f) {
+        auto field_source = fields_source.add(fs.createField<Value>());
+        auto source       = array::make_view<Value, 2>(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<Value>(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);
+}