From 74e19607c3b8bd5f51d3dfcb34e35a2b1a217461 Mon Sep 17 00:00:00 2001 From: Tiago Quintino Date: Wed, 24 Apr 2019 18:38:17 +0100 Subject: [PATCH 01/11] Fix bamboo --- bamboo/CLANG-env.sh | 1 + bamboo/GCC-env.sh | 3 ++- bamboo/INTEL-env.sh | 2 ++ 3 files changed, 5 insertions(+), 1 deletion(-) diff --git a/bamboo/CLANG-env.sh b/bamboo/CLANG-env.sh index c7bfda5aa..f16fb80cb 100644 --- a/bamboo/CLANG-env.sh +++ b/bamboo/CLANG-env.sh @@ -18,6 +18,7 @@ module unload emos module unload fftw module unload libemos module unload metview +module unload netcdf4 module load cmake/3.10.2 diff --git a/bamboo/GCC-env.sh b/bamboo/GCC-env.sh index 3d58400c9..94b67cf7c 100644 --- a/bamboo/GCC-env.sh +++ b/bamboo/GCC-env.sh @@ -1,6 +1,6 @@ #!/bin/bash -[[ $(uname) == "Darwin" ]] && return +[[ $(uname) == "Darwin" ]] && return # no module environment on the Mac # initialise module environment if it is not if [[ ! $(command -v module > /dev/null 2>&1) ]]; then @@ -13,5 +13,6 @@ module unload emos module unload fftw module unload libemos module unload metview +module unload netcdf4 module load cmake/3.10.2 diff --git a/bamboo/INTEL-env.sh b/bamboo/INTEL-env.sh index 3b9fdb20e..3f77166ef 100644 --- a/bamboo/INTEL-env.sh +++ b/bamboo/INTEL-env.sh @@ -1,5 +1,7 @@ #!/bin/bash +[[ $(uname) == "Darwin" ]] && return # no module environment on the Mac + # initialise module environment if it is not if [[ ! $(command -v module > /dev/null 2>&1) ]]; then . /usr/local/apps/module/init/bash From 66b5bfa47ca1ce80797961ffbec0a415168e1965 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 7 May 2019 11:44:49 +0100 Subject: [PATCH 02/11] ATLAS-230 Fix PGI-19.4 compilation of GridTools array (internal compiler bugs) --- src/atlas/array/gridtools/GridToolsArray.cc | 3 + .../array/gridtools/GridToolsArrayHelpers.h | 36 ++++++++--- .../array/gridtools/GridToolsArrayView.h | 4 +- .../array/gridtools/GridToolsMakeView.cc | 64 ++++++++++--------- src/atlas/array/gridtools/GridToolsMakeView.h | 10 ++- src/tests/array/test_array.cc | 5 ++ 6 files changed, 79 insertions(+), 43 deletions(-) diff --git a/src/atlas/array/gridtools/GridToolsArray.cc b/src/atlas/array/gridtools/GridToolsArray.cc index 9a472a63a..538b1b194 100644 --- a/src/atlas/array/gridtools/GridToolsArray.cc +++ b/src/atlas/array/gridtools/GridToolsArray.cc @@ -64,6 +64,7 @@ class ArrayT_impl { template > void construct( UInts... dims ) { + static_assert( sizeof...(UInts) > 0, "1" ); auto gt_storage = create_gt_storage::type>( dims... ); using data_store_t = typename std::remove_pointer::type; array_.data_store_ = std::unique_ptr( new GridToolsDataStore( gt_storage ) ); @@ -464,6 +465,7 @@ template ArrayT::ArrayT( idx_t dim0 ) { ArrayT_impl( *this ).construct( dim0 ); } + template ArrayT::ArrayT( idx_t dim0, idx_t dim1 ) { ArrayT_impl( *this ).construct( dim0, dim1 ); @@ -569,3 +571,4 @@ template Array* Array::wrap( long unsigned*, const ArraySpec& ); } // namespace array } // namespace atlas + diff --git a/src/atlas/array/gridtools/GridToolsArrayHelpers.h b/src/atlas/array/gridtools/GridToolsArrayHelpers.h index cff711378..4f4e064e3 100644 --- a/src/atlas/array/gridtools/GridToolsArrayHelpers.h +++ b/src/atlas/array/gridtools/GridToolsArrayHelpers.h @@ -148,10 +148,20 @@ struct get_pack_size { using type = ::gridtools::static_uint; }; +template +struct gt_storage_t { + using type = gridtools::storage_traits::data_store_t< + Value, + gridtools::storage_traits::custom_layout_storage_info_t< + 0, + LayoutMap, + ::gridtools::zero_halo + > + >; +}; + template -static gridtools::storage_traits::data_store_t< - Value, gridtools::storage_traits::custom_layout_storage_info_t< - 0, LayoutMap, ::gridtools::zero_halo::type::value>>>* +typename gt_storage_t::type::value>::type* create_gt_storage( UInts... dims ) { static_assert( ( sizeof...( dims ) > 0 ), "Error: can not create storages without any dimension" ); @@ -161,16 +171,24 @@ create_gt_storage( UInts... dims ) { typedef gridtools::storage_traits::data_store_t data_store_t; data_store_t* ds; - if (::gridtools::accumulate(::gridtools::multiplies(), dims... ) == 0 ) { ds = new data_store_t(); } - else { - storage_info_ty si( dims... ); - ds = new data_store_t( si ); - } + if (::gridtools::accumulate(::gridtools::multiplies(), dims... ) == 0 ) { ds = new data_store_t(); } + else { + storage_info_ty si( dims... ); + ds = new data_store_t( si ); + } return ds; } +template +struct gt_wrap_storage_t { + using type = gridtools::storage_traits::data_store_t< + Value, + gridtools::storage_traits::storage_info_t<0, Rank> + >; +}; + template -static gridtools::storage_traits::data_store_t>* +static typename gt_wrap_storage_t::type* wrap_gt_storage( Value* data, std::array&& shape, std::array&& strides ) { static_assert( ( Rank > 0 ), "Error: can not create storages without any dimension" ); typedef gridtools::storage_traits::storage_info_t<0, Rank, ::gridtools::zero_halo> storage_info_ty; diff --git a/src/atlas/array/gridtools/GridToolsArrayView.h b/src/atlas/array/gridtools/GridToolsArrayView.h index 65f38d73c..ad57269f3 100644 --- a/src/atlas/array/gridtools/GridToolsArrayView.h +++ b/src/atlas/array/gridtools/GridToolsArrayView.h @@ -14,10 +14,10 @@ #include #include +#include "atlas/array/Array.h" #include "atlas/array/ArrayUtil.h" #include "atlas/array/ArrayViewDefs.h" #include "atlas/array/LocalView.h" -#include "atlas/array/gridtools/GridToolsMakeView.h" #include "atlas/array/gridtools/GridToolsTraits.h" #include "atlas/library/config.h" @@ -26,6 +26,7 @@ namespace atlas { namespace array { + template class ArrayView { public: @@ -57,6 +58,7 @@ class ArrayView { ATLAS_HOST_DEVICE ArrayView( const ArrayView& other ); ArrayView( data_view_t data_view, const Array& array ); + value_type* data() { return gt_data_view_.data(); } value_type const* data() const { return gt_data_view_.data(); } diff --git a/src/atlas/array/gridtools/GridToolsMakeView.cc b/src/atlas/array/gridtools/GridToolsMakeView.cc index 3e9a6fe48..44dc17fb8 100644 --- a/src/atlas/array/gridtools/GridToolsMakeView.cc +++ b/src/atlas/array/gridtools/GridToolsMakeView.cc @@ -13,15 +13,14 @@ #include #include -#include "atlas/array.h" +#include "atlas/array/Array.h" +#include "atlas/array/ArrayViewDefs.h" #include "atlas/array/ArrayView.h" #include "atlas/array/IndexView.h" +#include "atlas/array/gridtools/GridToolsTraits.h" #include "atlas/runtime/Exception.h" - #include "atlas/library/config.h" -#if ATLAS_HAVE_GRIDTOOLS_STORAGE -#include "atlas/array/gridtools/GridToolsTraits.h" -#endif + //------------------------------------------------------------------------------ namespace atlas { @@ -45,20 +44,21 @@ static void check_metadata( const Array& array ) { namespace gridtools { + template -data_view_tt make_gt_host_view( const Array& array ) { +typename gt_view::type +make_gt_host_view( const Array& array ) { using storage_info_ty = storage_traits::storage_info_t<0, Rank>; using data_store_t = storage_traits::data_store_t; data_store_t* ds = reinterpret_cast( const_cast( array.storage() ) ); - return ::gridtools::make_host_view( *ds ); } template -data_view_tt make_gt_device_view( const Array& array ) { - typedef storage_traits::storage_info_t<0, Rank> storage_info_ty; - typedef storage_traits::data_store_t data_store_t; +typename gt_view::type make_gt_device_view( const Array& array ) { + using storage_info_ty = storage_traits::storage_info_t<0, Rank>; + using data_store_t = storage_traits::data_store_t; data_store_t* ds = reinterpret_cast( const_cast( array.storage() ) ); #if ATLAS_GRIDTOOLS_STORAGE_BACKEND_CUDA @@ -67,6 +67,7 @@ data_view_tt make_gt_device_view( co return ::gridtools::make_host_view( *ds ); #endif } + } // namespace gridtools template @@ -114,6 +115,7 @@ ArrayView make_view( const Array& array ) { // Explicit instantiation namespace atlas { namespace array { + #define EXPLICIT_TEMPLATE_INSTANTIATION( RANK ) \ template ArrayView make_view( const Array& ); \ template ArrayView make_view( const Array& ); \ @@ -169,48 +171,47 @@ namespace array { \ template IndexView make_host_indexview( const Array& ); \ template IndexView make_host_indexview( const Array& ); \ - \ namespace gridtools { \ - template data_view_tt \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ \ - template data_view_tt \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template data_view_tt \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ } @@ -226,5 +227,6 @@ EXPLICIT_TEMPLATE_INSTANTIATION( 8 ) EXPLICIT_TEMPLATE_INSTANTIATION( 9 ) #undef EXPLICIT_TEMPLATE_INSTANTIATION + } // namespace array } // namespace atlas diff --git a/src/atlas/array/gridtools/GridToolsMakeView.h b/src/atlas/array/gridtools/GridToolsMakeView.h index 134ce5a42..735563a7f 100644 --- a/src/atlas/array/gridtools/GridToolsMakeView.h +++ b/src/atlas/array/gridtools/GridToolsMakeView.h @@ -16,11 +16,17 @@ namespace atlas { namespace array { namespace gridtools { +template +struct gt_view { + using type = ::gridtools::data_view< + gridtools::storage_traits::data_store_t>, gridtools::get_access_mode(AccessMode) >; +}; + template -data_view_tt make_gt_host_view( const Array& array ); +typename gt_view::type make_gt_host_view( const Array& array ); template -data_view_tt make_gt_device_view( const Array& array ); +typename gt_view::type make_gt_device_view( const Array& array ); } // namespace gridtools } // namespace array diff --git a/src/tests/array/test_array.cc b/src/tests/array/test_array.cc index 5c6e5b4ab..4bf6394f1 100644 --- a/src/tests/array/test_array.cc +++ b/src/tests/array/test_array.cc @@ -15,6 +15,11 @@ #include "atlas/library/config.h" #include "tests/AtlasTestEnvironment.h" +#if ATLAS_HAVE_GRIDTOOLS_STORAGE +#include "atlas/array/gridtools/GridToolsMakeView.h" +#endif + + #ifdef ATLAS_HAVE_GRIDTOOLS_STORAGE #if ATLAS_GRIDTOOLS_STORAGE_BACKEND_CUDA #define PADDED 1 From 97f8d093cdc4aef0a5f6984865e60baf2d757d23 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 13 May 2019 18:32:08 +0100 Subject: [PATCH 03/11] Fix static linking for LAEA projection --- src/atlas/projection/detail/ProjectionFactory.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/atlas/projection/detail/ProjectionFactory.cc b/src/atlas/projection/detail/ProjectionFactory.cc index 4522fa71b..726faa8bd 100644 --- a/src/atlas/projection/detail/ProjectionFactory.cc +++ b/src/atlas/projection/detail/ProjectionFactory.cc @@ -14,6 +14,7 @@ #include "atlas/projection/detail/ProjectionFactory.h" #include "atlas/util/Config.h" +#include "atlas/projection/detail/LambertAzimuthalEqualAreaProjection.h" #include "atlas/projection/detail/LambertProjection.h" #include "atlas/projection/detail/LonLatProjection.h" #include "atlas/projection/detail/MercatorProjection.h" @@ -34,6 +35,7 @@ void force_link() { ProjectionBuilder(); ProjectionBuilder(); ProjectionBuilder(); + ProjectionBuilder(); } } link; } From 9e9e32eab608dc9f2dab5e13423453f3c32a322e Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 13 May 2019 18:32:21 +0100 Subject: [PATCH 04/11] Fix multiple defined symbol --- src/atlas/mesh/actions/BuildNode2CellConnectivity.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atlas/mesh/actions/BuildNode2CellConnectivity.h b/src/atlas/mesh/actions/BuildNode2CellConnectivity.h index 58f0d3ac2..636a145ea 100644 --- a/src/atlas/mesh/actions/BuildNode2CellConnectivity.h +++ b/src/atlas/mesh/actions/BuildNode2CellConnectivity.h @@ -27,7 +27,7 @@ class BuildNode2CellConnectivity { Mesh& mesh_; }; -void build_node_to_cell_connectivity( Mesh& mesh ) { +inline void build_node_to_cell_connectivity( Mesh& mesh ) { BuildNode2CellConnectivity{mesh}(); } From 34f16cf8e8f1b32be30226cc8ddd3dc10e66a3ba Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 31 May 2019 16:15:03 +0100 Subject: [PATCH 05/11] ATLAS-232 Add knowledge of 'west' to ZonalBandDomain and GlobalDomain --- src/atlas/domain/Domain.cc | 9 +++++---- src/atlas/domain/Domain.h | 2 +- src/atlas/domain/detail/GlobalDomain.cc | 15 +++++++++++++-- src/atlas/domain/detail/GlobalDomain.h | 2 +- src/atlas/domain/detail/RectangularDomain.cc | 2 +- src/atlas/domain/detail/ZonalBandDomain.cc | 13 +++++++++++-- src/atlas/domain/detail/ZonalBandDomain.h | 2 +- 7 files changed, 33 insertions(+), 12 deletions(-) diff --git a/src/atlas/domain/Domain.cc b/src/atlas/domain/Domain.cc index 3149f625a..9aad26338 100644 --- a/src/atlas/domain/Domain.cc +++ b/src/atlas/domain/Domain.cc @@ -26,7 +26,8 @@ RectangularDomain::RectangularDomain( const Interval& x, const Interval& y, cons Domain( ( RD::is_global( x, y, units ) ) ? new atlas::domain::GlobalDomain( x[0] ) : ( RD::is_zonal_band( x, units ) ? new atlas::domain::ZonalBandDomain( y, x[0] ) - : new atlas::domain::RectangularDomain( x, y, units ) ) ) {} + : new atlas::domain::RectangularDomain( x, y, units ) ) ) { +} RectangularDomain::RectangularDomain( const Domain& domain ) : Domain( domain ), @@ -56,9 +57,9 @@ double RectangularDomain::ymax() const { return domain_->ymax(); } -ZonalBandDomain::ZonalBandDomain( const Interval& y ) : - RectangularDomain( ( ZD::is_global( y ) ) ? new atlas::domain::GlobalDomain() - : new atlas::domain::ZonalBandDomain( y ) ) {} +ZonalBandDomain::ZonalBandDomain( const Interval& y, const double& west ) : + RectangularDomain( ( ZD::is_global( y ) ) ? new atlas::domain::GlobalDomain(west) + : new atlas::domain::ZonalBandDomain( y, west ) ) {} ZonalBandDomain::ZonalBandDomain( const Domain& domain ) : RectangularDomain( domain ), diff --git a/src/atlas/domain/Domain.h b/src/atlas/domain/Domain.h index d049d436e..60abbbab4 100644 --- a/src/atlas/domain/Domain.h +++ b/src/atlas/domain/Domain.h @@ -131,7 +131,7 @@ class ZonalBandDomain : public RectangularDomain { public: using RectangularDomain::RectangularDomain; - ZonalBandDomain( const Interval& y ); + ZonalBandDomain( const Interval& y, const double& west = 0. ); ZonalBandDomain( const Domain& ); operator bool() { return domain_; } diff --git a/src/atlas/domain/detail/GlobalDomain.cc b/src/atlas/domain/detail/GlobalDomain.cc index 0ca20f472..57a7137e9 100644 --- a/src/atlas/domain/detail/GlobalDomain.cc +++ b/src/atlas/domain/detail/GlobalDomain.cc @@ -22,17 +22,28 @@ namespace { constexpr std::array yrange() { return {-90., 90.}; } + +static double get_west( const eckit::Parametrisation& params ) { + double west = 0.; + params.get("west",west); + return west; +} + } // namespace GlobalDomain::GlobalDomain( const double west ) : ZonalBandDomain( yrange(), west ) {} GlobalDomain::GlobalDomain() : ZonalBandDomain( yrange() ) {} -GlobalDomain::GlobalDomain( const eckit::Parametrisation& ) : GlobalDomain() {} +GlobalDomain::GlobalDomain( const eckit::Parametrisation& params ) : + GlobalDomain( get_west( params ) ) {} GlobalDomain::Spec GlobalDomain::spec() const { Spec domain_spec; domain_spec.set( "type", type() ); + if( xmin() != 0. ) { + domain_spec.set( "west", xmin() ); + } return domain_spec; } @@ -41,7 +52,7 @@ void GlobalDomain::hash( eckit::Hash& h ) const { } void GlobalDomain::print( std::ostream& os ) const { - os << "GlobalDomain"; + os << "GlobalDomain[west=" << xmin() << "]" ; } namespace { diff --git a/src/atlas/domain/detail/GlobalDomain.h b/src/atlas/domain/detail/GlobalDomain.h index a11ad64cf..dffa6afec 100644 --- a/src/atlas/domain/detail/GlobalDomain.h +++ b/src/atlas/domain/detail/GlobalDomain.h @@ -18,6 +18,7 @@ namespace domain { class GlobalDomain : public ZonalBandDomain { public: GlobalDomain(); + GlobalDomain( const double west ); GlobalDomain( const eckit::Parametrisation& p ); static std::string static_type() { return "global"; } @@ -44,7 +45,6 @@ class GlobalDomain : public ZonalBandDomain { private: friend class ::atlas::RectangularDomain; - GlobalDomain( const double west ); }; } // namespace domain diff --git a/src/atlas/domain/detail/RectangularDomain.cc b/src/atlas/domain/detail/RectangularDomain.cc index 66e5e52a6..306d5a546 100644 --- a/src/atlas/domain/detail/RectangularDomain.cc +++ b/src/atlas/domain/detail/RectangularDomain.cc @@ -54,7 +54,7 @@ bool RectangularDomain::is_global( const Interval& x, const Interval& y, const s if ( units != "degrees" ) return false; const double eps = 1.e-12; - return std::abs( ( x[1] - x[0] ) - 360. ) < eps && std::abs( ( y[1] - y[0] ) - 180. ) < eps; + return std::abs( ( x[1] - x[0] ) - 360. ) < eps && std::abs( std::abs( y[1] - y[0] ) - 180. ) < eps; } bool RectangularDomain::is_zonal_band( const Interval& x, const std::string& units ) { diff --git a/src/atlas/domain/detail/ZonalBandDomain.cc b/src/atlas/domain/detail/ZonalBandDomain.cc index 410c4cc5d..507ae1e74 100644 --- a/src/atlas/domain/detail/ZonalBandDomain.cc +++ b/src/atlas/domain/detail/ZonalBandDomain.cc @@ -34,6 +34,12 @@ static std::array get_interval_y( const eckit::Parametrisation& param return {ymin, ymax}; } +static double get_west( const eckit::Parametrisation& params ) { + double west = 0.; + params.get("west",west); + return west; +} + /* constexpr std::array interval_x() { return {0., 360.}; @@ -49,7 +55,7 @@ bool ZonalBandDomain::is_global( const Interval& y ) { } ZonalBandDomain::ZonalBandDomain( const eckit::Parametrisation& params ) : - ZonalBandDomain( get_interval_y( params ) ) {} + ZonalBandDomain( get_interval_y( params ), get_west( params ) ) {} ZonalBandDomain::ZonalBandDomain( const Interval& interval_y ) : ZonalBandDomain( interval_y, /*west*/ 0. ) {} @@ -69,6 +75,9 @@ ZonalBandDomain::Spec ZonalBandDomain::spec() const { domain_spec.set( "type", type() ); domain_spec.set( "ymin", ymin() ); domain_spec.set( "ymax", ymax() ); + if( xmin() != 0. ) { + domain_spec.set( "west", xmin() ); + } return domain_spec; } @@ -78,7 +87,7 @@ void ZonalBandDomain::hash( eckit::Hash& h ) const { void ZonalBandDomain::print( std::ostream& os ) const { os << "ZonalBandDomain[" - << "ymin=" << ymin() << ",ymax=" << ymax() << "]"; + << "ymin=" << ymin() << ",ymax=" << ymax() << ",west=" << xmin() << "]"; } bool ZonalBandDomain::containsNorthPole() const { diff --git a/src/atlas/domain/detail/ZonalBandDomain.h b/src/atlas/domain/detail/ZonalBandDomain.h index e277112dd..4e8494edc 100644 --- a/src/atlas/domain/detail/ZonalBandDomain.h +++ b/src/atlas/domain/detail/ZonalBandDomain.h @@ -34,6 +34,7 @@ class ZonalBandDomain : public atlas::domain::RectangularDomain { // constructor ZonalBandDomain( const eckit::Parametrisation& ); ZonalBandDomain( const Interval& ); + ZonalBandDomain( const Interval&, const double west ); static std::string static_type() { return "zonal_band"; } virtual std::string type() const override { return static_type(); } @@ -60,7 +61,6 @@ class ZonalBandDomain : public atlas::domain::RectangularDomain { protected: friend class ::atlas::RectangularDomain; - ZonalBandDomain( const Interval&, const double west ); private: bool global_; From 05defdfa374fcd0240c08c33caaa5fb4da66348a Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 31 May 2019 16:23:56 +0100 Subject: [PATCH 06/11] ATLAS-232 Global grids cropped with zonal or global domains with non-zero west are now respected --- src/atlas/grid/detail/grid/Structured.cc | 136 +++++++++++------------ src/atlas/grid/detail/grid/Structured.h | 18 +++ 2 files changed, 81 insertions(+), 73 deletions(-) diff --git a/src/atlas/grid/detail/grid/Structured.cc b/src/atlas/grid/detail/grid/Structured.cc index 4250d765e..304107847 100644 --- a/src/atlas/grid/detail/grid/Structured.cc +++ b/src/atlas/grid/detail/grid/Structured.cc @@ -51,7 +51,7 @@ Structured::Structured( const std::string& name, XSpace xspace, YSpace yspace, P name_( name ), xspace_( xspace ), yspace_( yspace ) { - // Copry members + // Copy members if ( projection ) projection_ = projection; else @@ -87,32 +87,15 @@ Structured::Structured( const std::string& name, XSpace xspace, YSpace yspace, P computeTruePeriodicity(); - if ( domain && domain.global() ) - domain_ = Domain( Grid::Config( "type", "global" ) ); - else - computeDomain(); + if ( not domain ) { computeDomain(); } } void Structured::computeDomain() { if ( periodic() ) { - if ( yspace().max() - yspace().min() == 180. ) { domain_ = Domain( Grid::Config( "type", "global" ) ); } - else { - Grid::Config config; - config.set( "type", "zonal_band" ); - config.set( "ymin", yspace().min() ); - config.set( "ymax", yspace().max() ); - domain_ = Domain( config ); - } + domain_ = ZonalBandDomain( {yspace().min(),yspace().max()}, xspace().min() ); } - else if ( not domain_ ) { - Grid::Config config; - config.set( "type", "rectangular" ); - config.set( "xmin", xmin_[0] ); - config.set( "xmax", xmax_[0] ); - config.set( "ymin", yspace().min() ); - config.set( "ymax", yspace().max() ); - config.set( "units", projection_.units() ); - domain_ = Domain( config ); + else { + domain_ = RectangularDomain( {xspace().min(),xspace().max()}, {yspace().min(),yspace().max()}, projection_.units() ); } } @@ -133,6 +116,8 @@ Structured::XSpace::XSpace( const std::array& interval, std::initiali Structured::XSpace::XSpace( const Spacing& spacing ) : impl_( new Implementation( spacing ) ) {} +Structured::XSpace::XSpace( const std::vector& spacings ) : impl_( new Implementation( spacings ) ) {} + Structured::XSpace::XSpace( const Config& config ) : impl_( new Implementation( config ) ) {} Structured::XSpace::XSpace( const std::vector& config ) : impl_( new Implementation( config ) ) {} @@ -168,6 +153,8 @@ Structured::XSpace::Implementation::Implementation( const Config& config ) { nxmin_ = std::numeric_limits::max(); nxmax_ = 0; + min_ = std::numeric_limits::max(); + max_ = -std::numeric_limits::max(); for ( idx_t j = 0; j < ny; ++j ) { if ( not v_N.empty() ) config_xspace.set( "N", v_N[j] ); @@ -181,6 +168,8 @@ Structured::XSpace::Implementation::Implementation( const Config& config ) { dx_.push_back( xspace.step ); nxmin_ = std::min( nxmin_, nx_[j] ); nxmax_ = std::max( nxmax_, nx_[j] ); + min_ = std::min( min_, xspace.start ); + max_ = std::max( max_, xspace.end ); } } @@ -189,6 +178,8 @@ Structured::XSpace::Implementation::Implementation( const std::vector& c nxmin_ = std::numeric_limits::max(); nxmax_ = 0; + min_ = std::numeric_limits::max(); + max_ = -std::numeric_limits::max(); std::string xspace_type; for ( idx_t j = 0; j < ny(); ++j ) { @@ -201,6 +192,8 @@ Structured::XSpace::Implementation::Implementation( const std::vector& c dx_.push_back( xspace.step ); nxmin_ = std::min( nxmin_, nx_[j] ); nxmax_ = std::max( nxmax_, nx_[j] ); + min_ = std::min( min_, xspace.start ); + max_ = std::max( max_, xspace.end ); } } @@ -222,11 +215,15 @@ Structured::XSpace::Implementation::Implementation( const std::array& dx_( ny_ ) { nxmin_ = std::numeric_limits::max(); nxmax_ = 0; + min_ = std::numeric_limits::max(); + max_ = -std::numeric_limits::max(); double length = interval[1] - interval[0]; for ( idx_t j = 0; j < ny_; ++j ) { nxmin_ = std::min( nxmin_, nx_[j] ); nxmax_ = std::max( nxmax_, nx_[j] ); dx_[j] = endpoint ? length / double( nx_[j] - 1 ) : length / double( nx_[j] ); + min_ = std::min( min_, xmin_[j] ); + max_ = std::max( max_, xmax_[j] ); } } template Structured::XSpace::Implementation::Implementation( const std::array& interval, @@ -249,8 +246,35 @@ Structured::XSpace::Implementation::Implementation( const Spacing& spacing ) : dx_[0] = linspace.step(); nxmax_ = nx_[0]; nxmin_ = nx_[0]; + min_ = spacing.min(); + max_ = spacing.max(); +} + +Structured::XSpace::Implementation::Implementation( const std::vector& spacings ) : + ny_( spacings.size() ), + nx_( ny_ ), + xmin_( ny_ ), + xmax_( ny_ ), + dx_( ny_ ) { + nxmax_ = 0; + nxmin_ = std::numeric_limits::max(); + min_ = std::numeric_limits::max(); + max_ = -std::numeric_limits::max(); + for ( idx_t j = 0; j < ny_; ++j ) { + const spacing::LinearSpacing& linspace = dynamic_cast( *spacings[j].get() ); + + nx_[j] = linspace.size(); + xmin_[j] = linspace.min(); + xmax_[j] = linspace.max(); + dx_[j] = linspace.step(); + nxmin_ = std::min( nxmin_, nx_[j] ); + nxmax_ = std::max( nxmax_, nx_[j] ); + min_ = std::min( min_, xmin_[j] ); + max_ = std::max( max_, xmax_[j] ); + } } + std::string Structured::XSpace::Implementation::type() const { return "linear"; } @@ -302,7 +326,7 @@ class Normalise { public: Normalise( const RectangularDomain& domain ) : degrees_( domain.units() == "degrees" ), - normalise_( domain.xmin(), domain.xmax() ) {} + normalise_( domain.xmin() ) {} double operator()( double x ) const { if ( degrees_ ) { x = normalise_( x ); } @@ -316,59 +340,11 @@ class Normalise { } // namespace void Structured::crop( const Domain& dom ) { - if ( dom.global() ) return; - ATLAS_ASSERT( dom.units() == projection().units() ); - auto zonal_domain = ZonalBandDomain( dom ); - auto rect_domain = RectangularDomain( dom ); + auto rect_domain = RectangularDomain( dom ); - if ( zonal_domain ) { - const double cropped_ymin = zonal_domain.ymin(); - const double cropped_ymax = zonal_domain.ymax(); - - idx_t jmin = ny(); - idx_t jmax = 0; - for ( idx_t j = 0; j < ny(); ++j ) { - if ( zonal_domain.contains_y( y( j ) ) ) { - jmin = std::min( j, jmin ); - jmax = std::max( j, jmax ); - } - } - idx_t cropped_ny = jmax - jmin + 1; - std::vector cropped_y( y_.begin() + jmin, y_.begin() + jmin + cropped_ny ); - std::vector cropped_xmin( xmin_.begin() + jmin, xmin_.begin() + jmin + cropped_ny ); - std::vector cropped_xmax( xmax_.begin() + jmin, xmax_.begin() + jmin + cropped_ny ); - std::vector cropped_dx( dx_.begin() + jmin, dx_.begin() + jmin + cropped_ny ); - std::vector cropped_nx( nx_.begin() + jmin, nx_.begin() + jmin + cropped_ny ); - ATLAS_ASSERT( idx_t( cropped_nx.size() ) == cropped_ny ); - - idx_t cropped_nxmin, cropped_nxmax; - cropped_nxmin = cropped_nxmax = cropped_nx.front(); - for ( idx_t j = 1; j < cropped_ny; ++j ) { - cropped_nxmin = std::min( cropped_nx[j], cropped_nxmin ); - cropped_nxmax = std::max( cropped_nx[j], cropped_nxmax ); - } - idx_t cropped_npts = std::accumulate( cropped_nx.begin(), cropped_nx.end(), idx_t{0} ); - - Spacing cropped_yspace( - new spacing::CustomSpacing( cropped_ny, cropped_y.data(), {cropped_ymin, cropped_ymax} ) ); - - // Modify grid - { - domain_ = dom; - yspace_ = cropped_yspace; - xmin_ = cropped_xmin; - xmax_ = cropped_xmax; - dx_ = cropped_dx; - nx_ = cropped_nx; - nxmin_ = cropped_nxmin; - nxmax_ = cropped_nxmax; - npts_ = cropped_npts; - y_ = cropped_y; - } - } - else if ( rect_domain ) { + if ( rect_domain ) { const double cropped_ymin = rect_domain.ymin(); const double cropped_ymax = rect_domain.ymax(); @@ -420,10 +396,24 @@ void Structured::crop( const Domain& dom ) { Spacing cropped_yspace( new spacing::CustomSpacing( cropped_ny, cropped_y.data(), {cropped_ymin, cropped_ymax} ) ); + std::vector cropped_xspacings( cropped_ny ); + for ( idx_t j = 0; j < cropped_ny; ++j ) { + cropped_xspacings[j] = new spacing::LinearSpacing( cropped_xmin[j], cropped_xmax[j], cropped_nx[j], true ); + } + XSpace cropped_xspace( cropped_xspacings ); + + for ( idx_t j = 0; j < cropped_ny; ++j ) { + ATLAS_ASSERT( eckit::types::is_approximately_equal( cropped_xspace.xmin()[j], cropped_xmin[j] ) ); + ATLAS_ASSERT( eckit::types::is_approximately_equal( cropped_xspace.xmax()[j], cropped_xmax[j] ) ); + ATLAS_ASSERT( cropped_xspace.nxmin() == cropped_nxmin ); + ATLAS_ASSERT( cropped_xspace.nxmax() == cropped_nxmax ); + } + // Modify grid { domain_ = dom; yspace_ = cropped_yspace; + xspace_ = cropped_xspace; xmin_ = cropped_xmin; xmax_ = cropped_xmax; dx_ = cropped_dx; diff --git a/src/atlas/grid/detail/grid/Structured.h b/src/atlas/grid/detail/grid/Structured.h index 23c0699ba..95e126a8a 100644 --- a/src/atlas/grid/detail/grid/Structured.h +++ b/src/atlas/grid/detail/grid/Structured.h @@ -195,6 +195,8 @@ class Structured : public Grid { Implementation( const Spacing& ); + Implementation( const std::vector& ); + Implementation( const Config& ); Implementation( const std::vector& ); @@ -219,6 +221,12 @@ class Structured : public Grid { /// Value of longitude increment const std::vector& dx() const { return dx_; } + /// Value of minimum x over entire grid + double min() const { return min_; } + + /// Value of maximum x over entire grid + double max() const { return max_; } + Spec spec() const; std::string type() const; @@ -234,6 +242,8 @@ class Structured : public Grid { std::vector xmin_; std::vector xmax_; std::vector dx_; + double min_; + double max_; }; public: @@ -243,6 +253,8 @@ class Structured : public Grid { XSpace( const Spacing& ); + XSpace( const std::vector& ); + // Constructor NVector can be either std::vector or std::vector or initializer list template XSpace( const std::array& interval, const NVector& N, bool endpoint = true ); @@ -273,6 +285,12 @@ class Structured : public Grid { /// Value of longitude increment const std::vector& dx() const { return impl_->dx(); } + /// Value of minimum x over entire grid + double min() const { return impl_->min(); } + + /// Value of maximum x over entire grid + double max() const { return impl_->max(); } + Spec spec() const; std::string type() const { return impl_->type(); } From 1b80652fc28c9aae8e19d219d9f6e6081a8cbf64 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 31 May 2019 17:46:26 +0100 Subject: [PATCH 07/11] Nicer error message for unsupported TransLocal grids --- src/atlas/trans/local/TransLocal.cc | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/atlas/trans/local/TransLocal.cc b/src/atlas/trans/local/TransLocal.cc index 4f4802680..c8f214a32 100644 --- a/src/atlas/trans/local/TransLocal.cc +++ b/src/atlas/trans/local/TransLocal.cc @@ -324,8 +324,8 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma gridGlobal_ = grid; if ( not gridGlobal_.domain().global() ) { + // The grid is not a nest of a global grid if ( RegularGrid( grid_ ) ) { - // non-nested regular grid no_nest = true; no_symmetry_ = true; useFFT_ = false; @@ -336,7 +336,14 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma useGlobalLeg = false; } else { - ATLAS_NOTIMPLEMENTED; + Log::info() << "Global grid: " << gridGlobal_.spec() << std::endl; + Log::info() << "Actual grid: " << grid_.spec() << std::endl; + std::stringstream log; + log << "Transforms to non-regular regional grids is not supported, unless it defined as a cropping of a global grid." << std::endl; + log << " Input grid: " << grid.spec() << std::endl; + log << " Applied domain: " << domain << std::endl; + log << " Regional grid after Domain applied (the output grid): " << grid_.spec() << std::endl; + throw_NotImplemented( log.str(), Here() ); // non-nested reduced grids are not supported } } From f5dd1bed3ababbb072a97fba43b907d95d2ddbae Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 3 Jun 2019 17:37:04 +0100 Subject: [PATCH 08/11] ATLAS-232 Fixes to Structured grid creation --- src/atlas/grid/detail/grid/Structured.cc | 203 ++++++++++++----------- 1 file changed, 107 insertions(+), 96 deletions(-) diff --git a/src/atlas/grid/detail/grid/Structured.cc b/src/atlas/grid/detail/grid/Structured.cc index 304107847..e2260e27e 100644 --- a/src/atlas/grid/detail/grid/Structured.cc +++ b/src/atlas/grid/detail/grid/Structured.cc @@ -83,11 +83,10 @@ Structured::Structured( const std::string& name, XSpace xspace, YSpace yspace, P } npts_ = std::accumulate( nx_.begin(), nx_.end(), idx_t{0} ); - if ( domain ) { crop( domain ); } + crop( domain ); computeTruePeriodicity(); - if ( not domain ) { computeDomain(); } } void Structured::computeDomain() { @@ -340,123 +339,135 @@ class Normalise { } // namespace void Structured::crop( const Domain& dom ) { - ATLAS_ASSERT( dom.units() == projection().units() ); + if ( dom ) { + ATLAS_ASSERT( dom.units() == projection().units() ); - auto rect_domain = RectangularDomain( dom ); - if ( rect_domain ) { - const double cropped_ymin = rect_domain.ymin(); - const double cropped_ymax = rect_domain.ymax(); + auto rect_domain = RectangularDomain( dom ); - // Cropping in Y - idx_t jmin = ny(); - idx_t jmax = 0; - for ( idx_t j = 0; j < ny(); ++j ) { - if ( rect_domain.contains_y( y( j ) ) ) { - jmin = std::min( j, jmin ); - jmax = std::max( j, jmax ); + if ( rect_domain ) { + const double cropped_ymin = rect_domain.ymin(); + const double cropped_ymax = rect_domain.ymax(); + + // Cropping in Y + idx_t jmin = ny(); + idx_t jmax = 0; + for ( idx_t j = 0; j < ny(); ++j ) { + if ( rect_domain.contains_y( y( j ) ) ) { + jmin = std::min( j, jmin ); + jmax = std::max( j, jmax ); + } } - } - ATLAS_ASSERT( jmax >= jmin ); - - idx_t cropped_ny = jmax - jmin + 1; - std::vector cropped_y( y_.begin() + jmin, y_.begin() + jmin + cropped_ny ); - std::vector cropped_dx( dx_.begin() + jmin, dx_.begin() + jmin + cropped_ny ); - - std::vector cropped_xmin( cropped_ny, std::numeric_limits::max() ); - std::vector cropped_xmax( cropped_ny, -std::numeric_limits::max() ); - std::vector cropped_nx( cropped_ny ); - - // Cropping in X - Normalise normalise( rect_domain ); - for ( idx_t j = jmin, jcropped = 0; j <= jmax; ++j, ++jcropped ) { - idx_t n = 0; - for ( idx_t i = 0; i < nx( j ); ++i ) { - const double _x = normalise( x( i, j ) ); - if ( rect_domain.contains_x( _x ) ) { - cropped_xmin[jcropped] = std::min( cropped_xmin[jcropped], _x ); - cropped_xmax[jcropped] = std::max( cropped_xmax[jcropped], _x ); - ++n; + ATLAS_ASSERT( jmax >= jmin ); + + idx_t cropped_ny = jmax - jmin + 1; + std::vector cropped_y( y_.begin() + jmin, y_.begin() + jmin + cropped_ny ); + std::vector cropped_dx( dx_.begin() + jmin, dx_.begin() + jmin + cropped_ny ); + + std::vector cropped_xmin( cropped_ny, std::numeric_limits::max() ); + std::vector cropped_xmax( cropped_ny, -std::numeric_limits::max() ); + std::vector cropped_nx( cropped_ny ); + + // Cropping in X + Normalise normalise( rect_domain ); + for ( idx_t j = jmin, jcropped = 0; j <= jmax; ++j, ++jcropped ) { + idx_t n = 0; + for ( idx_t i = 0; i < nx( j ); ++i ) { + const double _x = normalise( x( i, j ) ); + if ( rect_domain.contains_x( _x ) ) { + cropped_xmin[jcropped] = std::min( cropped_xmin[jcropped], _x ); + cropped_xmax[jcropped] = std::max( cropped_xmax[jcropped], _x ); + ++n; + } + } + cropped_nx[jcropped] = n; + } + bool endpoint = true; + if ( ZonalBandDomain( rect_domain ) ) { + for ( idx_t j = 0; j < cropped_ny; ++j ) { + if ( eckit::types::is_approximately_equal( cropped_xmax[j] + cropped_dx[j], + cropped_xmin[j] + 360., 1.e-10 ) ) { + cropped_xmax[j] = cropped_xmin[j] + 360.; + endpoint = false; + } } } - cropped_nx[jcropped] = n; - } - // Complete structures + // Complete structures - idx_t cropped_nxmin, cropped_nxmax; - cropped_nxmin = cropped_nxmax = cropped_nx.front(); + idx_t cropped_nxmin, cropped_nxmax; + cropped_nxmin = cropped_nxmax = cropped_nx.front(); - for ( idx_t j = 1; j < cropped_ny; ++j ) { - cropped_nxmin = std::min( cropped_nx[j], cropped_nxmin ); - cropped_nxmax = std::max( cropped_nx[j], cropped_nxmax ); - } - idx_t cropped_npts = std::accumulate( cropped_nx.begin(), cropped_nx.end(), idx_t{0} ); + for ( idx_t j = 1; j < cropped_ny; ++j ) { + cropped_nxmin = std::min( cropped_nx[j], cropped_nxmin ); + cropped_nxmax = std::max( cropped_nx[j], cropped_nxmax ); + } + idx_t cropped_npts = std::accumulate( cropped_nx.begin(), cropped_nx.end(), idx_t{0} ); - Spacing cropped_yspace( - new spacing::CustomSpacing( cropped_ny, cropped_y.data(), {cropped_ymin, cropped_ymax} ) ); + Spacing cropped_yspace( + new spacing::CustomSpacing( cropped_ny, cropped_y.data(), {cropped_ymin, cropped_ymax} ) ); - std::vector cropped_xspacings( cropped_ny ); - for ( idx_t j = 0; j < cropped_ny; ++j ) { - cropped_xspacings[j] = new spacing::LinearSpacing( cropped_xmin[j], cropped_xmax[j], cropped_nx[j], true ); - } - XSpace cropped_xspace( cropped_xspacings ); + std::vector cropped_xspacings( cropped_ny ); + for ( idx_t j = 0; j < cropped_ny; ++j ) { + cropped_xspacings[j] = + new spacing::LinearSpacing( cropped_xmin[j], cropped_xmax[j], cropped_nx[j], endpoint ); + } + XSpace cropped_xspace( cropped_xspacings ); - for ( idx_t j = 0; j < cropped_ny; ++j ) { - ATLAS_ASSERT( eckit::types::is_approximately_equal( cropped_xspace.xmin()[j], cropped_xmin[j] ) ); - ATLAS_ASSERT( eckit::types::is_approximately_equal( cropped_xspace.xmax()[j], cropped_xmax[j] ) ); - ATLAS_ASSERT( cropped_xspace.nxmin() == cropped_nxmin ); - ATLAS_ASSERT( cropped_xspace.nxmax() == cropped_nxmax ); - } + for ( idx_t j = 0; j < cropped_ny; ++j ) { + ATLAS_ASSERT( eckit::types::is_approximately_equal( cropped_xspace.xmin()[j], cropped_xmin[j] ) ); + if ( cropped_nx[j] > 1 ) { + ATLAS_ASSERT( eckit::types::is_approximately_equal( cropped_xspace.dx()[j], cropped_dx[j], 1.e-10 ) ); + } + ATLAS_ASSERT( cropped_xspace.nx()[j] == cropped_nx[j] ); + ATLAS_ASSERT( cropped_xspace.nxmin() == cropped_nxmin ); + ATLAS_ASSERT( cropped_xspace.nxmax() == cropped_nxmax ); + } - // Modify grid - { - domain_ = dom; - yspace_ = cropped_yspace; - xspace_ = cropped_xspace; - xmin_ = cropped_xmin; - xmax_ = cropped_xmax; - dx_ = cropped_dx; - nx_ = cropped_nx; - nxmin_ = cropped_nxmin; - nxmax_ = cropped_nxmax; - npts_ = cropped_npts; - y_ = cropped_y; + // Modify grid + { + yspace_ = cropped_yspace; + xspace_ = cropped_xspace; + xmin_ = cropped_xmin; + xmax_ = cropped_xmax; + dx_ = cropped_dx; + nx_ = cropped_nx; + nxmin_ = cropped_nxmin; + nxmax_ = cropped_nxmax; + npts_ = cropped_npts; + y_ = cropped_y; + } + } + else { + std::stringstream errmsg; + errmsg << "Cannot crop the grid with domain " << dom; + throw_Exception( errmsg.str(), Here() ); } } - else { - std::stringstream errmsg; - errmsg << "Cannot crop the grid with domain " << dom; - throw_Exception( errmsg.str(), Here() ); - } + domain_ = RectangularDomain( {xspace_.min(),xspace_.max()}, {yspace_.min(),yspace_.max()}, projection_.units() ); } void Structured::computeTruePeriodicity() { - if ( projection_.strictlyRegional() ) { periodic_x_ = false; } - else if ( domain_ && domain_.global() ) { - periodic_x_ = true; + if ( projection_.strictlyRegional() ) { + periodic_x_ = false; + return; + } + if ( not ZonalBandDomain( domain_ ) ) { + periodic_x_ = false; + return; } - else { - // domain could be zonal band - idx_t j = ny() / 2; - if ( std::abs( xmin_[j] + ( nx_[j] - 1 ) * dx_[j] - xmax_[j] ) < 1.e-11 ) { - periodic_x_ = false; // This would lead to duplicated points - } - else { - // High chance to be periodic. Check anyway. - const PointLonLat Pllmin = projection().lonlat( PointXY( xmin_[j], y_[j] ) ); - const PointLonLat Pllmax = projection().lonlat( PointXY( xmax_[j], y_[j] ) ); + idx_t j = ny() / 2; + const PointLonLat Pllmin = projection().lonlat( PointXY( xmin_[j], y_[j] ) ); + const PointLonLat Pllmax = projection().lonlat( PointXY( xmax_[j], y_[j] ) ); - Point3 Pxmin; - util::UnitSphere::convertSphericalToCartesian( Pllmin, Pxmin ); + Point3 Pxmin; + util::UnitSphere::convertSphericalToCartesian( Pllmin, Pxmin ); - Point3 Pxmax; - util::UnitSphere::convertSphericalToCartesian( Pllmax, Pxmax ); + Point3 Pxmax; + util::UnitSphere::convertSphericalToCartesian( Pllmax, Pxmax ); - periodic_x_ = points_equal( Pxmin, Pxmax ); - } - } + periodic_x_ = points_equal( Pxmin, Pxmax ); } void Structured::print( std::ostream& os ) const { From 40aca6b462e18f5fcce0598e2d4dfa2d2e549c0e Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 4 Jun 2019 09:10:52 +0100 Subject: [PATCH 09/11] Apply clang-format --- src/atlas/array/gridtools/GridToolsArray.cc | 3 +- .../array/gridtools/GridToolsArrayHelpers.h | 37 ++++------ .../array/gridtools/GridToolsMakeView.cc | 73 +++++++++---------- src/atlas/array/gridtools/GridToolsMakeView.h | 3 +- src/atlas/domain/Domain.cc | 5 +- src/atlas/domain/detail/GlobalDomain.cc | 11 +-- src/atlas/domain/detail/ZonalBandDomain.cc | 6 +- src/atlas/functionspace/StructuredColumns.cc | 12 +-- .../detail/StructuredColumns_setup.cc | 12 +-- src/atlas/grid/detail/grid/Structured.cc | 17 ++--- src/atlas/trans/local/TransLocal.cc | 11 ++- 11 files changed, 84 insertions(+), 106 deletions(-) diff --git a/src/atlas/array/gridtools/GridToolsArray.cc b/src/atlas/array/gridtools/GridToolsArray.cc index 538b1b194..6da313159 100644 --- a/src/atlas/array/gridtools/GridToolsArray.cc +++ b/src/atlas/array/gridtools/GridToolsArray.cc @@ -64,7 +64,7 @@ class ArrayT_impl { template > void construct( UInts... dims ) { - static_assert( sizeof...(UInts) > 0, "1" ); + static_assert( sizeof...( UInts ) > 0, "1" ); auto gt_storage = create_gt_storage::type>( dims... ); using data_store_t = typename std::remove_pointer::type; array_.data_store_ = std::unique_ptr( new GridToolsDataStore( gt_storage ) ); @@ -571,4 +571,3 @@ template Array* Array::wrap( long unsigned*, const ArraySpec& ); } // namespace array } // namespace atlas - diff --git a/src/atlas/array/gridtools/GridToolsArrayHelpers.h b/src/atlas/array/gridtools/GridToolsArrayHelpers.h index 4f4e064e3..c058163cb 100644 --- a/src/atlas/array/gridtools/GridToolsArrayHelpers.h +++ b/src/atlas/array/gridtools/GridToolsArrayHelpers.h @@ -148,21 +148,15 @@ struct get_pack_size { using type = ::gridtools::static_uint; }; -template +template struct gt_storage_t { - using type = gridtools::storage_traits::data_store_t< - Value, - gridtools::storage_traits::custom_layout_storage_info_t< - 0, - LayoutMap, - ::gridtools::zero_halo - > - >; + using type = gridtools::storage_traits::data_store_t< + Value, gridtools::storage_traits::custom_layout_storage_info_t<0, LayoutMap, ::gridtools::zero_halo>>; }; template -typename gt_storage_t::type::value>::type* -create_gt_storage( UInts... dims ) { +typename gt_storage_t::type::value>::type* create_gt_storage( + UInts... dims ) { static_assert( ( sizeof...( dims ) > 0 ), "Error: can not create storages without any dimension" ); constexpr static unsigned int rank = get_pack_size::type::value; @@ -171,25 +165,22 @@ create_gt_storage( UInts... dims ) { typedef gridtools::storage_traits::data_store_t data_store_t; data_store_t* ds; - if (::gridtools::accumulate(::gridtools::multiplies(), dims... ) == 0 ) { ds = new data_store_t(); } - else { - storage_info_ty si( dims... ); - ds = new data_store_t( si ); - } + if (::gridtools::accumulate(::gridtools::multiplies(), dims... ) == 0 ) { ds = new data_store_t(); } + else { + storage_info_ty si( dims... ); + ds = new data_store_t( si ); + } return ds; } -template +template struct gt_wrap_storage_t { - using type = gridtools::storage_traits::data_store_t< - Value, - gridtools::storage_traits::storage_info_t<0, Rank> - >; + using type = gridtools::storage_traits::data_store_t>; }; template -static typename gt_wrap_storage_t::type* -wrap_gt_storage( Value* data, std::array&& shape, std::array&& strides ) { +static typename gt_wrap_storage_t::type* wrap_gt_storage( Value* data, std::array&& shape, + std::array&& strides ) { static_assert( ( Rank > 0 ), "Error: can not create storages without any dimension" ); typedef gridtools::storage_traits::storage_info_t<0, Rank, ::gridtools::zero_halo> storage_info_ty; typedef gridtools::storage_traits::data_store_t data_store_t; diff --git a/src/atlas/array/gridtools/GridToolsMakeView.cc b/src/atlas/array/gridtools/GridToolsMakeView.cc index 44dc17fb8..b6115a5dd 100644 --- a/src/atlas/array/gridtools/GridToolsMakeView.cc +++ b/src/atlas/array/gridtools/GridToolsMakeView.cc @@ -14,12 +14,12 @@ #include #include "atlas/array/Array.h" -#include "atlas/array/ArrayViewDefs.h" #include "atlas/array/ArrayView.h" +#include "atlas/array/ArrayViewDefs.h" #include "atlas/array/IndexView.h" #include "atlas/array/gridtools/GridToolsTraits.h" -#include "atlas/runtime/Exception.h" #include "atlas/library/config.h" +#include "atlas/runtime/Exception.h" //------------------------------------------------------------------------------ @@ -46,8 +46,7 @@ namespace gridtools { template -typename gt_view::type -make_gt_host_view( const Array& array ) { +typename gt_view::type make_gt_host_view( const Array& array ) { using storage_info_ty = storage_traits::storage_info_t<0, Rank>; using data_store_t = storage_traits::data_store_t; @@ -56,7 +55,7 @@ make_gt_host_view( const Array& array ) { } template -typename gt_view::type make_gt_device_view( const Array& array ) { +typename gt_view::type make_gt_device_view( const Array& array ) { using storage_info_ty = storage_traits::storage_info_t<0, Rank>; using data_store_t = storage_traits::data_store_t; @@ -172,46 +171,46 @@ namespace array { template IndexView make_host_indexview( const Array& ); \ template IndexView make_host_indexview( const Array& ); \ namespace gridtools { \ - template typename gt_view::type \ - make_gt_host_view( const Array& array ); \ - template typename gt_view::type \ - make_gt_host_view( const Array& array ); \ - template typename gt_view::type \ - make_gt_host_view( const Array& array ); \ - template typename gt_view::type \ - make_gt_host_view( const Array& array ); \ - template typename gt_view::type \ + template typename gt_view::type make_gt_host_view( \ + const Array& array ); \ + template typename gt_view::type make_gt_host_view( \ + const Array& array ); \ + template typename gt_view::type make_gt_host_view( \ + const Array& array ); \ + template typename gt_view::type make_gt_host_view( \ + const Array& array ); \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ - template typename gt_view::type \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ - template typename gt_view::type \ - make_gt_host_view( const Array& array ); \ - template typename gt_view::type \ - make_gt_host_view( const Array& array ); \ - template typename gt_view::type \ - make_gt_host_view( const Array& array ); \ - template typename gt_view::type \ + template typename gt_view::type make_gt_host_view( \ + const Array& array ); \ + template typename gt_view::type make_gt_host_view( \ + const Array& array ); \ + template typename gt_view::type make_gt_host_view( \ + const Array& array ); \ + template typename gt_view::type \ make_gt_host_view( const Array& array ); \ \ - template typename gt_view::type \ - make_gt_device_view( const Array& array ); \ - template typename gt_view::type \ - make_gt_device_view( const Array& array ); \ - template typename gt_view::type \ - make_gt_device_view( const Array& array ); \ - template typename gt_view::type \ - make_gt_device_view( const Array& array ); \ - template typename gt_view::type \ + template typename gt_view::type make_gt_device_view( \ + const Array& array ); \ + template typename gt_view::type make_gt_device_view( \ + const Array& array ); \ + template typename gt_view::type make_gt_device_view( \ + const Array& array ); \ + template typename gt_view::type make_gt_device_view( \ + const Array& array ); \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template typename gt_view::type \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template typename gt_view::type \ - make_gt_device_view( const Array& array ); \ - template typename gt_view::type \ + template typename gt_view::type make_gt_device_view( \ + const Array& array ); \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template typename gt_view::type \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ - template typename gt_view::type \ + template typename gt_view::type \ make_gt_device_view( const Array& array ); \ } diff --git a/src/atlas/array/gridtools/GridToolsMakeView.h b/src/atlas/array/gridtools/GridToolsMakeView.h index 735563a7f..325c0e0fc 100644 --- a/src/atlas/array/gridtools/GridToolsMakeView.h +++ b/src/atlas/array/gridtools/GridToolsMakeView.h @@ -19,7 +19,8 @@ namespace gridtools { template struct gt_view { using type = ::gridtools::data_view< - gridtools::storage_traits::data_store_t>, gridtools::get_access_mode(AccessMode) >; + gridtools::storage_traits::data_store_t>, + gridtools::get_access_mode( AccessMode )>; }; template diff --git a/src/atlas/domain/Domain.cc b/src/atlas/domain/Domain.cc index 9aad26338..bc2a869ae 100644 --- a/src/atlas/domain/Domain.cc +++ b/src/atlas/domain/Domain.cc @@ -26,8 +26,7 @@ RectangularDomain::RectangularDomain( const Interval& x, const Interval& y, cons Domain( ( RD::is_global( x, y, units ) ) ? new atlas::domain::GlobalDomain( x[0] ) : ( RD::is_zonal_band( x, units ) ? new atlas::domain::ZonalBandDomain( y, x[0] ) - : new atlas::domain::RectangularDomain( x, y, units ) ) ) { -} + : new atlas::domain::RectangularDomain( x, y, units ) ) ) {} RectangularDomain::RectangularDomain( const Domain& domain ) : Domain( domain ), @@ -58,7 +57,7 @@ double RectangularDomain::ymax() const { } ZonalBandDomain::ZonalBandDomain( const Interval& y, const double& west ) : - RectangularDomain( ( ZD::is_global( y ) ) ? new atlas::domain::GlobalDomain(west) + RectangularDomain( ( ZD::is_global( y ) ) ? new atlas::domain::GlobalDomain( west ) : new atlas::domain::ZonalBandDomain( y, west ) ) {} ZonalBandDomain::ZonalBandDomain( const Domain& domain ) : diff --git a/src/atlas/domain/detail/GlobalDomain.cc b/src/atlas/domain/detail/GlobalDomain.cc index 57a7137e9..b7eeffc01 100644 --- a/src/atlas/domain/detail/GlobalDomain.cc +++ b/src/atlas/domain/detail/GlobalDomain.cc @@ -25,7 +25,7 @@ constexpr std::array yrange() { static double get_west( const eckit::Parametrisation& params ) { double west = 0.; - params.get("west",west); + params.get( "west", west ); return west; } @@ -35,15 +35,12 @@ GlobalDomain::GlobalDomain( const double west ) : ZonalBandDomain( yrange(), wes GlobalDomain::GlobalDomain() : ZonalBandDomain( yrange() ) {} -GlobalDomain::GlobalDomain( const eckit::Parametrisation& params ) : - GlobalDomain( get_west( params ) ) {} +GlobalDomain::GlobalDomain( const eckit::Parametrisation& params ) : GlobalDomain( get_west( params ) ) {} GlobalDomain::Spec GlobalDomain::spec() const { Spec domain_spec; domain_spec.set( "type", type() ); - if( xmin() != 0. ) { - domain_spec.set( "west", xmin() ); - } + if ( xmin() != 0. ) { domain_spec.set( "west", xmin() ); } return domain_spec; } @@ -52,7 +49,7 @@ void GlobalDomain::hash( eckit::Hash& h ) const { } void GlobalDomain::print( std::ostream& os ) const { - os << "GlobalDomain[west=" << xmin() << "]" ; + os << "GlobalDomain[west=" << xmin() << "]"; } namespace { diff --git a/src/atlas/domain/detail/ZonalBandDomain.cc b/src/atlas/domain/detail/ZonalBandDomain.cc index 507ae1e74..902e5c84f 100644 --- a/src/atlas/domain/detail/ZonalBandDomain.cc +++ b/src/atlas/domain/detail/ZonalBandDomain.cc @@ -36,7 +36,7 @@ static std::array get_interval_y( const eckit::Parametrisation& param static double get_west( const eckit::Parametrisation& params ) { double west = 0.; - params.get("west",west); + params.get( "west", west ); return west; } @@ -75,9 +75,7 @@ ZonalBandDomain::Spec ZonalBandDomain::spec() const { domain_spec.set( "type", type() ); domain_spec.set( "ymin", ymin() ); domain_spec.set( "ymax", ymax() ); - if( xmin() != 0. ) { - domain_spec.set( "west", xmin() ); - } + if ( xmin() != 0. ) { domain_spec.set( "west", xmin() ); } return domain_spec; } diff --git a/src/atlas/functionspace/StructuredColumns.cc b/src/atlas/functionspace/StructuredColumns.cc index 08eb84aff..5f05a31c0 100644 --- a/src/atlas/functionspace/StructuredColumns.cc +++ b/src/atlas/functionspace/StructuredColumns.cc @@ -236,9 +236,7 @@ void StructuredColumns::set_field_metadata( const eckit::Configuration& config, config.get( "variables", variables ); field.set_variables( variables ); - if( config.has("type") ){ - field.metadata().set("type",config.getString("type")); - } + if ( config.has( "type" ) ) { field.metadata().set( "type", config.getString( "type" ) ); } } array::DataType StructuredColumns::config_datatype( const eckit::Configuration& config ) const { @@ -416,11 +414,9 @@ Field StructuredColumns::createField( const eckit::Configuration& options ) cons } Field StructuredColumns::createField( const Field& other, const eckit::Configuration& config ) const { - return createField( option::datatype( other.datatype() ) | - option::levels( other.levels() ) | - option::variables( other.variables() ) | - option::type( other.metadata().getString("type","scalar") ) | - config ); + return createField( option::datatype( other.datatype() ) | option::levels( other.levels() ) | + option::variables( other.variables() ) | + option::type( other.metadata().getString( "type", "scalar" ) ) | config ); } // ---------------------------------------------------------------------------- diff --git a/src/atlas/functionspace/detail/StructuredColumns_setup.cc b/src/atlas/functionspace/detail/StructuredColumns_setup.cc index 7395d2759..899dfe0f8 100644 --- a/src/atlas/functionspace/detail/StructuredColumns_setup.cc +++ b/src/atlas/functionspace/detail/StructuredColumns_setup.cc @@ -338,11 +338,11 @@ void StructuredColumns::setup( const grid::Distribution& distribution, const eck ATLAS_TRACE_SCOPE( "Create fields" ) { size_halo_ = gridpoints.size(); field_partition_ = Field( "partition", array::make_datatype(), array::make_shape( size_halo_ ) ); - field_ghost_ = Field( "ghost", array::make_datatype(), array::make_shape( size_halo_ ) ); - field_global_index_ = Field( "glb_idx", array::make_datatype(), array::make_shape( size_halo_ ) ); - field_index_i_ = Field( "index_i", array::make_datatype(), array::make_shape( size_halo_ ) ); - field_index_j_ = Field( "index_j", array::make_datatype(), array::make_shape( size_halo_ ) ); - field_xy_ = Field( "xy", array::make_datatype(), array::make_shape( size_halo_, 2 ) ); + field_ghost_ = Field( "ghost", array::make_datatype(), array::make_shape( size_halo_ ) ); + field_global_index_ = Field( "glb_idx", array::make_datatype(), array::make_shape( size_halo_ ) ); + field_index_i_ = Field( "index_i", array::make_datatype(), array::make_shape( size_halo_ ) ); + field_index_j_ = Field( "index_j", array::make_datatype(), array::make_shape( size_halo_ ) ); + field_xy_ = Field( "xy", array::make_datatype(), array::make_shape( size_halo_, 2 ) ); auto xy = array::make_view( field_xy_ ); auto part = array::make_view( field_partition_ ); @@ -371,7 +371,7 @@ void StructuredColumns::setup( const grid::Distribution& distribution, const eck if ( not in_domain ) { global_idx( gp.r ) = compute_g( gp.i, gp.j ); part( gp.r ) = compute_p( gp.i, gp.j ); - ghost( gp.r ) = 1; + ghost( gp.r ) = 1; } index_i( gp.r ) = gp.i; index_j( gp.r ) = gp.j; diff --git a/src/atlas/grid/detail/grid/Structured.cc b/src/atlas/grid/detail/grid/Structured.cc index e2260e27e..2192b75c8 100644 --- a/src/atlas/grid/detail/grid/Structured.cc +++ b/src/atlas/grid/detail/grid/Structured.cc @@ -86,15 +86,13 @@ Structured::Structured( const std::string& name, XSpace xspace, YSpace yspace, P crop( domain ); computeTruePeriodicity(); - } void Structured::computeDomain() { - if ( periodic() ) { - domain_ = ZonalBandDomain( {yspace().min(),yspace().max()}, xspace().min() ); - } + if ( periodic() ) { domain_ = ZonalBandDomain( {yspace().min(), yspace().max()}, xspace().min() ); } else { - domain_ = RectangularDomain( {xspace().min(),xspace().max()}, {yspace().min(),yspace().max()}, projection_.units() ); + domain_ = RectangularDomain( {xspace().min(), xspace().max()}, {yspace().min(), yspace().max()}, + projection_.units() ); } } @@ -385,8 +383,8 @@ void Structured::crop( const Domain& dom ) { bool endpoint = true; if ( ZonalBandDomain( rect_domain ) ) { for ( idx_t j = 0; j < cropped_ny; ++j ) { - if ( eckit::types::is_approximately_equal( cropped_xmax[j] + cropped_dx[j], - cropped_xmin[j] + 360., 1.e-10 ) ) { + if ( eckit::types::is_approximately_equal( cropped_xmax[j] + cropped_dx[j], cropped_xmin[j] + 360., + 1.e-10 ) ) { cropped_xmax[j] = cropped_xmin[j] + 360.; endpoint = false; } @@ -417,7 +415,8 @@ void Structured::crop( const Domain& dom ) { for ( idx_t j = 0; j < cropped_ny; ++j ) { ATLAS_ASSERT( eckit::types::is_approximately_equal( cropped_xspace.xmin()[j], cropped_xmin[j] ) ); if ( cropped_nx[j] > 1 ) { - ATLAS_ASSERT( eckit::types::is_approximately_equal( cropped_xspace.dx()[j], cropped_dx[j], 1.e-10 ) ); + ATLAS_ASSERT( + eckit::types::is_approximately_equal( cropped_xspace.dx()[j], cropped_dx[j], 1.e-10 ) ); } ATLAS_ASSERT( cropped_xspace.nx()[j] == cropped_nx[j] ); ATLAS_ASSERT( cropped_xspace.nxmin() == cropped_nxmin ); @@ -444,7 +443,7 @@ void Structured::crop( const Domain& dom ) { throw_Exception( errmsg.str(), Here() ); } } - domain_ = RectangularDomain( {xspace_.min(),xspace_.max()}, {yspace_.min(),yspace_.max()}, projection_.units() ); + domain_ = RectangularDomain( {xspace_.min(), xspace_.max()}, {yspace_.min(), yspace_.max()}, projection_.units() ); } void Structured::computeTruePeriodicity() { diff --git a/src/atlas/trans/local/TransLocal.cc b/src/atlas/trans/local/TransLocal.cc index c8f214a32..2b1ec9707 100644 --- a/src/atlas/trans/local/TransLocal.cc +++ b/src/atlas/trans/local/TransLocal.cc @@ -335,16 +335,15 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma gridGlobal_ = grid_; useGlobalLeg = false; } - else { - Log::info() << "Global grid: " << gridGlobal_.spec() << std::endl; - Log::info() << "Actual grid: " << grid_.spec() << std::endl; - std::stringstream log; - log << "Transforms to non-regular regional grids is not supported, unless it defined as a cropping of a global grid." << std::endl; + else { // non-nested reduced grids are not supported + std::ostringstream log; + log << "Transforms to non-regular regional grids is not supported, unless it defined as a cropping of " + "a global grid." + << std::endl; log << " Input grid: " << grid.spec() << std::endl; log << " Applied domain: " << domain << std::endl; log << " Regional grid after Domain applied (the output grid): " << grid_.spec() << std::endl; throw_NotImplemented( log.str(), Here() ); - // non-nested reduced grids are not supported } } From 235bb2509f844a76ae2e098270d9eb48eefb8af0 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 4 Jun 2019 11:42:37 +0100 Subject: [PATCH 10/11] ATLAS-232 Add unit-tests, adapted from MIR (issue MIR-374) --- src/atlas/domain/Domain.cc | 4 + src/atlas/domain/Domain.h | 2 + src/atlas/domain/detail/RectangularDomain.h | 2 + src/tests/grid/CMakeLists.txt | 1 + src/tests/grid/test_grid_cropping.cc | 205 ++++++++++++++++++++ 5 files changed, 214 insertions(+) create mode 100644 src/tests/grid/test_grid_cropping.cc diff --git a/src/atlas/domain/Domain.cc b/src/atlas/domain/Domain.cc index bc2a869ae..9d3475c7f 100644 --- a/src/atlas/domain/Domain.cc +++ b/src/atlas/domain/Domain.cc @@ -40,6 +40,10 @@ bool RectangularDomain::contains_y( double y ) const { return domain_->contains_y( y ); } +bool RectangularDomain::zonal_band() const { + return domain_->zonal_band(); +} + double RectangularDomain::xmin() const { return domain_->xmin(); } diff --git a/src/atlas/domain/Domain.h b/src/atlas/domain/Domain.h index 60abbbab4..685864ef5 100644 --- a/src/atlas/domain/Domain.h +++ b/src/atlas/domain/Domain.h @@ -110,6 +110,8 @@ class RectangularDomain : public Domain { /// Checks if the y-value is contained in the domain bool contains_y( double y ) const; + bool zonal_band() const; + double xmin() const; double xmax() const; double ymin() const; diff --git a/src/atlas/domain/detail/RectangularDomain.h b/src/atlas/domain/detail/RectangularDomain.h index 010d383f2..cccacfd80 100644 --- a/src/atlas/domain/detail/RectangularDomain.h +++ b/src/atlas/domain/detail/RectangularDomain.h @@ -50,6 +50,8 @@ class RectangularDomain : public Domain { virtual bool containsSouthPole() const override; virtual bool global() const override { return global_; } + bool zonal_band() const { return is_zonal_band( {xmin_, xmax_}, units_ ); } + virtual bool empty() const override { return ( xmin_ == xmax_ ) or ( ymin_ == ymax_ ); } virtual Spec spec() const override; diff --git a/src/tests/grid/CMakeLists.txt b/src/tests/grid/CMakeLists.txt index 6cd8e76ae..d3a80914e 100644 --- a/src/tests/grid/CMakeLists.txt +++ b/src/tests/grid/CMakeLists.txt @@ -21,6 +21,7 @@ foreach(test test_field test_grid_ptr test_grids + test_grid_cropping test_vertical test_state test_grid_hash) diff --git a/src/tests/grid/test_grid_cropping.cc b/src/tests/grid/test_grid_cropping.cc new file mode 100644 index 000000000..8cf31136d --- /dev/null +++ b/src/tests/grid/test_grid_cropping.cc @@ -0,0 +1,205 @@ +/* + * (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 +#include +#include + +#include "atlas/grid/Iterator.h" +#include "atlas/grid/StructuredGrid.h" +#include "atlas/mesh/Mesh.h" +#include "atlas/meshgenerator.h" +#include "atlas/output/Gmsh.h" +#include "atlas/runtime/Log.h" +#include "atlas/util/Config.h" + +#include "tests/AtlasTestEnvironment.h" + +using Config = atlas::util::Config; + +namespace atlas { +namespace test { + +//----------------------------------------------------------------------------- + +namespace { +struct Increments { + double dx; + double dy; + Increments( double _dx, double _dy ) : dx( _dx ), dy( _dy ) {} +}; + +static Grid RegularLL( const Increments& increments, const Domain& _domain ) { + // As MIR defines a "RegularLL" grid: + // - If the domain is global or zonal band, the grid is periodic, with the domain's West defining the West boundary + // and the East-most point is at (East-dx) where (East = West+360) + // - If the domain is not a zonal band or global, the East-most point is at (East) + double dx = increments.dx; + double dy = increments.dy; + + RectangularDomain domain{_domain}; + double n = domain.ymax(); + double s = domain.ymin(); + double w = domain.xmin(); + double e = domain.xmax(); + long nj = long( ( n - s ) / dy ) + 1; + long ni = long( ( e - w ) / dx ) + long( domain.zonal_band() ? 0 : 1 ); + + using atlas::grid::LinearSpacing; + StructuredGrid::XSpace xspace( LinearSpacing( w, e, ni, not domain.zonal_band() ) ); + StructuredGrid::YSpace yspace( LinearSpacing( n, s, nj ) ); + + return StructuredGrid( xspace, yspace, StructuredGrid::Projection(), domain ); +} + +static Grid RegularLL( const Increments& increments ) { + return RegularLL( increments, RectangularDomain{{0., 360.}, {-90., 90.}} ); +} + +static Grid RegularLL11( const Domain& domain ) { + return RegularLL( Increments{1., 1.}, domain ); +} + +static Grid RegularLL11() { + return RegularLL( Increments{1., 1.}, RectangularDomain{{0., 360.}, {-90., 90.}} ); +} + +} // namespace + +CASE( "Test number of points for global grids" ) { + // --- Global grids --- // + EXPECT( Grid( "O16" ).size() == 1600 ); + + EXPECT( Grid( "O1280" ).size() == 6599680 ); + + EXPECT( RegularLL( Increments{1., 1.} ).size() == 65160 ); +} + +CASE( "Test number of points for non-global grids" ) { + Grid global = RegularLL11(); + + auto test_regional_size = [&]( const Domain& domain, idx_t size ) -> bool { + if ( RegularLL11( domain ).size() != size ) { + Log::error() << "RegularLL11( domain ).size() != size ---> [ " << RegularLL11( domain ).size() + << " != " << size << " ]" << std::endl; + Log::error() << " with Domain = " << domain << std::endl; + return false; + } + return true; + }; + + auto test_cropping_size = [&]( const Domain& domain, idx_t size ) -> bool { + if ( Grid( global, domain ).size() != size ) { + Log::error() << "Grid( global, domain ).size() != size ---> [ " << Grid( global, domain ).size() + << " != " << size << " ]" << std::endl; + Log::error() << " with Domain = " << domain << std::endl; + return false; + } + return true; + }; + + // -- Parallels + { + // North Pole latitude + EXPECT( test_regional_size( RectangularDomain{{0., 360.}, {90., 90.}}, 360 ) ); + EXPECT( test_cropping_size( RectangularDomain{{0., 360.}, {90., 90.}}, 360 ) ); + + // North Pole latitude plus one extra + EXPECT( test_regional_size( RectangularDomain{{0., 360.}, {89., 90.}}, 720 ) ); + EXPECT( test_cropping_size( RectangularDomain{{0., 360.}, {89., 90.}}, 720 ) ); + + // Equator latitude + EXPECT( test_regional_size( RectangularDomain{{0., 360.}, {0., 0.}}, 360 ) ); + EXPECT( test_cropping_size( RectangularDomain{{0., 360.}, {0., 0.}}, 360 ) ); + + // South Pole latitude + EXPECT( test_regional_size( RectangularDomain{{0., 360.}, {-90., -90.}}, 360 ) ); + EXPECT( test_cropping_size( RectangularDomain{{0., 360.}, {-90., -90.}}, 360 ) ); + + // South Pole latitude plus one extra + EXPECT( test_regional_size( RectangularDomain{{0., 360.}, {-90., -89.}}, 720 ) ); + EXPECT( test_cropping_size( RectangularDomain{{0., 360.}, {-90., -89.}}, 720 ) ); + } + + // -- Meridians + { + // Greenwich Meridian + EXPECT( test_regional_size( RectangularDomain{{0., 0.}, {-90., 90.}}, 181 ) ); + EXPECT( test_cropping_size( RectangularDomain{{0., 0.}, {-90., 90.}}, 181 ) ); + + // Greenwich Meridian + one to East + EXPECT( test_regional_size( RectangularDomain{{0., 1.}, {-90., 90.}}, 2 * 181 ) ); + EXPECT( test_cropping_size( RectangularDomain{{0., 1.}, {-90., 90.}}, 2 * 181 ) ); + + // Greenwich Meridian + one to West + EXPECT( test_regional_size( RectangularDomain{{-1., 0.}, {-90., 90.}}, 2 * 181 ) ); + EXPECT( test_cropping_size( RectangularDomain{{-1., 0.}, {-90., 90.}}, 2 * 181 ) ); + + // Meridian, one degree to the East of Greenwich + EXPECT( test_regional_size( RectangularDomain{{-1., -1.}, {-90., 90.}}, 181 ) ); + EXPECT( test_cropping_size( RectangularDomain{{-1., -1.}, {-90., 90.}}, 181 ) ); + + // Date line at 180 degrees East + EXPECT( test_regional_size( RectangularDomain{{180., 180.}, {-90., 90.}}, 181 ) ); + EXPECT( test_cropping_size( RectangularDomain{{180., 180.}, {-90., 90.}}, 181 ) ); + + // Date line at 180 degrees West + EXPECT( test_regional_size( RectangularDomain{{-180., -180.}, {-90., 90.}}, 181 ) ); + EXPECT( test_cropping_size( RectangularDomain{{-180., -180.}, {-90., 90.}}, 181 ) ); + } + + // Crop to single point {0.,0.} + EXPECT( test_cropping_size( RectangularDomain{{-0.5, 0.5}, {-0.5, 0.5}}, 1 ) ); + EXPECT( test_cropping_size( RectangularDomain{{0., 0.}, {0., 0.}}, 1 ) ); + + // Regional grid with same domain leads to bounds given by domain --> 4 points + EXPECT( test_regional_size( RectangularDomain{{-0.5, 0.5}, {-0.5, 0.5}}, 4 ) ); + EXPECT( test_regional_size( RectangularDomain{{0., 0.}, {0., 0.}}, 1 ) ); +} + +CASE( "Test first point of cropped global grids (MIR-374)" ) { + auto test_first_point = []( const std::string& name, const Domain& domain, const PointLonLat& p ) -> bool { + double eps = 1.e-6; + Grid global_grid{name}; + Grid cropped_grid{global_grid, domain}; + + PointLonLat first_point = *cropped_grid.lonlat().begin(); + if ( not is_approximately_equal( first_point.lon(), p.lon(), eps ) || + not is_approximately_equal( first_point.lat(), p.lat(), eps ) ) { + auto old_precision = Log::error().precision( 16 ); + Log::error() << "First point doesn't match for Grid " << name << " with domain " << domain << std::endl; + Log::error() << " first_point = " << first_point << std::endl; + Log::error() << " expected = " << p << std::endl; + Log::error().precision( old_precision ); + return false; + } + return true; + }; + + // clang-format off + EXPECT( test_first_point( "O16", RectangularDomain{ {-180.,180.}, {-90.,90.} }, PointLonLat{-180.,85.76058712044382} ) ); + EXPECT( test_first_point( "O640", RectangularDomain{ {-180.,180.}, {-90.,90.} }, PointLonLat{-180.,89.89239644558999} ) ); + EXPECT( test_first_point( "O16", RectangularDomain{ {-180.,180.}, {-60.,90.} }, PointLonLat{-180.,85.76058712044382} ) ); + EXPECT( test_first_point( "O640", RectangularDomain{ {-180.,180.}, {-60.,90.} }, PointLonLat{-180.,89.89239644558999} ) ); + EXPECT( test_first_point( "O16", RectangularDomain{ {-180.,180.}, {-90.,60.} }, PointLonLat{-180.,58.14295404920328} ) ); + EXPECT( test_first_point( "O640", RectangularDomain{ {-180.,180.}, {-90.,60.} }, PointLonLat{-180.,59.953135752239} ) ); + EXPECT( test_first_point( "O16", RectangularDomain{ {-10. ,10. }, {-10.,10.} }, PointLonLat{-9.473684210526358,8.306702856518804} ) ); + EXPECT( test_first_point( "O640", RectangularDomain{ {-10. ,10. }, {-10.,10.} }, PointLonLat{-9.878048780487802,9.910190568389} ) ); + // clang-format on +} + +//----------------------------------------------------------------------------- + +} // namespace test +} // namespace atlas + +int main( int argc, char** argv ) { + return atlas::test::run( argc, argv ); +} From 8fdef00caa927029c19d87c6b3f3a4fc07b0cd18 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 4 Jun 2019 15:18:20 +0100 Subject: [PATCH 11/11] Version 0.17.2 --- CHANGELOG.md | 8 ++++++++ VERSION.cmake | 2 +- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 10b68ca50..994fad0a9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,12 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html ## [Unreleased] +## [0.17.2] - 2019-06-04 +### Fixed +- Compilation with PGI 19.4 +- Structured Grid creation for periodic domains that do not start at 0 degrees longitude (Greenwich) + + ## [0.17.1] - 2019-04-22 ### Added - Option to declaration of field type (vector/scalar) when creating field in FunctionSpace @@ -107,6 +113,8 @@ 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.17.2]: https://github.com/ecmwf/atlas/compare/0.17.1...0.17.2 +[0.17.1]: https://github.com/ecmwf/atlas/compare/0.17.0...0.17.1 [0.17.0]: https://github.com/ecmwf/atlas/compare/0.16.0...0.17.0 [0.16.0]: https://github.com/ecmwf/atlas/compare/0.15.2...0.16.0 [0.15.2]: https://github.com/ecmwf/atlas/compare/0.15.1...0.15.2 diff --git a/VERSION.cmake b/VERSION.cmake index f8f5f8c36..f146c319a 100644 --- a/VERSION.cmake +++ b/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.17.1" ) +set ( ${PROJECT_NAME}_VERSION_STR "0.17.2" )