Merge branch 'release/0.15.2'

* release/0.15.2: (21 commits)
  Update changelog
  version bump
  Too strict tolerance for test_trans_pole
  fixing L-grid domains on southern hemisphere (MIR-283)
  Fix TransIFS with Gridtools-CUDA backend
  removed float comparison which might have caused issues at the equator
  fixed missing equator (MIR-285) and issue at poles also updated the unit tests
  Silenced warning
  Running cuda 9.1 on lxg
  Fix wrong mod from previous commit
  ATLAS-164 fix failing unit tests on cray
  Remove PGI compiler warnings
  ATLAS-162 Enable unstructured test
  ATLAS-171 Another unit-test
  ATLAS-171 Implement cropped unstructured grids
  ATLAS-173 Work on structuredcolumns indexing
  correct typo
  MIR-282 Massage default options for TransLocal to deal better with unstructured grids
  MIR-283 Add unit test
  MIR-283 Fix remaining issues
  ...

CHANGELOG.md

@@ -7,6 +7,19 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html
 
 ## [Unreleased]
 
+## [0.15.2] - 2018-08-31
+### Changed
+- Initialisation of Fields to signalling NaN in debug builds, uninitialised in
+  non-debug builds (used to be initialised to zero as part of std::vector construction)
+
+### Added
+- Implementation of cropped unstructured grids so that spectral transforms to
+  unstructured grids are allowed
+
+### Fixed
+- Spectral transforms to grids including pole and equator
+- Build with gridtools CUDA backend
+
 ## [0.15.1] - 2018-07-17
 ### Fixed
 - Compilation for Intel 18 debug
@@ -45,6 +58,7 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html
 ## 0.13.0 - 2018-02-16
 
 [Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop
+[0.15.2]: https://github.com/ecmwf/atlas/compare/0.15.1...0.15.2
 [0.15.1]: https://github.com/ecmwf/atlas/compare/0.15.0...0.15.1
 [0.15.0]: https://github.com/ecmwf/atlas/compare/0.14.0...0.15.0
 [0.14.0]: https://github.com/ecmwf/atlas/compare/0.13.2...0.14.0

CMakeLists.txt

@@ -169,7 +169,7 @@ ecbuild_add_option( FEATURE ACC
 
 if( ATLAS_HAVE_ACC )
   if( CMAKE_Fortran_COMPILER_ID MATCHES "PGI" )
-    set( ACC_Fortran_FLAGS -acc -ta=tesla,cuda${CUDA_VERSION},nordc )
+    set( ACC_Fortran_FLAGS -acc -ta=tesla,nordc )
     set( ACC_C_FLAGS ${ACC_Fortran_FLAGS} )
     find_program( ACC_C_COMPILER NAMES pgcc HINTS ${PGI_DIR} ENV PGI_DIR PATH_SUFFIXES bin )
     if( NOT ACC_C_COMPILER )

VERSION.cmake

@@ -6,5 +6,5 @@
 # granted to it by virtue of its status as an intergovernmental organisation nor
 # does it submit to any jurisdiction.
 
-set  ( ${PROJECT_NAME}_VERSION_STR  "0.15.1" )
+set  ( ${PROJECT_NAME}_VERSION_STR  "0.15.2" )
 

src/atlas/domain/Domain.cc

@@ -20,7 +20,7 @@ using ZD = atlas::domain::ZonalBandDomain;
 
 namespace atlas {
 
-Domain::Domain() : domain_( new domain::EmptyDomain() ) {}
+Domain::Domain() : domain_() {}
 
 Domain::Domain( const Domain& domain ) : domain_( domain.domain_ ) {}
 

src/atlas/domain/Domain.h

@@ -44,7 +44,7 @@ class Domain {
     Domain( const Implementation* );
     Domain( const eckit::Parametrisation& );
 
-    operator bool() { return true; }
+    operator bool() const { return domain_; }
 
     /// Type of the domain
     std::string type() const;

src/atlas/functionspace/StructuredColumns.cc

@@ -412,10 +412,14 @@ StructuredColumns::StructuredColumns( const Grid& grid, const grid::Partitioner&
     field_partition_    = Field( "partition", array::make_datatype<int>(), array::make_shape( size_halo_ ) );
     field_global_index_ = Field( "glb_idx", array::make_datatype<gidx_t>(), array::make_shape( size_halo_ ) );
     field_remote_index_ = Field( "remote_idx", array::make_datatype<idx_t>(), array::make_shape( size_halo_ ) );
+    field_index_i_      = Field( "index_i", array::make_datatype<idx_t>(), array::make_shape( size_halo_ ) );
+    field_index_j_      = Field( "index_j", array::make_datatype<idx_t>(), array::make_shape( size_halo_ ) );
 
     auto part       = array::make_view<int, 1>( field_partition_ );
     auto global_idx = array::make_view<gidx_t, 1>( field_global_index_ );
     auto remote_idx = array::make_view<idx_t, 1>( field_remote_index_ );
+    auto index_i    = array::make_indexview<idx_t, 1>( field_index_i_ );
+    auto index_j    = array::make_indexview<idx_t, 1>( field_index_j_ );
 
     for ( const GridPoint& gp : gridpoints ) {
         bool in_domain( false );
@@ -432,6 +436,8 @@ StructuredColumns::StructuredColumns( const Grid& grid, const grid::Partitioner&
             global_idx( gp.r ) = compute_g( gp.i, gp.j );
             part( gp.r )       = compute_p( gp.i, gp.j );
         }
+        index_i( gp.r ) = gp.i;
+        index_j( gp.r ) = gp.j;
     }
 
     ATLAS_TRACE_SCOPE( "Parallelisation ..." ) {
@@ -944,6 +950,14 @@ field::FieldImpl* atlas__fs__StructuredColumns__partition( const detail::Structu
 field::FieldImpl* atlas__fs__StructuredColumns__global_index( const detail::StructuredColumns* This ) {
     return This->global_index().get();
 }
+
+field::FieldImpl* atlas__fs__StructuredColumns__index_i( const detail::StructuredColumns* This ) {
+    return This->index_i().get();
+}
+
+field::FieldImpl* atlas__fs__StructuredColumns__index_j( const detail::StructuredColumns* This ) {
+    return This->index_j().get();
+}
 }
 // ----------------------------------------------------------------------------
 

src/atlas/functionspace/StructuredColumns.h

@@ -95,6 +95,8 @@ class StructuredColumns : public FunctionSpaceImpl {
     Field partition() const { return field_partition_; }
     Field global_index() const { return field_global_index_; }
     Field remote_index() const { return field_remote_index_; }
+    Field index_i() const { return field_index_i_; }
+    Field index_j() const { return field_index_j_; }
 
 private:  // methods
     size_t config_size( const eckit::Configuration& config ) const;
@@ -121,6 +123,8 @@ class StructuredColumns : public FunctionSpaceImpl {
     Field field_partition_;
     Field field_global_index_;
     Field field_remote_index_;
+    Field field_index_i_;
+    Field field_index_j_;
 
     class Map2to1 {
     public:
@@ -303,6 +307,8 @@ int atlas__fs__StructuredColumns__i_end_halo( const detail::StructuredColumns* T
 field::FieldImpl* atlas__fs__StructuredColumns__xy( const detail::StructuredColumns* This );
 field::FieldImpl* atlas__fs__StructuredColumns__partition( const detail::StructuredColumns* This );
 field::FieldImpl* atlas__fs__StructuredColumns__global_index( const detail::StructuredColumns* This );
+field::FieldImpl* atlas__fs__StructuredColumns__index_i( const detail::StructuredColumns* This );
+field::FieldImpl* atlas__fs__StructuredColumns__index_j( const detail::StructuredColumns* This );
 }
 
 }  // namespace functionspace

src/atlas/grid/Grid.cc

@@ -34,7 +34,9 @@ Grid::Grid( const Grid& grid ) : grid_( grid.grid_ ) {}
 Grid::Grid( const Grid::Implementation* grid ) : grid_( grid ) {}
 
 Grid::Grid( const std::string& shortname, const Domain& domain ) {
-    grid_ = Grid::Implementation::create( shortname, Config( "domain", domain.spec() ) );
+    Config config;
+    if ( domain ) config.set( "domain", domain.spec() );
+    grid_ = Grid::Implementation::create( shortname, config );
 }
 
 Grid::Grid( const Grid& grid, const Grid::Domain& domain ) {
@@ -74,6 +76,10 @@ UnstructuredGrid::UnstructuredGrid( std::initializer_list<PointXY> xy ) :
     Grid( new UnstructuredGrid::grid_t( xy ) ),
     grid_( unstructured_grid( get() ) ) {}
 
+UnstructuredGrid::UnstructuredGrid( const Grid& grid, const Grid::Domain& domain ) :
+    Grid( new UnstructuredGrid::grid_t( *grid.get(), domain ) ),
+    grid_( unstructured_grid( get() ) ) {}
+
 inline const StructuredGrid::grid_t* structured_grid( const Grid::Implementation* grid ) {
     return dynamic_cast<const StructuredGrid::grid_t*>( grid );
 }

src/atlas/grid/Grid.h

@@ -152,6 +152,7 @@ class UnstructuredGrid : public Grid {
     UnstructuredGrid( std::vector<PointXY>* );   // takes ownership
     UnstructuredGrid( std::vector<PointXY>&& );  // move constructor
     UnstructuredGrid( std::initializer_list<PointXY> );
+    UnstructuredGrid( const Grid&, const Domain& );  // Create a new unstructured grid!
 
     operator bool() const { return valid(); }
 

src/atlas/grid/detail/grid/Grid.cc

@@ -82,7 +82,7 @@ const Grid* Grid::create( const Grid& grid, const Domain& domain ) {
         return new Structured( g.name(), g.xspace(), g.yspace(), g.projection(), domain );
     }
     else {
-        NOTIMP;
+        return new Unstructured( grid, domain );
     }
 }
 

src/atlas/grid/detail/grid/Grid.h

@@ -48,6 +48,7 @@ class Grid : public eckit::Owned {
         virtual const IteratorXY& operator++()                   = 0;
         virtual bool operator==( const IteratorXY& other ) const = 0;
         virtual bool operator!=( const IteratorXY& other ) const = 0;
+        virtual ~IteratorXY() {}
     };
 
     class IteratorLonLat {
@@ -57,6 +58,7 @@ class Grid : public eckit::Owned {
         virtual const IteratorLonLat& operator++()                   = 0;
         virtual bool operator==( const IteratorLonLat& other ) const = 0;
         virtual bool operator!=( const IteratorLonLat& other ) const = 0;
+        virtual ~IteratorLonLat() {}
     };
 
 public:  // methods

src/atlas/grid/detail/grid/Structured.cc

@@ -78,11 +78,11 @@ Structured::Structured( const std::string& name, XSpace xspace, YSpace yspace, P
     }
     npts_ = size_t( std::accumulate( nx_.begin(), nx_.end(), 0 ) );
 
-    if ( not domain.empty() ) { crop( domain ); }
+    if ( domain ) { crop( domain ); }
 
     computeTruePeriodicity();
 
-    if ( domain.global() )
+    if ( domain && domain.global() )
         domain_ = Domain( Grid::Config( "type", "global" ) );
     else
         computeDomain();
@@ -99,7 +99,7 @@ void Structured::computeDomain() {
             domain_ = Domain( config );
         }
     }
-    else if ( domain_.empty() ) {
+    else if ( not domain_ ) {
         Grid::Config config;
         config.set( "type", "rectangular" );
         config.set( "xmin", xmin_[0] );
@@ -429,7 +429,7 @@ void Structured::crop( const Domain& dom ) {
 
 void Structured::computeTruePeriodicity() {
     if ( projection_.strictlyRegional() ) { periodic_x_ = false; }
-    else if ( domain_.global() ) {
+    else if ( domain_ && domain_.global() ) {
         periodic_x_ = true;
     }
     else {

src/atlas/grid/detail/grid/Unstructured.cc

@@ -11,13 +11,17 @@
 #include "atlas/grid/detail/grid/Unstructured.h"
 
 #include <limits>
+#include <memory>
 
 #include "eckit/memory/Builder.h"
+#include "eckit/types/FloatCompare.h"
 
 #include "atlas/array/ArrayView.h"
 #include "atlas/field/Field.h"
+#include "atlas/grid/Iterator.h"
 #include "atlas/mesh/Mesh.h"
 #include "atlas/mesh/Nodes.h"
+#include "atlas/option.h"
 #include "atlas/runtime/Log.h"
 #include "atlas/util/CoordinateEnums.h"
 
@@ -42,6 +46,64 @@ Unstructured::Unstructured( const Mesh& m ) : Grid(), points_( new std::vector<P
     }
 }
 
+
+namespace {
+class Normalise {
+public:
+    Normalise( const RectangularDomain& domain ) :
+        degrees_( domain.units() == "degrees" ),
+        xmin_( domain.xmin() ),
+        xmax_( domain.xmax() ),
+        eps_( 1e-11 ) {}
+
+    double operator()( double x ) const {
+        if ( degrees_ ) {
+            while ( eckit::types::is_strictly_greater<double>( xmin_, x, eps_ ) ) {
+                x += 360.;
+            }
+            while ( eckit::types::is_strictly_greater<double>( x, xmax_, eps_ ) ) {
+                x -= 360.;
+            }
+        }
+        return x;
+    }
+
+private:
+    const bool degrees_;
+    const double xmin_;
+    const double xmax_;
+    const double eps_;
+};
+}  // namespace
+
+
+Unstructured::Unstructured( const Grid& grid, Domain domain ) : Grid() {
+    domain_ = domain;
+    points_.reset( new std::vector<PointXY> );
+    points_->reserve( grid.size() );
+    if ( not domain_ ) { domain_ = Domain( option::type( "global" ) ); }
+    atlas::grid::IteratorXY it( grid.xy_begin() );
+    PointXY p;
+    if ( RectangularDomain( domain_ ) ) {
+        auto normalise = Normalise( RectangularDomain( domain_ ) );
+        while ( it.next( p ) ) {
+            p.x() = normalise( p.x() );
+            if ( domain_.contains( p ) ) { points_->emplace_back( p ); }
+        }
+    }
+    else if ( ZonalBandDomain( domain_ ) ) {
+        while ( it.next( p ) ) {
+            if ( domain_.contains( p ) ) { points_->emplace_back( p ); }
+        }
+    }
+    else {
+        while ( it.next( p ) ) {
+            points_->emplace_back( p );
+        }
+    }
+    points_->shrink_to_fit();
+}
+
 Unstructured::Unstructured( const util::Config& p ) : Grid() {
     util::Config config_domain;
     config_domain.set( "type", "global" );

src/atlas/grid/detail/grid/Unstructured.h

@@ -16,6 +16,7 @@
 #pragma once
 
 #include <cstddef>
+#include <functional>
 #include <memory>
 #include <vector>
 
@@ -154,11 +155,15 @@ class Unstructured : public Grid {
         size_t n_;
     };
 
+
 public:  // methods
     static std::string static_type() { return "unstructured"; }
     virtual std::string name() const;
     virtual std::string type() const { return static_type(); }
 
+    /// Constructor converting any Grid with domain to an unstructured grid
+    Unstructured( const Grid&, Domain );
+
     /// Constructor taking a list of parameters
     Unstructured( const Config& );
 

src/atlas/option/TransOptions.h

@@ -113,7 +113,7 @@ class nproma : public util::Config {
 
 class warning : public util::Config {
 public:
-    warning( int );
+    warning( int = 1 );
 };
 
 // ----------------------------------------------------------------------------