From f81de650ff61462c9b337409447e95d7d9f925bd Mon Sep 17 00:00:00 2001 From: Pedro Maciel Date: Thu, 19 Jul 2018 15:05:24 +0100 Subject: [PATCH 01/21] MIR-283 Fix for almost empty domains --- src/atlas/domain/Domain.cc | 2 +- src/atlas/domain/Domain.h | 2 +- src/atlas/grid/detail/grid/Structured.cc | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/atlas/domain/Domain.cc b/src/atlas/domain/Domain.cc index 8f5a62861..ea4b9583d 100644 --- a/src/atlas/domain/Domain.cc +++ b/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_ ) {} diff --git a/src/atlas/domain/Domain.h b/src/atlas/domain/Domain.h index 07f3550fe..4df3284e4 100644 --- a/src/atlas/domain/Domain.h +++ b/src/atlas/domain/Domain.h @@ -44,7 +44,7 @@ class Domain { Domain( const Implementation* ); Domain( const eckit::Parametrisation& ); - operator bool() { return true; } + operator bool() { return domain_; } /// Type of the domain std::string type() const; diff --git a/src/atlas/grid/detail/grid/Structured.cc b/src/atlas/grid/detail/grid/Structured.cc index 6cb838a1e..f82f5071a 100644 --- a/src/atlas/grid/detail/grid/Structured.cc +++ b/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(); @@ -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 { From 395f837c9e8f64178d29ff681be3d863f0b23882 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 19 Jul 2018 18:07:06 +0100 Subject: [PATCH 02/21] MIR-283 Fix remaining issues --- src/atlas/domain/Domain.h | 2 +- src/atlas/grid/Grid.cc | 4 +++- src/atlas/grid/detail/grid/Structured.cc | 2 +- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/atlas/domain/Domain.h b/src/atlas/domain/Domain.h index 4df3284e4..832a3bc32 100644 --- a/src/atlas/domain/Domain.h +++ b/src/atlas/domain/Domain.h @@ -44,7 +44,7 @@ class Domain { Domain( const Implementation* ); Domain( const eckit::Parametrisation& ); - operator bool() { return domain_; } + operator bool() const { return domain_; } /// Type of the domain std::string type() const; diff --git a/src/atlas/grid/Grid.cc b/src/atlas/grid/Grid.cc index 2a130116b..f2ec85a4b 100644 --- a/src/atlas/grid/Grid.cc +++ b/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 ) { diff --git a/src/atlas/grid/detail/grid/Structured.cc b/src/atlas/grid/detail/grid/Structured.cc index f82f5071a..d491834b5 100644 --- a/src/atlas/grid/detail/grid/Structured.cc +++ b/src/atlas/grid/detail/grid/Structured.cc @@ -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] ); From 0516f7dd24f1311752ebdc0efb8b104c9f4ecb9b Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 19 Jul 2018 18:15:50 +0100 Subject: [PATCH 03/21] MIR-283 Add unit test --- src/tests/grid/test_grids.cc | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/src/tests/grid/test_grids.cc b/src/tests/grid/test_grids.cc index 64f87095d..594124af7 100644 --- a/src/tests/grid/test_grids.cc +++ b/src/tests/grid/test_grids.cc @@ -147,6 +147,34 @@ CASE( "test_cropping previous case" ) { EXPECT( cropped.size() == 267 ); } +CASE( "cropping with line at north pole" ) { + StructuredGrid grid( "L16", RectangularDomain( {0, 360}, {90, 90} ) ); + EXPECT( grid.ny() == 1 ); + EXPECT( grid.nx( 0 ) == 64 ); + EXPECT( grid.size() == 64 ); +} + +CASE( "cropping with line at south pole" ) { + StructuredGrid grid( "L16", RectangularDomain( {0, 360}, {-90, -90} ) ); + EXPECT( grid.ny() == 1 ); + EXPECT( grid.nx( 0 ) == 64 ); + EXPECT( grid.size() == 64 ); +} + +CASE( "cropping with line at equator" ) { + StructuredGrid grid( "L16", RectangularDomain( {0, 360}, {0, 0} ) ); + EXPECT( grid.ny() == 1 ); + EXPECT( grid.nx( 0 ) == 64 ); + EXPECT( grid.size() == 64 ); +} + +CASE( "cropping single point at equator" ) { + StructuredGrid grid( "L16", RectangularDomain( {0, 0}, {0, 0} ) ); + EXPECT( grid.ny() == 1 ); + EXPECT( grid.nx( 0 ) == 1 ); + EXPECT( grid.size() == 1 ); +} + //----------------------------------------------------------------------------- From 282cba886acd7e9ef57d58fda716eca74adf866c Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 20 Jul 2018 10:41:04 +0100 Subject: [PATCH 04/21] MIR-282 Massage default options for TransLocal to deal better with unstructured grids --- src/atlas/option/TransOptions.h | 2 +- src/atlas/trans/local/TransLocal.cc | 36 ++++++++++++++++++++++++----- 2 files changed, 31 insertions(+), 7 deletions(-) diff --git a/src/atlas/option/TransOptions.h b/src/atlas/option/TransOptions.h index bb5c0636c..74e71bd86 100644 --- a/src/atlas/option/TransOptions.h +++ b/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 ); }; // ---------------------------------------------------------------------------- diff --git a/src/atlas/trans/local/TransLocal.cc b/src/atlas/trans/local/TransLocal.cc index b0c62d8de..a683047ba 100644 --- a/src/atlas/trans/local/TransLocal.cc +++ b/src/atlas/trans/local/TransLocal.cc @@ -62,7 +62,7 @@ class TransParameters { bool global() const { return config_.getBool( "global", false ); } - int warning() const { return config_.getLong( "warning", 0 ); } + int warning() const { return config_.getLong( "warning", 1 ); } int fft() const { static const std::map string_to_FFT = {{"OFF", (int)option::FFT::OFF}, @@ -176,9 +176,27 @@ int num_n( const int truncation, const int m, const bool symmetric ) { return len; } +class AllocationFailed : public eckit::Exception { +public: + AllocationFailed( size_t bytes, const eckit::CodeLocation& ); + +private: + static std::string error_message( size_t bytes ) { + std::stringstream ss; + ss << "AllocationFailed: Could not allocate " << eckit::Bytes( bytes ); + return ss.str(); + } +}; + +AllocationFailed::AllocationFailed( size_t bytes, const eckit::CodeLocation& loc ) : + Exception( error_message( bytes ), loc ) {} + + void alloc_aligned( double*& ptr, size_t n ) { const size_t alignment = 64 * sizeof( double ); - posix_memalign( (void**)&ptr, alignment, sizeof( double ) * n ); + size_t bytes = sizeof( double ) * n; + int err = posix_memalign( (void**)&ptr, alignment, bytes ); + if ( err ) { throw AllocationFailed( bytes, Here() ); } } void free_aligned( double*& ptr ) { @@ -272,7 +290,7 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma int nlonsMax = 0; int neqtr = 0; useFFT_ = TransParameters( config ).fft(); - unstruct_precomp_ = true; + unstruct_precomp_ = ( config.has( "precompute" ) ? precompute_ : false ); no_symmetry_ = false; nlatsNH_ = 0; nlatsSH_ = 0; @@ -601,9 +619,15 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma ATLAS_TRACE( "Legendre precomputations (unstructured)" ); if ( warning() ) { - Log::warning() << "WARNING: Precomputations for spectral transforms could take a long time and consume " - "a lot of memory (unstructured grid approach)! Results may contain aliasing errors." - << std::endl; + Log::warning() + << "WARNING: Precomputations for spectral transforms could take a long time as there's no structure" + " to take advantage of!!!" + << std::endl + << "The precomputed values may and consume at least " + << eckit::Bytes( sizeof( double ) * legendre_size( truncation_ ) * grid_.size() ) << " (" + << eckit::Bytes( sizeof( double ) * legendre_size( truncation_ ) ) << " for each of " + << grid_.size() << " grid points )" << std::endl + << "Furthermore, results may contain aliasing errors." << std::endl; } std::vector lats( grid_.size() ); From 432aef886c4fad1d50735db38b8774b34441ddd7 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 20 Jul 2018 16:55:08 +0100 Subject: [PATCH 05/21] correct typo --- src/atlas/trans/local/TransLocal.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/atlas/trans/local/TransLocal.cc b/src/atlas/trans/local/TransLocal.cc index a683047ba..e8b697df6 100644 --- a/src/atlas/trans/local/TransLocal.cc +++ b/src/atlas/trans/local/TransLocal.cc @@ -623,7 +623,7 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma << "WARNING: Precomputations for spectral transforms could take a long time as there's no structure" " to take advantage of!!!" << std::endl - << "The precomputed values may and consume at least " + << "The precomputed values consume at least " << eckit::Bytes( sizeof( double ) * legendre_size( truncation_ ) * grid_.size() ) << " (" << eckit::Bytes( sizeof( double ) * legendre_size( truncation_ ) ) << " for each of " << grid_.size() << " grid points )" << std::endl From a621a08551b81951c9bde62577419bdd30c18c21 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 20 Jul 2018 19:42:05 +0100 Subject: [PATCH 06/21] ATLAS-173 Work on structuredcolumns indexing --- src/atlas/functionspace/StructuredColumns.cc | 14 +++++++++++ src/atlas/functionspace/StructuredColumns.h | 6 +++++ ...functionspace_StructuredColumns_module.F90 | 24 ++++++++++++++++++- .../functionspace/fctest_functionspace.F90 | 22 ++++++++++++----- 4 files changed, 59 insertions(+), 7 deletions(-) diff --git a/src/atlas/functionspace/StructuredColumns.cc b/src/atlas/functionspace/StructuredColumns.cc index 3a10f19ba..d3c392052 100644 --- a/src/atlas/functionspace/StructuredColumns.cc +++ b/src/atlas/functionspace/StructuredColumns.cc @@ -412,10 +412,14 @@ StructuredColumns::StructuredColumns( const Grid& grid, const grid::Partitioner& field_partition_ = Field( "partition", array::make_datatype(), array::make_shape( size_halo_ ) ); field_global_index_ = Field( "glb_idx", array::make_datatype(), array::make_shape( size_halo_ ) ); field_remote_index_ = Field( "remote_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_ ) ); auto part = array::make_view( field_partition_ ); auto global_idx = array::make_view( field_global_index_ ); auto remote_idx = array::make_view( field_remote_index_ ); + auto index_i = array::make_indexview( field_index_i_ ); + auto index_j = array::make_indexview( 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(); +} } // ---------------------------------------------------------------------------- diff --git a/src/atlas/functionspace/StructuredColumns.h b/src/atlas/functionspace/StructuredColumns.h index 1416591b4..9e4859bec 100644 --- a/src/atlas/functionspace/StructuredColumns.h +++ b/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 diff --git a/src/atlas_f/functionspace/atlas_functionspace_StructuredColumns_module.F90 b/src/atlas_f/functionspace/atlas_functionspace_StructuredColumns_module.F90 index 7af294483..47a06f50f 100644 --- a/src/atlas_f/functionspace/atlas_functionspace_StructuredColumns_module.F90 +++ b/src/atlas_f/functionspace/atlas_functionspace_StructuredColumns_module.F90 @@ -68,8 +68,15 @@ module atlas_functionspace_StructuredColumns_module procedure :: i_end_halo procedure :: xy + !! Return xy coordinate field procedure :: partition + !! Return partition field procedure :: global_index + !! Return global_index field + procedure :: index_i + !! Return index_i field + procedure :: index_j + !! Return index_j field procedure, private :: set_index @@ -273,7 +280,6 @@ function partition(this) result(field) call field%return() end function - function global_index(this) result(field) use atlas_functionspace_StructuredColumns_c_binding type(atlas_Field) :: field @@ -282,6 +288,22 @@ function global_index(this) result(field) call field%return() end function +function index_i(this) result(field) + use atlas_functionspace_StructuredColumns_c_binding + type(atlas_Field) :: field + class(atlas_functionspace_StructuredColumns), intent(in) :: this + field = atlas_Field( atlas__fs__StructuredColumns__index_i(this%c_ptr()) ) + call field%return() +end function + +function index_j(this) result(field) + use atlas_functionspace_StructuredColumns_c_binding + type(atlas_Field) :: field + class(atlas_functionspace_StructuredColumns), intent(in) :: this + field = atlas_Field( atlas__fs__StructuredColumns__index_j(this%c_ptr()) ) + call field%return() +end function + subroutine halo_exchange_fieldset(this,fieldset) use atlas_functionspace_StructuredColumns_c_binding class(atlas_functionspace_StructuredColumns), intent(in) :: this diff --git a/src/tests/functionspace/fctest_functionspace.F90 b/src/tests/functionspace/fctest_functionspace.F90 index 67796c11a..294c91627 100644 --- a/src/tests/functionspace/fctest_functionspace.F90 +++ b/src/tests/functionspace/fctest_functionspace.F90 @@ -499,9 +499,13 @@ module fcta_FunctionSpace_fxt integer :: i, j character(len=10) str -type(atlas_Field) field -type(atlas_Field) field_xy -real(8), pointer :: xy(:,:), x(:) +type(atlas_Field) :: field +type(atlas_Field) :: field_xy +type(atlas_Field) :: field_global_index +type(atlas_Field) :: field_index_j +real(8), pointer :: xy(:,:), x(:) +integer(ATLAS_KIND_GIDX), pointer :: global_index(:) +integer(ATLAS_KIND_IDX), pointer :: index_j(:) integer, parameter :: XX=1 grid = atlas_StructuredGrid("O8") @@ -512,14 +516,19 @@ module fcta_FunctionSpace_fxt FCTEST_CHECK_EQUAL( field%owners(), 1 ) field_xy = fs%xy() FCTEST_CHECK_EQUAL( field_xy%owners(), 2 ) -call field%host_data(x) -call field_xy%host_data(xy) +call field%data(x) +call field_xy%data(xy) +field_global_index = fs%global_index() +call field_global_index%data(global_index) + +field_index_j = fs%index_j() +call field_index_j%data(index_j) do j=fs%j_begin_halo(),fs%j_end_halo() write(str,'(I4,A)') j, ' : ' call atlas_log%info(str,newl=.false.) do i=fs%i_begin_halo(j),fs%i_end_halo(j) - write(str,'(I4)') fs%index(i,j) + write(str,'(I4)') global_index( fs%index(i,j) ) call atlas_log%info(str,newl=.false.) x(fs%index(i,j)) = xy(XX,fs%index(i,j)) enddo @@ -547,6 +556,7 @@ module fcta_FunctionSpace_fxt call field%final() FCTEST_CHECK_EQUAL( field_xy%owners(), 1 ) call field_xy%final() +call field_global_index%final() call fs%final() call fs_base%final() call grid%final() From ecfee08b77c59e771db25606f4a1f48b18c2880c Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 23 Jul 2018 13:58:30 +0100 Subject: [PATCH 07/21] ATLAS-171 Implement cropped unstructured grids --- src/atlas/grid/Grid.cc | 4 ++ src/atlas/grid/Grid.h | 1 + src/atlas/grid/detail/grid/Grid.cc | 2 +- src/atlas/grid/detail/grid/Unstructured.cc | 62 ++++++++++++++++++++++ src/atlas/grid/detail/grid/Unstructured.h | 5 ++ src/tests/grid/test_grids.cc | 17 ++++++ 6 files changed, 90 insertions(+), 1 deletion(-) diff --git a/src/atlas/grid/Grid.cc b/src/atlas/grid/Grid.cc index f2ec85a4b..aba813870 100644 --- a/src/atlas/grid/Grid.cc +++ b/src/atlas/grid/Grid.cc @@ -76,6 +76,10 @@ UnstructuredGrid::UnstructuredGrid( std::initializer_list 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( grid ); } diff --git a/src/atlas/grid/Grid.h b/src/atlas/grid/Grid.h index 69916b108..b3e12fd37 100644 --- a/src/atlas/grid/Grid.h +++ b/src/atlas/grid/Grid.h @@ -152,6 +152,7 @@ class UnstructuredGrid : public Grid { UnstructuredGrid( std::vector* ); // takes ownership UnstructuredGrid( std::vector&& ); // move constructor UnstructuredGrid( std::initializer_list ); + UnstructuredGrid( const Grid&, const Domain& ); // Create a new unstructured grid! operator bool() const { return valid(); } diff --git a/src/atlas/grid/detail/grid/Grid.cc b/src/atlas/grid/detail/grid/Grid.cc index 1b887421c..e65a65dc3 100644 --- a/src/atlas/grid/detail/grid/Grid.cc +++ b/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 ); } } diff --git a/src/atlas/grid/detail/grid/Unstructured.cc b/src/atlas/grid/detail/grid/Unstructured.cc index bdb5bf208..e9b867011 100644 --- a/src/atlas/grid/detail/grid/Unstructured.cc +++ b/src/atlas/grid/detail/grid/Unstructured.cc @@ -11,13 +11,17 @@ #include "atlas/grid/detail/grid/Unstructured.h" #include +#include #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

( xmin_, x, eps_ ) ) { + x += 360.; + } + while ( eckit::types::is_strictly_greater( 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 ); + 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" ); diff --git a/src/atlas/grid/detail/grid/Unstructured.h b/src/atlas/grid/detail/grid/Unstructured.h index 135b6d726..4b985e92b 100644 --- a/src/atlas/grid/detail/grid/Unstructured.h +++ b/src/atlas/grid/detail/grid/Unstructured.h @@ -16,6 +16,7 @@ #pragma once #include +#include #include #include @@ -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& ); diff --git a/src/tests/grid/test_grids.cc b/src/tests/grid/test_grids.cc index 594124af7..9d96cd777 100644 --- a/src/tests/grid/test_grids.cc +++ b/src/tests/grid/test_grids.cc @@ -25,6 +25,7 @@ #include "tests/AtlasTestEnvironment.h" using StructuredGrid = atlas::grid::StructuredGrid; +using UnstructuredGrid = atlas::grid::UnstructuredGrid; using Grid = atlas::Grid; using Regular = atlas::grid::RegularGrid; using ReducedGaussianGrid = atlas::grid::ReducedGaussianGrid; @@ -175,6 +176,22 @@ CASE( "cropping single point at equator" ) { EXPECT( grid.size() == 1 ); } +CASE( "Create cropped unstructured grid using rectangular domain" ) { + StructuredGrid agrid( "L8" ); + auto domain = RectangularDomain( {-27, 45}, {33, 73} ); + StructuredGrid sgrid( agrid, domain ); + UnstructuredGrid ugrid( agrid, domain ); + EXPECT( ugrid.size() == sgrid.size() ); +} + +CASE( "Create cropped unstructured grid using zonal domain" ) { + StructuredGrid agrid( "L8" ); + auto domain = ZonalBandDomain( {33, 73} ); + StructuredGrid sgrid( agrid, domain ); + UnstructuredGrid ugrid( agrid, domain ); + EXPECT( ugrid.size() == sgrid.size() ); +} + //----------------------------------------------------------------------------- From 500fc0743a8cc16b2367b69b653d75ccfd38fde1 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 23 Jul 2018 14:29:09 +0100 Subject: [PATCH 08/21] ATLAS-171 Another unit-test --- src/tests/grid/test_grids.cc | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/tests/grid/test_grids.cc b/src/tests/grid/test_grids.cc index 9d96cd777..fdc049704 100644 --- a/src/tests/grid/test_grids.cc +++ b/src/tests/grid/test_grids.cc @@ -192,6 +192,16 @@ CASE( "Create cropped unstructured grid using zonal domain" ) { EXPECT( ugrid.size() == sgrid.size() ); } +CASE( "Create unstructured from unstructured" ) { + StructuredGrid agrid( "L8" ); + UnstructuredGrid global_unstructured( agrid, Domain() ); + EXPECT( UnstructuredGrid( global_unstructured ) ); + EXPECT( global_unstructured.size() == agrid.size() ); + auto domain = ZonalBandDomain( {33, 73} ); + UnstructuredGrid ugrid( global_unstructured, domain ); + EXPECT( ugrid.size() == StructuredGrid( agrid, domain ).size() ); +} + //----------------------------------------------------------------------------- From 6fe83410bb37c8a458a5ea7eab9809044a9a16d4 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 23 Jul 2018 14:45:25 +0100 Subject: [PATCH 09/21] ATLAS-162 Enable unstructured test --- src/tests/trans/test_transgeneral.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/tests/trans/test_transgeneral.cc b/src/tests/trans/test_transgeneral.cc index e2ccea7de..73c3942a3 100644 --- a/src/tests/trans/test_transgeneral.cc +++ b/src/tests/trans/test_transgeneral.cc @@ -794,7 +794,7 @@ CASE( "test_trans_domain" ) { #endif //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- -#if 0 +#if 1 CASE( "test_trans_unstructured" ) { Log::info() << "test_trans_unstructured" << std::endl; // test transgeneral by comparing with analytic solution on an unstructured grid @@ -838,7 +838,9 @@ CASE( "test_trans_unstructured" ) { std::vector rgp_analytic2( gu.size() ); trans::Trans transLocal1( grid_global, testdomain, trc, option::type( "local" ) ); - trans::Trans transLocal2( gu, trc, util::Config( "type", "local" ) ); + + // ATLAS-173 : This should also work with precompute = true, ans should give same. + trans::Trans transLocal2( gu, trc, option::type( "local" ) | util::Config( "precompute", false ) ); int icase = 0; for ( int ivar_in = 2; ivar_in < 3; ivar_in++ ) { // vorticity, divergence, scalar From 3af551e0cd67b5369298cacbfef359cdb359f7c2 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 24 Jul 2018 19:31:35 +0100 Subject: [PATCH 10/21] Remove PGI compiler warnings --- src/atlas/trans/ifs/TransIFS.cc | 12 ++++----- src/atlas/trans/local/TransLocal.cc | 42 ++++++++++++++--------------- src/tests/mesh/fctest_meshgen.F90 | 1 + 3 files changed, 27 insertions(+), 28 deletions(-) diff --git a/src/atlas/trans/ifs/TransIFS.cc b/src/atlas/trans/ifs/TransIFS.cc index f49dce544..feefd3c20 100644 --- a/src/atlas/trans/ifs/TransIFS.cc +++ b/src/atlas/trans/ifs/TransIFS.cc @@ -311,7 +311,7 @@ struct PackNodeColumns { default: ATLAS_DEBUG_VAR( field.rank() ); NOTIMP; - break; + //break; } } @@ -375,7 +375,7 @@ struct PackStructuredColumns { default: ATLAS_DEBUG_VAR( field.rank() ); NOTIMP; - break; + //break; } } @@ -418,7 +418,7 @@ struct PackSpectral { default: ATLAS_DEBUG_VAR( field.rank() ); NOTIMP; - break; + //break; } } @@ -468,7 +468,7 @@ struct UnpackNodeColumns { default: ATLAS_DEBUG_VAR( field.rank() ); NOTIMP; - break; + //break; } } @@ -532,7 +532,7 @@ struct UnpackStructuredColumns { default: ATLAS_DEBUG_VAR( field.rank() ); NOTIMP; - break; + //break; } } @@ -575,7 +575,7 @@ struct UnpackSpectral { default: ATLAS_DEBUG_VAR( field.rank() ); NOTIMP; - break; + //break; } } diff --git a/src/atlas/trans/local/TransLocal.cc b/src/atlas/trans/local/TransLocal.cc index e8b697df6..6856c0a9a 100644 --- a/src/atlas/trans/local/TransLocal.cc +++ b/src/atlas/trans/local/TransLocal.cc @@ -44,23 +44,27 @@ class TransParameters { TransParameters( const eckit::Configuration& config ) : config_( config ) {} ~TransParameters() {} - bool scalar_derivatives() const { return config_.getBool( "scalar_derivatives", false ); } + /* + * For the future + */ + // bool scalar_derivatives() const { return config_.getBool( "scalar_derivatives", false ); } - bool wind_EW_derivatives() const { return config_.getBool( "wind_EW_derivatives", false ); } + // bool wind_EW_derivatives() const { return config_.getBool( "wind_EW_derivatives", false ); } - bool vorticity_divergence_fields() const { return config_.getBool( "vorticity_divergence_fields", false ); } + // bool vorticity_divergence_fields() const { return config_.getBool( "vorticity_divergence_fields", false ); } - std::string read_legendre() const { return config_.getString( "read_legendre", "" ); } + // std::string read_legendre() const { return config_.getString( "read_legendre", "" ); } + // bool global() const { return config_.getBool( "global", false ); } - std::string write_legendre() const { return config_.getString( "write_legendre", "" ); } + // std::string read_fft() const { return config_.getString( "read_fft", "" ); } - bool export_legendre() const { return config_.getBool( "export_legendre", false ); } - std::string read_fft() const { return config_.getString( "read_fft", "" ); } + std::string write_legendre() const { return config_.getString( "write_legendre", "" ); } std::string write_fft() const { return config_.getString( "write_fft", "" ); } - bool global() const { return config_.getBool( "global", false ); } + + bool export_legendre() const { return config_.getBool( "export_legendre", false ); } int warning() const { return config_.getLong( "warning", 1 ); } @@ -91,12 +95,6 @@ struct ReadCache { return v; } - Grid read_grid() { - long& size = *read( 1 ); - char* json = read( size ); - return Grid( eckit::YAMLConfiguration( std::string( json, size ) ) ); - } - char* begin; size_t pos; }; @@ -145,10 +143,10 @@ struct FFTW_Wisdom { FFTW_Wisdom() { wisdom = fftw_export_wisdom_to_string(); } ~FFTW_Wisdom() { free( wisdom ); } }; -std::ostream& operator<<( std::ostream& out, const FFTW_Wisdom& w ) { - out << w.wisdom; - return out; -} +//std::ostream& operator<<( std::ostream& out, const FFTW_Wisdom& w ) { +// out << w.wisdom; +// return out; +//} #endif } // namespace @@ -162,10 +160,10 @@ size_t legendre_size( const size_t truncation ) { return ( truncation + 2 ) * ( truncation + 1 ) / 2; } -int nlats_northernHemisphere( const int nlats ) { - return ceil( nlats / 2. ); - // using ceil here should make it possible to have odd number of latitudes (with the centre latitude being the equator) -} +//int nlats_northernHemisphere( const int nlats ) { +// return ceil( nlats / 2. ); +// // using ceil here should make it possible to have odd number of latitudes (with the centre latitude being the equator) +//} int num_n( const int truncation, const int m, const bool symmetric ) { int len = 0; diff --git a/src/tests/mesh/fctest_meshgen.F90 b/src/tests/mesh/fctest_meshgen.F90 index 50446c5f2..b0b67c071 100644 --- a/src/tests/mesh/fctest_meshgen.F90 +++ b/src/tests/mesh/fctest_meshgen.F90 @@ -125,6 +125,7 @@ TEST( test_meshgen ) + use, intrinsic :: iso_c_binding use atlas_module implicit none type(atlas_StructuredGrid) :: grid From 123f5d17b581302736fa51131e34b843346ecc0d Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 30 Jul 2018 14:05:40 +0000 Subject: [PATCH 11/21] ATLAS-164 fix failing unit tests on cray --- src/atlas/trans/ifs/TransIFS.cc | 3 ++- src/tests/trans/fctest_trans.F90 | 5 +++++ src/tests/trans/test_trans.cc | 10 +++++++++- 3 files changed, 16 insertions(+), 2 deletions(-) diff --git a/src/atlas/trans/ifs/TransIFS.cc b/src/atlas/trans/ifs/TransIFS.cc index feefd3c20..c1ec7c786 100644 --- a/src/atlas/trans/ifs/TransIFS.cc +++ b/src/atlas/trans/ifs/TransIFS.cc @@ -796,7 +796,8 @@ void TransIFS::__dirtrans( const functionspace::NodeColumns& gp, const FieldSet& transform.nscalar = nfld; transform.rgp = rgp.data(); transform.rspscalar = rspec.data(); - + transform.ngpblks = ngptot(); + transform.nproma = 1; TRANS_CHECK(::trans_dirtrans( &transform ) ); } diff --git a/src/tests/trans/fctest_trans.F90 b/src/tests/trans/fctest_trans.F90 index 77d47049e..d5b118c20 100644 --- a/src/tests/trans/fctest_trans.F90 +++ b/src/tests/trans/fctest_trans.F90 @@ -160,6 +160,7 @@ end module fctest_atlas_trans_fixture windfield = nodes_fs%create_field(name="wind",kind=atlas_real(c_double),levels=nlev,variables=3) call windfield%data(wind) + wind(:,:,:) = 0._c_double write(0,*) "nodes_fs%owners()",nodes_fs%owners() vorfield = spectral_fs%create_field(name="vorticity",kind=atlas_real(c_double),levels=nlev) @@ -313,6 +314,7 @@ end module fctest_atlas_trans_fixture integer :: jfld, nfld character(len=10) :: fieldname real(c_double) :: norm +real(c_double), pointer :: gvar(:) grid = atlas_StructuredGrid("O24") trans = atlas_Trans(grid,23) @@ -329,6 +331,9 @@ end module fctest_atlas_trans_fixture ! Read global field data ! ... + FCTEST_CHECK_EQUAL( fieldg%rank(), 1 ) + call fieldg%data(gvar) + gvar(:) = 0. call gridpoints%scatter(fieldg,field) diff --git a/src/tests/trans/test_trans.cc b/src/tests/trans/test_trans.cc index b81037a30..4c2e3b03a 100644 --- a/src/tests/trans/test_trans.cc +++ b/src/tests/trans/test_trans.cc @@ -292,6 +292,8 @@ CASE( "test_spectral_fields" ) { Field spf = spectral.createField( option::name( "spf" ) ); Field gpf = nodal.createField( option::name( "gpf" ) ); + + array::make_view( gpf ).assign(0); EXPECT_NO_THROW( trans.dirtrans( gpf, spf ) ); EXPECT_NO_THROW( trans.invtrans( spf, gpf ) ); @@ -386,6 +388,8 @@ CASE( "test_trans_using_grid" ) { Field spf = sp.createField( option::name( "spf" ) ); Field gpf = gp.createField( option::name( "gpf" ) ); + array::make_view( gpf ).assign(0); + EXPECT_NO_THROW( trans.dirtrans( gpf, spf ) ); EXPECT_NO_THROW( trans.invtrans( spf, gpf ) ); @@ -412,6 +416,8 @@ CASE( "test_trans_using_functionspace_NodeColumns" ) { Field spf = sp.createField( option::name( "spf" ) ); Field gpf = gp.createField( option::name( "gpf" ) ); + array::make_view( gpf ).assign(0); + EXPECT_NO_THROW( trans.dirtrans( gpf, spf ) ); EXPECT_NO_THROW( trans.invtrans( spf, gpf ) ); @@ -438,6 +444,8 @@ CASE( "test_trans_using_functionspace_StructuredColumns" ) { Field spf = sp.createField( option::name( "spf" ) ); Field gpf = gp.createField( option::name( "gpf" ) ); + array::make_view( gpf ).assign(0); + EXPECT_NO_THROW( trans.dirtrans( gpf, spf ) ); EXPECT_NO_THROW( trans.invtrans( spf, gpf ) ); @@ -460,7 +468,7 @@ CASE( "test_trans_MIR_lonlat" ) { trans::Trans trans( grid, 47 ); // global fields - std::vector spf( trans.spectralCoefficients() ); + std::vector spf( trans.spectralCoefficients(), 0. ); std::vector gpf( grid.size() ); if ( mpi::comm().size() == 1 ) { From 412563608c4c32b1c9af7f8bc8c55a9d88d22e0c Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 31 Jul 2018 15:32:22 +0000 Subject: [PATCH 12/21] Fix wrong mod from previous commit --- src/atlas/trans/ifs/TransIFS.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/atlas/trans/ifs/TransIFS.cc b/src/atlas/trans/ifs/TransIFS.cc index c1ec7c786..34c36d605 100644 --- a/src/atlas/trans/ifs/TransIFS.cc +++ b/src/atlas/trans/ifs/TransIFS.cc @@ -796,8 +796,6 @@ void TransIFS::__dirtrans( const functionspace::NodeColumns& gp, const FieldSet& transform.nscalar = nfld; transform.rgp = rgp.data(); transform.rspscalar = rspec.data(); - transform.ngpblks = ngptot(); - transform.nproma = 1; TRANS_CHECK(::trans_dirtrans( &transform ) ); } From 1f1d11ea805b0d79ffcf3954a96ff26698428473 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 7 Aug 2018 17:41:10 +0000 Subject: [PATCH 13/21] Running cuda 9.1 on lxg --- CMakeLists.txt | 2 +- src/tests/CMakeLists.txt | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index bd0231591..657d07dfe 100644 --- a/CMakeLists.txt +++ b/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 ) diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 8e6406ff4..b82fecc35 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -53,8 +53,11 @@ macro( atlas_add_cuda_test ) set_property( TARGET ${_PAR_TARGET} PROPERTY BUILD_WITH_INSTALL_RPATH FALSE ) set_property( TARGET ${_PAR_TARGET} PROPERTY SKIP_BUILD_RPATH FALSE ) - add_test (${_PAR_TARGET} ${_PAR_TARGET}) - + if( CMAKE_CROSSCOMPILING_EMULATOR ) + add_test ( NAME ${_PAR_TARGET} COMMAND ${CMAKE_CROSSCOMPILING_EMULATOR} ${_PAR_TARGET} ) + else() + add_test ( NAME ${_PAR_TARGET} COMMAND ${_PAR_TARGET} ) + endif() endif() endmacro() From a8ace6fa1033a207e723b686f40a821a13c15ad2 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Wed, 8 Aug 2018 15:35:13 +0100 Subject: [PATCH 14/21] Silenced warning --- src/atlas/grid/detail/grid/Grid.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/atlas/grid/detail/grid/Grid.h b/src/atlas/grid/detail/grid/Grid.h index 3b2011b28..8f15dc2d4 100644 --- a/src/atlas/grid/detail/grid/Grid.h +++ b/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 From 58b361feb4ece09818cab2564b96b4398c59ddff Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Thu, 9 Aug 2018 11:03:43 +0100 Subject: [PATCH 15/21] fixed missing equator (MIR-285) and issue at poles also updated the unit tests --- src/atlas/trans/local/TransLocal.cc | 35 ++-- src/tests/trans/test_transgeneral.cc | 229 +++++++++++++++++++++++++-- 2 files changed, 239 insertions(+), 25 deletions(-) diff --git a/src/atlas/trans/local/TransLocal.cc b/src/atlas/trans/local/TransLocal.cc index 6856c0a9a..40231dc6d 100644 --- a/src/atlas/trans/local/TransLocal.cc +++ b/src/atlas/trans/local/TransLocal.cc @@ -31,6 +31,12 @@ #include #endif +// move latitudes at the poles to the following latitude: +// (otherwise we would divide by zero when computing u,v from U,V) +double latPole = 89.9999999; +// (latPole=89.9999999 seems to produce the best accuracy. Moving it further away +// or closer to the pole both increase the errors!) + namespace atlas { namespace trans { @@ -351,7 +357,7 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma else { Log::debug() << "Grid has " << nlats << " latitudes. Global grid has " << nlatsGlobal_ << std::endl; } - if ( useGlobalLeg ) { nlatsLeg_ = nlatsGlobal_ / 2; } + if ( useGlobalLeg ) { nlatsLeg_ = ( nlatsGlobal_ + 1 ) / 2; } else { nlatsLeg_ = nlatsLegDomain_; nlatsLegReduced_ = nlatsLeg_; @@ -432,12 +438,18 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma std::vector lons( nlonsMax ); if ( nlatsNH_ >= nlatsSH_ || useGlobalLeg ) { for ( size_t j = 0; j < nlatsLeg_; ++j ) { - lats[j] = gsLeg.y( j ) * util::Constants::degreesToRadians(); + double lat = gsLeg.y( j ); + if ( lat > latPole ) { lat = latPole; } + if ( lat < -latPole ) { lat = -latPole; } + lats[j] = lat * util::Constants::degreesToRadians(); } } else { for ( size_t j = nlats - 1, idx = 0; idx < nlatsLeg_; --j, ++idx ) { - lats[idx] = -gsLeg.y( j ) * util::Constants::degreesToRadians(); + double lat = gsLeg.y( j ); + if ( lat > latPole ) { lat = latPole; } + if ( lat < -latPole ) { lat = -latPole; } + lats[idx] = -lat * util::Constants::degreesToRadians(); } } for ( size_t j = 0; j < nlonsMax; ++j ) { @@ -562,7 +574,7 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma // write.close(); // } } - // other FFT implementations should be added with #elif statements + // other FFT implementations should be added with #elif statements #else useFFT_ = false; // no FFT implemented => default to dgemm std::string file_path = TransParameters( config ).write_fft(); @@ -880,8 +892,8 @@ void TransLocal::invtrans_legendre( const int truncation, const int nlats, const } } } - } // namespace trans - } // namespace atlas + } + } } // -------------------------------------------------------------------------------------------------------------------- @@ -1287,15 +1299,20 @@ void TransLocal::invtrans_uv( const int truncation, const int nb_scalar_fields, { if ( nb_vordiv_fields > 0 ) { ATLAS_TRACE( "compute u,v from U,V" ); - std::vector coslats( nlats ); + std::vector coslatinvs( nlats ); for ( size_t j = 0; j < nlats; ++j ) { - coslats[j] = std::cos( g.y( j ) * util::Constants::degreesToRadians() ); + double lat = g.y( j ); + if ( lat > latPole ) { lat = latPole; } + if ( lat < -latPole ) { lat = -latPole; } + double coslat = std::cos( lat * util::Constants::degreesToRadians() ); + coslatinvs[j] = 1. / coslat; + //Log::info() << "lat=" << g.y( j ) << " coslat=" << coslat << std::endl; } int idx = 0; for ( int jfld = 0; jfld < 2 * nb_vordiv_fields && jfld < nb_fields; jfld++ ) { for ( int jlat = 0; jlat < g.ny(); jlat++ ) { for ( int jlon = 0; jlon < g.nx( jlat ); jlon++ ) { - gp_fields[idx] /= coslats[jlat]; + gp_fields[idx] *= coslatinvs[jlat]; idx++; } } diff --git a/src/tests/trans/test_transgeneral.cc b/src/tests/trans/test_transgeneral.cc index 73c3942a3..6bd1f9a11 100644 --- a/src/tests/trans/test_transgeneral.cc +++ b/src/tests/trans/test_transgeneral.cc @@ -276,8 +276,8 @@ double sphericalharmonics_analytic_point( void spectral_transform_grid_analytic( const size_t trc, // truncation (in) bool trcFT, // truncation for Fourier transformation (in) - const double n, // total wave number (implemented so far for n<4 - const double m, // zonal wave number (implemented so far for m<4, m rgp1_analytic( g1.size() ); std::vector rgp2_analytic( g2.size() ); + StructuredGrid gs1( g1 ); + StructuredGrid gs2( g2 ); + double latPole = 89.9999999; + bool gridHasPole = gs1.y( 0 ) > latPole || gs2.y( 0 ) > latPole; + double tolerance; int icase = 0; for ( int ivar_in = 0; ivar_in < 3; ivar_in++ ) { // vorticity, divergence, scalar for ( int ivar_out = 0; ivar_out < 3; ivar_out++ ) { // u, v, scalar int nb_fld = 1; - if ( ivar_out == 2 ) { + if ( ivar_out == 2 && !gridHasPole ) { tolerance = 1.e-13; nb_fld = nb_scalar; } @@ -749,14 +758,201 @@ CASE( "test_trans_domain" ) { //for ( int j = 0; j < g1.size(); j++ ) { // Log::info() << rgp1[pos * g1.size() + j] << " "; //}; + //Log::info() << std::endl << "analytic1:"; + //for ( int j = 0; j < g1.size(); j++ ) { + // Log::info() << rgp1_analytic[j] << " "; + //}; //Log::info() << std::endl << "rgp2:"; //for ( int j = 0; j < g2.size(); j++ ) { // Log::info() << rgp2[pos * g2.size() + j] << " "; //}; + //Log::info() << std::endl << "analytic2:"; + //for ( int j = 0; j < g2.size(); j++ ) { + // Log::info() << rgp2_analytic[j] << " "; + //}; + //Log::info() << std::endl; + rav1 += rms_gen1; + rav2 += rms_gen2; + if ( !( rms_gen1 < tolerance ) || !( rms_gen2 < tolerance ) ) { + Log::info() + << "Case " << icase << " ivar_in=" << ivar_in << " ivar_out=" << ivar_out + << " m=" << m << " n=" << n << " imag=" << imag << " k=" << k << std::endl; + ATLAS_DEBUG_VAR( rms_gen1 ); + ATLAS_DEBUG_VAR( rms_gen2 ); + ATLAS_DEBUG_VAR( tolerance ); + } + EXPECT( rms_gen1 < tolerance ); + EXPECT( rms_gen2 < tolerance ); + icase++; + auto end = std::chrono::system_clock::now(); // + std::chrono::duration elapsed_seconds = end - start; + std::time_t end_time = std::chrono::system_clock::to_time_t( end ); + std::string time_str = std::ctime( &end_time ); + //Log::info() << "case " << icase << ", elapsed time: " << elapsed_seconds.count() + // << "s. Now: " << time_str.substr( 0, time_str.length() - 1 ) << std::endl; + } + k++; + } + } + } + } + } + } + Log::info() << "Vordiv+scalar comparison with trans: all " << icase << " cases successfully passed!" << std::endl; + rav1 /= icase; + Log::info() << "average RMS error of transLocal1: " << rav1 << std::endl; + rav2 /= icase; + Log::info() << "average RMS error of transLocal2: " << rav2 << std::endl; +} +//----------------------------------------------------------------------------- +CASE( "test_trans_pole" ) { + Log::info() << "test_trans_pole" << std::endl; + // test transform at the pole and with LinearSpacing grids + // not using caching in this test because not useful for LinearSpacing grids + + std::ostream& out = Log::info(); + + //Domain testdomain = ZonalBandDomain( {-90., 90.} ); + //Domain testdomain = ZonalBandDomain( {-.5, .5} ); + //Domain testdomain = RectangularDomain( {0., 30.}, {-.05, .05} ); + //Domain testdomain1 = ZonalBandDomain( {-10., 5.} ); + Domain testdomain1 = RectangularDomain( {-5., 5.}, {-2.5, 0.} ); + //Domain testdomain1 = RectangularDomain( {-1., 1.}, {50., 55.} ); + Domain testdomain2 = RectangularDomain( {-5., 5.}, {-2.5, 0.} ); + // Grid: (Adjust the following line if the test takes too long!) + + Grid global_grid1( "L3" ); + Grid global_grid2( "L3" ); + //Grid g1( global_grid, testdomain1 ); + //Grid g2( global_grid, testdomain2 ); + Grid g1( global_grid1 ); + //Grid g2( global_grid2 ); + + bool fourierTrc1 = true; + bool fourierTrc2 = false; + // fourierTrc1, fourierTrc2: need to be false if no global grid can be constructed + // (like for grids created with LinearSpacing) + + using grid::StructuredGrid; + using LinearSpacing = grid::LinearSpacing; + StructuredGrid g2( LinearSpacing( {0., 180.}, 3 ), LinearSpacing( {89., 90.}, 2 ) ); + // when using LinearSpacing: set fourierTrc2 to false + + int trc = 2; + Trace t1( Here(), "translocal1 construction" ); + trans::Trans transLocal1( global_grid1, g1.domain(), trc, option::type( "local" ) ); + t1.stop(); + Trace t2( Here(), "translocal2 construction" ); + trans::Trans transLocal2( g2, trc, option::type( "local" ) ); + t2.stop(); + + double rav1 = 0., rav2 = 0.; // compute average rms errors of transLocal1 and transLocal2 + + functionspace::Spectral spectral( trc ); + + int nb_scalar = 1, nb_vordiv = 1; + int N = ( trc + 2 ) * ( trc + 1 ) / 2, nb_all = nb_scalar + 2 * nb_vordiv; + std::vector sp( 2 * N * nb_scalar ); + std::vector vor( 2 * N * nb_vordiv ); + std::vector div( 2 * N * nb_vordiv ); + std::vector rspecg( 2 * N ); + std::vector rgp1( nb_all * g1.size() ); + std::vector rgp2( nb_all * g2.size() ); + std::vector rgp1_analytic( g1.size() ); + std::vector rgp2_analytic( g2.size() ); + + StructuredGrid gs1( g1 ); + StructuredGrid gs2( g2 ); + double latPole = 89.9999999; + bool gridHasPole = gs1.y( 0 ) > latPole || gs2.y( 0 ) > latPole; + double tolerance; + int icase = 0; + for ( int ivar_in = 0; ivar_in < 3; ivar_in++ ) { // vorticity, divergence, scalar + for ( int ivar_out = 0; ivar_out < 3; ivar_out++ ) { // u, v, scalar + int nb_fld = 1; + if ( ivar_out == 2 && !gridHasPole ) { + tolerance = 1.e-13; + nb_fld = nb_scalar; + } + else { + tolerance = 2.e-6; + nb_fld = nb_vordiv; + } + for ( int jfld = 0; jfld < nb_fld; jfld++ ) { // multiple fields + int k = 0; + for ( int m = 0; m <= trc; m++ ) { // zonal wavenumber + for ( int n = m; n <= trc; n++ ) { // total wavenumber + for ( int imag = 0; imag <= 1; imag++ ) { // real and imaginary part + + if ( sphericalharmonics_analytic_point( n, m, true, 0., 0., ivar_in, ivar_in ) == 0. && + icase < 1000 ) { + auto start = std::chrono::system_clock::now(); + for ( int j = 0; j < 2 * N * nb_scalar; j++ ) { + sp[j] = 0.; + } + for ( int j = 0; j < 2 * N * nb_vordiv; j++ ) { + vor[j] = 0.; + div[j] = 0.; + } + if ( ivar_in == 0 ) vor[k * nb_vordiv + jfld] = 1.; + if ( ivar_in == 1 ) div[k * nb_vordiv + jfld] = 1.; + if ( ivar_in == 2 ) sp[k * nb_scalar + jfld] = 1.; + + for ( int j = 0; j < nb_all * g1.size(); j++ ) { + rgp1[j] = 0.; + } + for ( int j = 0; j < nb_all * g2.size(); j++ ) { + rgp2[j] = 0.; + } + for ( int j = 0; j < g1.size(); j++ ) { + rgp1_analytic[j] = 0.; + } + for ( int j = 0; j < g2.size(); j++ ) { + rgp2_analytic[j] = 0.; + } + + spectral_transform_grid_analytic( trc, fourierTrc1, n, m, imag, g1, rspecg.data(), + rgp1_analytic.data(), ivar_in, ivar_out ); + + spectral_transform_grid_analytic( trc, fourierTrc2, n, m, imag, g2, rspecg.data(), + rgp2_analytic.data(), ivar_in, ivar_out ); + + //Log::info() << std::endl << "rgp1:"; + ATLAS_TRACE_SCOPE( "translocal1" ) + EXPECT_NO_THROW( transLocal1.invtrans( nb_scalar, sp.data(), nb_vordiv, vor.data(), + div.data(), rgp1.data() ) ); + + //Log::info() << std::endl << "rgp2:"; + ATLAS_TRACE_SCOPE( "translocal2" ) + EXPECT_NO_THROW( transLocal2.invtrans( nb_scalar, sp.data(), nb_vordiv, vor.data(), + div.data(), rgp2.data() ) ); + + int pos = ( ivar_out * nb_vordiv + jfld ); + + double rms_gen1 = + compute_rms( g1.size(), rgp1.data() + pos * g1.size(), rgp1_analytic.data() ); + + double rms_gen2 = + compute_rms( g2.size(), rgp2.data() + pos * g2.size(), rgp2_analytic.data() ); + + //Log::info() << "Case " << icase << " ivar_in=" << ivar_in << " ivar_out=" << ivar_out + // << " m=" << m << " n=" << n << " imag=" << imag << " k=" << k << std::endl + // << "rgp1:"; + //for ( int j = 0; j < g1.size(); j++ ) { + // Log::info() << rgp1[pos * g1.size() + j] << " "; + //}; //Log::info() << std::endl << "analytic1:"; //for ( int j = 0; j < g1.size(); j++ ) { // Log::info() << rgp1_analytic[j] << " "; //}; + //Log::info() << std::endl << "rgp2:"; + //for ( int j = 0; j < g2.size(); j++ ) { + // Log::info() << rgp2[pos * g2.size() + j] << " "; + //}; + //Log::info() << std::endl << "analytic2:"; + //for ( int j = 0; j < g2.size(); j++ ) { + // Log::info() << rgp2_analytic[j] << " "; + //}; //Log::info() << std::endl; rav1 += rms_gen1; rav2 += rms_gen2; @@ -793,7 +989,6 @@ CASE( "test_trans_domain" ) { } #endif //----------------------------------------------------------------------------- -//----------------------------------------------------------------------------- #if 1 CASE( "test_trans_unstructured" ) { Log::info() << "test_trans_unstructured" << std::endl; @@ -805,9 +1000,9 @@ CASE( "test_trans_unstructured" ) { //Domain testdomain = RectangularDomain( {20., 25.}, {40., 60.} ); Domain testdomain = RectangularDomain( {0., 90.}, {0., 90.} ); // Grid: (Adjust the following line if the test takes too long!) - Grid grid_global( "F120" ); + Grid grid_global( "F32" ); Grid g( grid_global, testdomain ); - int trc = 120; + int trc = 31; grid::StructuredGrid gs( g ); std::vector pts( g.size() ); int idx( 0 ); @@ -843,8 +1038,8 @@ CASE( "test_trans_unstructured" ) { trans::Trans transLocal2( gu, trc, option::type( "local" ) | util::Config( "precompute", false ) ); int icase = 0; - for ( int ivar_in = 2; ivar_in < 3; ivar_in++ ) { // vorticity, divergence, scalar - for ( int ivar_out = 2; ivar_out < 3; ivar_out++ ) { // u, v, scalar + for ( int ivar_in = 0; ivar_in < 3; ivar_in++ ) { // vorticity, divergence, scalar + for ( int ivar_out = 0; ivar_out < 3; ivar_out++ ) { // u, v, scalar int nb_fld = 1; if ( ivar_out == 2 ) { tolerance = 1.e-13; @@ -861,7 +1056,7 @@ CASE( "test_trans_unstructured" ) { for ( int imag = 0; imag <= 1; imag++ ) { // real and imaginary part if ( sphericalharmonics_analytic_point( n, m, true, 0., 0., ivar_in, ivar_in ) == 0. && - icase < 100 ) { + icase < 1000 ) { auto start = std::chrono::system_clock::now(); for ( int j = 0; j < 2 * N * nb_scalar; j++ ) { sp[j] = 0.; @@ -887,17 +1082,19 @@ CASE( "test_trans_unstructured" ) { rgp_analytic2[j] = 0.; } - spectral_transform_grid_analytic( trc, trc, n, m, imag, g, rspecg.data(), + spectral_transform_grid_analytic( trc, false, n, m, imag, g, rspecg.data(), rgp_analytic1.data(), ivar_in, ivar_out ); //Log::info() << icase << " m=" << m << " n=" << n << " imag=" << imag << " structured: "; + ATLAS_TRACE_SCOPE( "structured" ) EXPECT_NO_THROW( transLocal1.invtrans( nb_scalar, sp.data(), nb_vordiv, vor.data(), div.data(), rgp1.data() ) ); - spectral_transform_grid_analytic( trc, trc, n, m, imag, gu, rspecg.data(), + spectral_transform_grid_analytic( trc, false, n, m, imag, gu, rspecg.data(), rgp_analytic2.data(), ivar_in, ivar_out ); //Log::info() << icase << " m=" << m << " n=" << n << " imag=" << imag << " unstructured: "; + ATLAS_TRACE_SCOPE( "unstructured" ) EXPECT_NO_THROW( transLocal2.invtrans( nb_scalar, sp.data(), nb_vordiv, vor.data(), div.data(), rgp2.data() ) ); @@ -944,7 +1141,7 @@ CASE( "test_trans_unstructured" ) { } #endif - //----------------------------------------------------------------------------- +//----------------------------------------------------------------------------- #if 0 CASE( "test_trans_fourier_truncation" ) { From 23cada92e03a175fe5d9157c149d2de22285c82f Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Thu, 9 Aug 2018 15:18:02 +0100 Subject: [PATCH 16/21] removed float comparison which might have caused issues at the equator --- src/atlas/trans/local/TransLocal.cc | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/atlas/trans/local/TransLocal.cc b/src/atlas/trans/local/TransLocal.cc index 40231dc6d..9269a5e4a 100644 --- a/src/atlas/trans/local/TransLocal.cc +++ b/src/atlas/trans/local/TransLocal.cc @@ -25,6 +25,7 @@ #include "eckit/linalg/Matrix.h" #include "eckit/log/Bytes.h" #include "eckit/parser/JSON.h" +#include "eckit/types/FloatCompare.h" #include "atlas/library/defines.h" #if ATLAS_HAVE_FFTW @@ -313,9 +314,9 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma // assumptions: latitudes in g.y(j) are monotone and decreasing // no assumption on whether we have 0, 1 or 2 latitudes at the equator double lat = g.y( j ); - if ( lat > 0. ) { nlatsNH_++; } - if ( lat == 0. ) { neqtr++; } - if ( lat < 0. ) { nlatsSH_++; } + if ( eckit::types::is_strictly_greater( lat, 0. ) ) { nlatsNH_++; } + if ( eckit::types::is_approximately_equal( lat, 0. ) ) { neqtr++; } + if ( eckit::types::is_strictly_greater( 0., lat ) ) { nlatsSH_++; } } if ( neqtr > 0 ) { nlatsNH_++; From 720034f10bf0683fa65877c7613a68a8cb996651 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 9 Aug 2018 15:02:05 +0000 Subject: [PATCH 17/21] Fix TransIFS with Gridtools-CUDA backend --- src/atlas/trans/ifs/TransIFS.cc | 324 ++++++++++++-------- src/tests/trans/fctest_trans.F90 | 2 + src/tests/trans/test_trans_invtrans_grad.cc | 1 + 3 files changed, 191 insertions(+), 136 deletions(-) diff --git a/src/atlas/trans/ifs/TransIFS.cc b/src/atlas/trans/ifs/TransIFS.cc index 34c36d605..791a47017 100644 --- a/src/atlas/trans/ifs/TransIFS.cc +++ b/src/atlas/trans/ifs/TransIFS.cc @@ -25,7 +25,9 @@ using Topology = atlas::mesh::Nodes::Topology; using atlas::Field; using atlas::FunctionSpace; using atlas::array::ArrayView; +using atlas::array::LocalView; using atlas::array::make_view; +using atlas::array::make_shape; using atlas::functionspace::NodeColumns; using atlas::functionspace::Spectral; using atlas::functionspace::StructuredColumns; @@ -105,6 +107,25 @@ void trans_check( const int code, const char* msg, const eckit::CodeLocation& lo } #define TRANS_CHECK( CALL ) trans_check( CALL, #CALL, Here() ) + +static int compute_nfld( const Field& f ) { + int nfld=1; + for( int i=1; i& rgpview_; + LocalView& rgpview_; IsGhostNode is_ghost; size_t f; - PackNodeColumns( ArrayView& rgpview, const NodeColumns& fs ) : + PackNodeColumns( LocalView& rgpview, const NodeColumns& fs ) : rgpview_( rgpview ), is_ghost( fs.nodes() ), f( 0 ) {} @@ -359,10 +380,10 @@ struct PackNodeColumns { }; struct PackStructuredColumns { - ArrayView& rgpview_; + LocalView& rgpview_; size_t f; - PackStructuredColumns( ArrayView& rgpview ) : rgpview_( rgpview ), f( 0 ) {} + PackStructuredColumns( LocalView& rgpview ) : rgpview_( rgpview ), f( 0 ) {} void operator()( const Field& field ) { switch ( field.rank() ) { @@ -403,9 +424,9 @@ struct PackStructuredColumns { }; struct PackSpectral { - ArrayView& rspecview_; + LocalView& rspecview_; size_t f; - PackSpectral( ArrayView& rspecview ) : rspecview_( rspecview ), f( 0 ) {} + PackSpectral( LocalView& rspecview ) : rspecview_( rspecview ), f( 0 ) {} void operator()( const Field& field ) { switch ( field.rank() ) { @@ -445,11 +466,11 @@ struct PackSpectral { }; struct UnpackNodeColumns { - const ArrayView& rgpview_; + const LocalView& rgpview_; IsGhostNode is_ghost; size_t f; - UnpackNodeColumns( const ArrayView& rgpview, const NodeColumns& fs ) : + UnpackNodeColumns( const LocalView& rgpview, const NodeColumns& fs ) : rgpview_( rgpview ), is_ghost( fs.nodes() ), f( 0 ) {} @@ -516,10 +537,10 @@ struct UnpackNodeColumns { }; struct UnpackStructuredColumns { - const ArrayView& rgpview_; + const LocalView& rgpview_; size_t f; - UnpackStructuredColumns( const ArrayView& rgpview ) : rgpview_( rgpview ), f( 0 ) {} + UnpackStructuredColumns( const LocalView& rgpview ) : rgpview_( rgpview ), f( 0 ) {} void operator()( Field& field ) { switch ( field.rank() ) { @@ -560,9 +581,9 @@ struct UnpackStructuredColumns { }; struct UnpackSpectral { - const ArrayView& rspecview_; + const LocalView& rspecview_; size_t f; - UnpackSpectral( const ArrayView& rspecview ) : rspecview_( rspecview ), f( 0 ) {} + UnpackSpectral( const LocalView& rspecview ) : rspecview_( rspecview ), f( 0 ) {} void operator()( Field& field ) { switch ( field.rank() ) { @@ -758,30 +779,22 @@ void TransIFS::__dirtrans( const functionspace::NodeColumns& gp, const Field& gp void TransIFS::__dirtrans( const functionspace::NodeColumns& gp, const FieldSet& gpfields, const Spectral& sp, FieldSet& spfields, const eckit::Configuration& ) const { + assertCompatibleDistributions( gp, sp ); // Count total number of fields and do sanity checks - int nfld( 0 ); - for ( size_t jfld = 0; jfld < gpfields.size(); ++jfld ) { - const Field& f = gpfields[jfld]; - nfld += f.stride( 0 ); - } + const int nfld = compute_nfld( gpfields ); + const int trans_sp_nfld = compute_nfld( spfields ); - int trans_spnfld( 0 ); - for ( size_t jfld = 0; jfld < spfields.size(); ++jfld ) { - const Field& f = spfields[jfld]; - trans_spnfld += f.stride( 0 ); - } - - if ( nfld != trans_spnfld ) { + if ( nfld != trans_sp_nfld ) { throw eckit::SeriousBug( "dirtrans: different number of gridpoint fields than spectral fields", Here() ); } - // Arrays Trans expects - array::ArrayT rgp( nfld, ngptot() ); - array::ArrayT rspec( nspec2(), nfld ); - array::ArrayView rgpview = array::make_view( rgp ); - array::ArrayView rspecview = array::make_view( rspec ); + // Arrays Trans expects + std::vector rgp( nfld * ngptot() ); + std::vector rsp( nspec2() * nfld ); + auto rgpview = LocalView( rgp.data(), make_shape( nfld, ngptot() ) ); + auto rspview = LocalView( rsp.data(), make_shape( nspec2(), nfld ) ); // Pack gridpoints { @@ -794,14 +807,14 @@ void TransIFS::__dirtrans( const functionspace::NodeColumns& gp, const FieldSet& { struct ::DirTrans_t transform = ::new_dirtrans( trans_.get() ); transform.nscalar = nfld; - transform.rgp = rgp.data(); - transform.rspscalar = rspec.data(); + transform.rgp = rgp.data(); + transform.rspscalar = rsp.data(); TRANS_CHECK(::trans_dirtrans( &transform ) ); } // Unpack the spectral fields { - UnpackSpectral unpack( rspecview ); + UnpackSpectral unpack( rspview ); for ( size_t jfld = 0; jfld < spfields.size(); ++jfld ) unpack( spfields[jfld] ); } @@ -811,58 +824,73 @@ void TransIFS::__dirtrans( const functionspace::NodeColumns& gp, const FieldSet& void TransIFS::__dirtrans( const StructuredColumns& gp, const Field& gpfield, const Spectral& sp, Field& spfield, const eckit::Configuration& ) const { + ASSERT( gpfield.functionspace() == 0 || functionspace::StructuredColumns( gpfield.functionspace() ) ); ASSERT( spfield.functionspace() == 0 || functionspace::Spectral( spfield.functionspace() ) ); assertCompatibleDistributions( gp, sp ); - if ( gpfield.stride( 0 ) != spfield.stride( 0 ) ) { + if ( compute_nfld( gpfield ) != compute_nfld( spfield ) ) { throw eckit::SeriousBug( "dirtrans: different number of gridpoint fields than spectral fields", Here() ); } if ( (int)gpfield.shape( 0 ) != ngptot() ) { throw eckit::SeriousBug( "dirtrans: slowest moving index must be ngptot", Here() ); } - const int nfld = gpfield.stride( 0 ); + const int nfld = compute_nfld( gpfield ); + + // Arrays Trans expects + std::vector rgp( nfld * ngptot() ); + std::vector rsp( nspec2() * nfld ); + auto rgpview = LocalView( rgp.data(), make_shape( nfld, ngptot() ) ); + auto rspview = LocalView( rsp.data(), make_shape( nspec2(), nfld ) ); + + // Pack gridpoints + { + PackStructuredColumns pack( rgpview ); + pack( gpfield ); + } // Do transform { struct ::DirTrans_t transform = ::new_dirtrans( trans_.get() ); transform.nscalar = nfld; - transform.rgp = gpfield.data(); - transform.rspscalar = spfield.data(); - transform.ngpblks = gpfield.shape( 0 ); + transform.rgp = rgp.data(); + transform.rspscalar = rsp.data(); + transform.ngpblks = ngptot(); transform.nproma = 1; TRANS_CHECK(::trans_dirtrans( &transform ) ); } + + // Unpack spectral + { + UnpackSpectral unpack( rspview ); + unpack( spfield ); + } + } void TransIFS::__dirtrans( const StructuredColumns& gp, const FieldSet& gpfields, const Spectral& sp, FieldSet& spfields, const eckit::Configuration& ) const { + assertCompatibleDistributions( gp, sp ); // Count total number of fields and do sanity checks - int nfld( 0 ); + const int nfld = compute_nfld( gpfields ); for ( size_t jfld = 0; jfld < gpfields.size(); ++jfld ) { const Field& f = gpfields[jfld]; - nfld += f.stride( 0 ); ASSERT( f.functionspace() == 0 || functionspace::StructuredColumns( f.functionspace() ) ); } - int trans_spnfld( 0 ); - for ( size_t jfld = 0; jfld < spfields.size(); ++jfld ) { - const Field& f = spfields[jfld]; - trans_spnfld += f.stride( 0 ); - } + const int trans_sp_nfld = compute_nfld( spfields ); - if ( nfld != trans_spnfld ) { + if ( nfld != trans_sp_nfld ) { throw eckit::SeriousBug( "dirtrans: different number of gridpoint fields than spectral fields", Here() ); } // Arrays Trans expects - array::ArrayT rgp( nfld, ngptot() ); - array::ArrayT rspec( nspec2(), nfld ); - - array::ArrayView rgpview = array::make_view( rgp ); - array::ArrayView rspecview = array::make_view( rspec ); + std::vector rgp( nfld * ngptot() ); + std::vector rsp( nspec2() * nfld ); + auto rgpview = LocalView( rgp.data(), make_shape( nfld, ngptot() ) ); + auto rspview = LocalView( rsp.data(), make_shape( nspec2(), nfld ) ); // Pack gridpoints { @@ -875,15 +903,15 @@ void TransIFS::__dirtrans( const StructuredColumns& gp, const FieldSet& gpfields { struct ::DirTrans_t transform = ::new_dirtrans( trans_.get() ); transform.nscalar = nfld; - transform.rgp = rgp.data(); - transform.rspscalar = rspec.data(); + transform.rgp = rgp.data(); + transform.rspscalar = rsp.data(); TRANS_CHECK(::trans_dirtrans( &transform ) ); } // Unpack the spectral fields { - UnpackSpectral unpack( rspecview ); + UnpackSpectral unpack( rspview ); for ( size_t jfld = 0; jfld < spfields.size(); ++jfld ) unpack( spfields[jfld] ); } @@ -902,21 +930,12 @@ void TransIFS::__invtrans_grad( const Spectral& sp, const Field& spfield, const void TransIFS::__invtrans_grad( const Spectral& sp, const FieldSet& spfields, const functionspace::NodeColumns& gp, FieldSet& gradfields, const eckit::Configuration& config ) const { + assertCompatibleDistributions( gp, sp ); // Count total number of fields and do sanity checks - int nb_gridpoint_field( 0 ); - for ( size_t jfld = 0; jfld < gradfields.size(); ++jfld ) { - const Field& f = gradfields[jfld]; - nb_gridpoint_field += f.stride( 0 ); - } - - int nfld( 0 ); - for ( size_t jfld = 0; jfld < spfields.size(); ++jfld ) { - const Field& f = spfields[jfld]; - nfld += f.stride( 0 ); - ASSERT( std::max( 1, f.levels() ) == f.stride( 0 ) ); - } + const int nb_gridpoint_field = compute_nfld( gradfields ); + const int nfld = compute_nfld( spfields ); if ( nb_gridpoint_field != 2 * nfld ) // factor 2 because N-S and E-W derivatives throw eckit::SeriousBug( @@ -926,16 +945,15 @@ void TransIFS::__invtrans_grad( const Spectral& sp, const FieldSet& spfields, co // Arrays Trans expects // Allocate space for - array::ArrayT rgp( 3 * nfld, - ngptot() ); // (scalars) + (NS ders) + (EW ders) - array::ArrayT rspec( nspec2(), nfld ); + std::vector rgp( 3 * nfld * ngptot() ); // (scalars) + (NS ders) + (EW ders) + std::vector rsp( nspec2() * nfld ); + auto rgpview = LocalView( rgp.data(), make_shape(3 * nfld, ngptot() ) ); + auto rspview = LocalView( rsp.data(), make_shape( nspec2(), nfld ) ); - array::ArrayView rgpview = array::make_view( rgp ); - array::ArrayView rspecview = array::make_view( rspec ); // Pack spectral fields { - PackSpectral pack( rspecview ); + PackSpectral pack( rspview ); for ( size_t jfld = 0; jfld < spfields.size(); ++jfld ) pack( spfields[jfld] ); } @@ -944,8 +962,8 @@ void TransIFS::__invtrans_grad( const Spectral& sp, const FieldSet& spfields, co { struct ::InvTrans_t transform = ::new_invtrans( trans_.get() ); transform.nscalar = nfld; - transform.rgp = rgp.data(); - transform.rspscalar = rspec.data(); + transform.rgp = rgp.data(); + transform.rspscalar = rsp.data(); transform.lscalarders = true; TRANS_CHECK(::trans_invtrans( &transform ) ); @@ -957,23 +975,33 @@ void TransIFS::__invtrans_grad( const Spectral& sp, const FieldSet& spfields, co int f = nfld; // skip to where derivatives start for ( size_t dim = 0; dim < 2; ++dim ) { for ( size_t jfld = 0; jfld < gradfields.size(); ++jfld ) { - const size_t nlev = std::max( 1, gradfields[jfld].levels() ); const size_t nb_nodes = gradfields[jfld].shape( 0 ); - - array::LocalView field( gradfields[jfld].data(), - array::make_shape( nb_nodes, nlev, 2 ) ); - - for ( size_t jlev = 0; jlev < nlev; ++jlev ) { + const size_t nlev = gradfields[jfld].levels(); + if( nlev ) { + auto field = make_view( gradfields[jfld] ); + for ( size_t jlev = 0; jlev < nlev; ++jlev ) { + int n = 0; + for ( size_t jnode = 0; jnode < nb_nodes; ++jnode ) { + if ( !is_ghost( jnode ) ) { + field( jnode, jlev, 1 - dim ) = rgpview( f, n ); + ++n; + } + } + ASSERT( n == ngptot() ); + } + } + else { + auto field = make_view( gradfields[jfld] ); int n = 0; for ( size_t jnode = 0; jnode < nb_nodes; ++jnode ) { if ( !is_ghost( jnode ) ) { - field( jnode, jlev, 1 - dim ) = rgpview( f, n ); + field( jnode, 1 - dim ) = rgpview( f, n ); ++n; } } ASSERT( n == ngptot() ); - ++f; } + ++f; } } } @@ -994,34 +1022,26 @@ void TransIFS::__invtrans( const Spectral& sp, const Field& spfield, const funct void TransIFS::__invtrans( const Spectral& sp, const FieldSet& spfields, const functionspace::NodeColumns& gp, FieldSet& gpfields, const eckit::Configuration& config ) const { + assertCompatibleDistributions( gp, sp ); - // Count total number of fields and do sanity checks - int nfld( 0 ); - for ( size_t jfld = 0; jfld < gpfields.size(); ++jfld ) { - const Field& f = gpfields[jfld]; - nfld += f.stride( 0 ); - } - int nb_spectral_fields( 0 ); - for ( size_t jfld = 0; jfld < spfields.size(); ++jfld ) { - const Field& f = spfields[jfld]; - nb_spectral_fields += f.stride( 0 ); - } + // Count total number of fields and do sanity checks + const int nfld = compute_nfld( gpfields ); + const int nb_spectral_fields = compute_nfld( spfields ); if ( nfld != nb_spectral_fields ) throw eckit::SeriousBug( "invtrans: different number of gridpoint fields than spectral fields", Here() ); // Arrays Trans expects - array::ArrayT rgp( nfld, ngptot() ); - array::ArrayT rspec( nspec2(), nfld ); - - array::ArrayView rgpview = array::make_view( rgp ); - array::ArrayView rspecview = array::make_view( rspec ); + std::vector rgp( nfld * ngptot() ); + std::vector rsp( nspec2() * nfld ); + auto rgpview = LocalView( rgp.data(), make_shape( nfld, ngptot() ) ); + auto rspview = LocalView( rsp.data(), make_shape( nspec2(), nfld ) ); // Pack spectral fields { - PackSpectral pack( rspecview ); + PackSpectral pack( rspview ); for ( size_t jfld = 0; jfld < spfields.size(); ++jfld ) pack( spfields[jfld] ); } @@ -1030,9 +1050,8 @@ void TransIFS::__invtrans( const Spectral& sp, const FieldSet& spfields, const f { struct ::InvTrans_t transform = ::new_invtrans( trans_.get() ); transform.nscalar = nfld; - transform.rgp = rgp.data(); - transform.rspscalar = rspec.data(); - + transform.rgp = rgp.data(); + transform.rspscalar = rsp.data(); TRANS_CHECK(::trans_invtrans( &transform ) ); } @@ -1049,28 +1068,47 @@ void TransIFS::__invtrans( const Spectral& sp, const FieldSet& spfields, const f void TransIFS::__invtrans( const functionspace::Spectral& sp, const Field& spfield, const functionspace::StructuredColumns& gp, Field& gpfield, const eckit::Configuration& config ) const { + assertCompatibleDistributions( gp, sp ); ASSERT( gpfield.functionspace() == 0 || functionspace::StructuredColumns( gpfield.functionspace() ) ); ASSERT( spfield.functionspace() == 0 || functionspace::Spectral( spfield.functionspace() ) ); - if ( gpfield.stride( 0 ) != spfield.stride( 0 ) ) { + if ( compute_nfld( gpfield ) != compute_nfld( spfield ) ) { throw eckit::SeriousBug( "dirtrans: different number of gridpoint fields than spectral fields", Here() ); } if ( (int)gpfield.shape( 0 ) != ngptot() ) { throw eckit::SeriousBug( "dirtrans: slowest moving index must be ngptot", Here() ); } - const int nfld = gpfield.stride( 0 ); + const int nfld = compute_nfld( gpfield ); + + // Arrays Trans expects + std::vector rgp( nfld * ngptot() ); + std::vector rsp( nspec2() * nfld ); + auto rgpview = LocalView( rgp.data(), make_shape( nfld, ngptot() ) ); + auto rspview = LocalView( rsp.data(), make_shape( nspec2(), nfld ) ); + + // Pack spectral fields + { + PackSpectral pack( rspview ); + pack( spfield ); + } // Do transform { struct ::InvTrans_t transform = ::new_invtrans( trans_.get() ); transform.nscalar = nfld; - transform.rgp = gpfield.data(); - transform.rspscalar = spfield.data(); - transform.ngpblks = gpfield.shape( 0 ); + transform.rgp = rgp.data(); + transform.rspscalar = rsp.data(); + transform.ngpblks = ngptot(); transform.nproma = 1; TRANS_CHECK(::trans_invtrans( &transform ) ); } + + // Unpack gridpoint fields + { + UnpackStructuredColumns unpack( rgpview ); + unpack( gpfield ); + } } // -------------------------------------------------------------------------------------------- @@ -1078,21 +1116,17 @@ void TransIFS::__invtrans( const functionspace::Spectral& sp, const Field& spfie void TransIFS::__invtrans( const functionspace::Spectral& sp, const FieldSet& spfields, const functionspace::StructuredColumns& gp, FieldSet& gpfields, const eckit::Configuration& config ) const { + assertCompatibleDistributions( gp, sp ); // Count total number of fields and do sanity checks - int nfld( 0 ); + const int nfld = compute_nfld( gpfields ); for ( size_t jfld = 0; jfld < gpfields.size(); ++jfld ) { const Field& f = gpfields[jfld]; - nfld += f.stride( 0 ); ASSERT( f.functionspace() == 0 || functionspace::StructuredColumns( f.functionspace() ) ); } - int nb_spectral_fields( 0 ); - for ( size_t jfld = 0; jfld < spfields.size(); ++jfld ) { - const Field& f = spfields[jfld]; - nb_spectral_fields += f.stride( 0 ); - } + const int nb_spectral_fields = compute_nfld( spfields ); if ( nfld != nb_spectral_fields ) { std::stringstream msg; @@ -1102,15 +1136,14 @@ void TransIFS::__invtrans( const functionspace::Spectral& sp, const FieldSet& sp } // Arrays Trans expects - array::ArrayT rgp( nfld, ngptot() ); - array::ArrayT rspec( nspec2(), nfld ); - - array::ArrayView rgpview = array::make_view( rgp ); - array::ArrayView rspecview = array::make_view( rspec ); + std::vector rgp( nfld * ngptot() ); + std::vector rsp( nspec2() * nfld ); + auto rgpview = LocalView( rgp.data(), make_shape( nfld, ngptot() ) ); + auto rspview = LocalView( rsp.data(), make_shape( nspec2(), nfld ) ); // Pack spectral fields { - PackSpectral pack( rspecview ); + PackSpectral pack( rspview ); for ( size_t jfld = 0; jfld < spfields.size(); ++jfld ) pack( spfields[jfld] ); } @@ -1119,8 +1152,8 @@ void TransIFS::__invtrans( const functionspace::Spectral& sp, const FieldSet& sp { struct ::InvTrans_t transform = ::new_invtrans( trans_.get() ); transform.nscalar = nfld; - transform.rgp = rgp.data(); - transform.rspscalar = rspec.data(); + transform.rgp = rgp.data(); + transform.rspscalar = rsp.data(); TRANS_CHECK(::trans_invtrans( &transform ) ); } @@ -1137,15 +1170,16 @@ void TransIFS::__invtrans( const functionspace::Spectral& sp, const FieldSet& sp void TransIFS::__dirtrans_wind2vordiv( const functionspace::NodeColumns& gp, const Field& gpwind, const Spectral& sp, Field& spvor, Field& spdiv, const eckit::Configuration& ) const { + assertCompatibleDistributions( gp, sp ); // Count total number of fields and do sanity checks - size_t nfld = spvor.stride( 0 ); + const size_t nfld = compute_nfld( spvor ); if ( spdiv.shape( 0 ) != spvor.shape( 0 ) ) throw eckit::SeriousBug( "invtrans: vorticity not compatible with divergence.", Here() ); if ( spdiv.shape( 1 ) != spvor.shape( 1 ) ) throw eckit::SeriousBug( "invtrans: vorticity not compatible with divergence.", Here() ); - size_t nwindfld = gpwind.stride( 0 ); + const size_t nwindfld = compute_nfld( gpwind); if ( nwindfld != 2 * nfld && nwindfld != 3 * nfld ) throw eckit::SeriousBug( "dirtrans: wind field is not compatible with vorticity, divergence.", Here() ); @@ -1161,8 +1195,12 @@ void TransIFS::__dirtrans_wind2vordiv( const functionspace::NodeColumns& gp, con if ( spdiv.size() == 0 ) throw eckit::SeriousBug( "dirtrans: spectral divergence field is empty." ); // Arrays Trans expects - array::ArrayT rgp( 2 * nfld, size_t( ngptot() ) ); - array::ArrayView rgpview = array::make_view( rgp ); + std::vector rgp( 2 * nfld * ngptot() ); + std::vector rspvor( nspec2() * nfld ); + std::vector rspdiv( nspec2() * nfld ); + auto rgpview = LocalView( rgp.data(), make_shape( 2 * nfld, ngptot() ) ); + auto rspvorview = LocalView( rspvor.data(), make_shape( nspec2(), nfld ) ); + auto rspdivview = LocalView( rspdiv.data(), make_shape( nspec2(), nfld ) ); // Pack gridpoints { @@ -1175,28 +1213,34 @@ void TransIFS::__dirtrans_wind2vordiv( const functionspace::NodeColumns& gp, con { struct ::DirTrans_t transform = ::new_dirtrans( trans_.get() ); transform.nvordiv = nfld; - transform.rgp = rgp.data(); - transform.rspvor = spvor.data(); - transform.rspdiv = spdiv.data(); + transform.rgp = rgp.data(); + transform.rspvor = rspvor.data(); + transform.rspdiv = rspdiv.data(); ASSERT( transform.rspvor ); ASSERT( transform.rspdiv ); TRANS_CHECK(::trans_dirtrans( &transform ) ); } + + // Pack spectral fields + UnpackSpectral unpack_vor( rspvorview ); unpack_vor( spvor ); + UnpackSpectral unpack_div( rspdivview ); unpack_div( spdiv ); + } void TransIFS::__invtrans_vordiv2wind( const Spectral& sp, const Field& spvor, const Field& spdiv, const functionspace::NodeColumns& gp, Field& gpwind, const eckit::Configuration& config ) const { + assertCompatibleDistributions( gp, sp ); // Count total number of fields and do sanity checks - size_t nfld = spvor.stride( 0 ); + const int nfld = compute_nfld( spvor ); if ( spdiv.shape( 0 ) != spvor.shape( 0 ) ) throw eckit::SeriousBug( "invtrans: vorticity not compatible with divergence.", Here() ); if ( spdiv.shape( 1 ) != spvor.shape( 1 ) ) throw eckit::SeriousBug( "invtrans: vorticity not compatible with divergence.", Here() ); - size_t nwindfld = gpwind.stride( 0 ); + const int nwindfld = compute_nfld( gpwind ); if ( nwindfld != 2 * nfld && nwindfld != 3 * nfld ) throw eckit::SeriousBug( "invtrans: wind field is not compatible with vorticity, divergence.", Here() ); @@ -1214,16 +1258,24 @@ void TransIFS::__invtrans_vordiv2wind( const Spectral& sp, const Field& spvor, c if ( spdiv.size() == 0 ) throw eckit::SeriousBug( "invtrans: spectral divergence field is empty." ); // Arrays Trans expects - array::ArrayT rgp( 2 * nfld, size_t( ngptot() ) ); - array::ArrayView rgpview = array::make_view( rgp ); + std::vector rgp( 2 * nfld * ngptot() ); + std::vector rspvor( nspec2() * nfld ); + std::vector rspdiv( nspec2() * nfld ); + auto rgpview = LocalView( rgp.data(), make_shape( 2 * nfld, ngptot() ) ); + auto rspvorview = LocalView( rspvor.data(), make_shape( nspec2(), nfld ) ); + auto rspdivview = LocalView( rspdiv.data(), make_shape( nspec2(), nfld ) ); + + // Pack spectral fields + PackSpectral pack_vor( rspvorview ); pack_vor( spvor ); + PackSpectral pack_div( rspdivview ); pack_div( spdiv ); // Do transform { struct ::InvTrans_t transform = ::new_invtrans( trans_.get() ); transform.nvordiv = nfld; - transform.rgp = rgp.data(); - transform.rspvor = spvor.data(); - transform.rspdiv = spdiv.data(); + transform.rgp = rgp.data(); + transform.rspvor = rspvor.data(); + transform.rspdiv = rspdiv.data(); ASSERT( transform.rspvor ); ASSERT( transform.rspdiv ); @@ -1397,7 +1449,7 @@ int atlas__Trans__truncation( const TransIFS* This ) { } const Grid::Implementation* atlas__Trans__grid( const TransIFS* This ) { - ATLAS_ERROR_HANDLING( ASSERT( This ); ASSERT( This->grid() ); ATLAS_DEBUG_VAR( This->grid().get()->owners() ); + ATLAS_ERROR_HANDLING( ASSERT( This ); ASSERT( This->grid() ); return This->grid().get(); ); return nullptr; } diff --git a/src/tests/trans/fctest_trans.F90 b/src/tests/trans/fctest_trans.F90 index d5b118c20..3c17d3649 100644 --- a/src/tests/trans/fctest_trans.F90 +++ b/src/tests/trans/fctest_trans.F90 @@ -13,6 +13,7 @@ module fctest_atlas_trans_fixture use atlas_module +use fckit_module use iso_c_binding implicit none character(len=1024) :: msg @@ -25,6 +26,7 @@ end module fctest_atlas_trans_fixture ! ----------------------------------------------------------------------------- TESTSUITE_INIT + call fckit_main%init() call atlas_library%initialise() END_TESTSUITE_INIT diff --git a/src/tests/trans/test_trans_invtrans_grad.cc b/src/tests/trans/test_trans_invtrans_grad.cc index f20376769..ac5bff210 100644 --- a/src/tests/trans/test_trans_invtrans_grad.cc +++ b/src/tests/trans/test_trans_invtrans_grad.cc @@ -27,6 +27,7 @@ #include "atlas/trans/Trans.h" #include "atlas/util/CoordinateEnums.h" #include "atlas/util/Earth.h" +#include "atlas/array.h" #include "tests/AtlasTestEnvironment.h" From 23e2e52c529a28204841155e1bd9d02e62e242ac Mon Sep 17 00:00:00 2001 From: Andreas Mueller Date: Mon, 20 Aug 2018 15:08:55 +0100 Subject: [PATCH 18/21] fixing L-grid domains on southern hemisphere (MIR-283) --- src/atlas/trans/local/TransLocal.cc | 9 +- src/tests/trans/test_transgeneral.cc | 188 +++++++++++++++++++++++++++ 2 files changed, 193 insertions(+), 4 deletions(-) diff --git a/src/atlas/trans/local/TransLocal.cc b/src/atlas/trans/local/TransLocal.cc index 9269a5e4a..962e10db4 100644 --- a/src/atlas/trans/local/TransLocal.cc +++ b/src/atlas/trans/local/TransLocal.cc @@ -314,9 +314,7 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma // assumptions: latitudes in g.y(j) are monotone and decreasing // no assumption on whether we have 0, 1 or 2 latitudes at the equator double lat = g.y( j ); - if ( eckit::types::is_strictly_greater( lat, 0. ) ) { nlatsNH_++; } - if ( eckit::types::is_approximately_equal( lat, 0. ) ) { neqtr++; } - if ( eckit::types::is_strictly_greater( 0., lat ) ) { nlatsSH_++; } + ( eckit::types::is_approximately_equal( lat, 0. ) ? neqtr : ( lat < 0 ? nlatsSH_ : nlatsNH_ ) )++; } if ( neqtr > 0 ) { nlatsNH_++; @@ -372,7 +370,10 @@ TransLocal::TransLocal( const Cache& cache, const Grid& grid, const Domain& doma //Log::info() << std::endl; int jlatMinLeg_ = jlatMin_; if ( nlatsNH_ < nlatsSH_ ) { jlatMinLeg_ += nlatsNH_ - nlatsSH_; }; - if ( jlatMin_ > nlatsGlobal_ / 2 ) { jlatMinLeg_ -= 2 * ( jlatMin_ - nlatsGlobal_ / 2 ); }; + if ( jlatMin_ >= ( nlatsGlobal_ + 1 ) / 2 ) { + jlatMinLeg_ -= 2 * ( jlatMin_ - ( nlatsGlobal_ + 1 ) / 2 ); + if ( nlatsGlobal_ % 2 == 1 ) { jlatMinLeg_--; } + }; if ( useGlobalLeg ) { nlatsLegReduced_ = jlatMinLeg_ + nlatsLegDomain_; } // reduce truncation towards the pole for reduced meshes: diff --git a/src/tests/trans/test_transgeneral.cc b/src/tests/trans/test_transgeneral.cc index 6bd1f9a11..5af4fafcd 100644 --- a/src/tests/trans/test_transgeneral.cc +++ b/src/tests/trans/test_transgeneral.cc @@ -990,6 +990,194 @@ CASE( "test_trans_pole" ) { #endif //----------------------------------------------------------------------------- #if 1 +CASE( "test_trans_southpole" ) { + Log::info() << "test_trans_southpole" << std::endl; + // test created for MIR-283 (limited area domain on the southern hemisphere with L-grid) + + std::ostream& out = Log::info(); + + //Domain testdomain = ZonalBandDomain( {-90., 90.} ); + //Domain testdomain = ZonalBandDomain( {-.5, .5} ); + //Domain testdomain = RectangularDomain( {0., 30.}, {-.05, .05} ); + //Domain testdomain1 = ZonalBandDomain( {-10., 5.} ); + //Domain testdomain1 = RectangularDomain( {-5., 5.}, {-2.5, 0.} ); + Domain testdomain1 = RectangularDomain( {0., 10.}, {-90., -10.} ); + //Domain testdomain1 = RectangularDomain( {-1., 1.}, {50., 55.} ); + //Domain testdomain2 = RectangularDomain( {-5., 5.}, {-2.5, 0.} ); + Domain testdomain2 = RectangularDomain( {0., 10.}, {10., 90.} ); + // Grid: (Adjust the following line if the test takes too long!) + + Grid global_grid1( "L9" ); + Grid global_grid2( "L9" ); + Grid g1( global_grid1, testdomain1 ); + Grid g2( global_grid2, testdomain2 ); + //Grid g1( global_grid1 ); + //Grid g2( global_grid2 ); + + bool fourierTrc1 = true; + bool fourierTrc2 = true; + // fourierTrc1, fourierTrc2: need to be false if no global grid can be constructed + // (like for grids created with LinearSpacing) + + using grid::StructuredGrid; + using LinearSpacing = grid::LinearSpacing; + //StructuredGrid g2( LinearSpacing( {0., 10.}, 2 ), LinearSpacing( {-10., -90.}, 9 ) ); + // when using LinearSpacing: set fourierTrc2 to false + + int trc = 8; + Trace t1( Here(), "translocal1 construction" ); + trans::Trans transLocal1( global_grid1, g1.domain(), trc, option::type( "local" ) ); + t1.stop(); + Trace t2( Here(), "translocal2 construction" ); + //trans::Trans transLocal2( g2, trc, option::type( "local" ) ); + trans::Trans transLocal2( global_grid2, g2.domain(), trc, option::type( "local" ) ); + t2.stop(); + + double rav1 = 0., rav2 = 0.; // compute average rms errors of transLocal1 and transLocal2 + + functionspace::Spectral spectral( trc ); + + int nb_scalar = 1, nb_vordiv = 1; + int N = ( trc + 2 ) * ( trc + 1 ) / 2, nb_all = nb_scalar + 2 * nb_vordiv; + std::vector sp( 2 * N * nb_scalar ); + std::vector vor( 2 * N * nb_vordiv ); + std::vector div( 2 * N * nb_vordiv ); + std::vector rspecg( 2 * N ); + std::vector rgp1( nb_all * g1.size() ); + std::vector rgp2( nb_all * g2.size() ); + std::vector rgp1_analytic( g1.size() ); + std::vector rgp2_analytic( g2.size() ); + + StructuredGrid gs1( g1 ); + StructuredGrid gs2( g2 ); + double latPole = 89.9999999; + bool gridHasPole = gs1.y( 0 ) > latPole || gs2.y( 0 ) > latPole; + double tolerance; + int icase = 0; + for ( int ivar_in = 0; ivar_in < 3; ivar_in++ ) { // vorticity, divergence, scalar + for ( int ivar_out = 0; ivar_out < 3; ivar_out++ ) { // u, v, scalar + int nb_fld = 1; + if ( ivar_out == 2 && !gridHasPole ) { + tolerance = 1.e-13; + nb_fld = nb_scalar; + } + else { + tolerance = 2.e-6; + nb_fld = nb_vordiv; + } + for ( int jfld = 0; jfld < nb_fld; jfld++ ) { // multiple fields + int k = 0; + for ( int m = 0; m <= trc; m++ ) { // zonal wavenumber + for ( int n = m; n <= trc; n++ ) { // total wavenumber + for ( int imag = 0; imag <= 1; imag++ ) { // real and imaginary part + + if ( sphericalharmonics_analytic_point( n, m, true, 0., 0., ivar_in, ivar_in ) == 0. && + icase < 1000 ) { + auto start = std::chrono::system_clock::now(); + for ( int j = 0; j < 2 * N * nb_scalar; j++ ) { + sp[j] = 0.; + } + for ( int j = 0; j < 2 * N * nb_vordiv; j++ ) { + vor[j] = 0.; + div[j] = 0.; + } + if ( ivar_in == 0 ) vor[k * nb_vordiv + jfld] = 1.; + if ( ivar_in == 1 ) div[k * nb_vordiv + jfld] = 1.; + if ( ivar_in == 2 ) sp[k * nb_scalar + jfld] = 1.; + + for ( int j = 0; j < nb_all * g1.size(); j++ ) { + rgp1[j] = 0.; + } + for ( int j = 0; j < nb_all * g2.size(); j++ ) { + rgp2[j] = 0.; + } + for ( int j = 0; j < g1.size(); j++ ) { + rgp1_analytic[j] = 0.; + } + for ( int j = 0; j < g2.size(); j++ ) { + rgp2_analytic[j] = 0.; + } + + spectral_transform_grid_analytic( trc, fourierTrc1, n, m, imag, g1, rspecg.data(), + rgp1_analytic.data(), ivar_in, ivar_out ); + + spectral_transform_grid_analytic( trc, fourierTrc2, n, m, imag, g2, rspecg.data(), + rgp2_analytic.data(), ivar_in, ivar_out ); + + //Log::info() << "Case " << icase << " ivar_in=" << ivar_in << " ivar_out=" << ivar_out + // << " m=" << m << " n=" << n << " imag=" << imag << " k=" << k << std::endl; + + ATLAS_TRACE_SCOPE( "translocal1" ) + EXPECT_NO_THROW( transLocal1.invtrans( nb_scalar, sp.data(), nb_vordiv, vor.data(), + div.data(), rgp1.data() ) ); + + int pos = ( ivar_out * nb_vordiv + jfld ); + + double rms_gen1 = + compute_rms( g1.size(), rgp1.data() + pos * g1.size(), rgp1_analytic.data() ); + + //Log::info() << "rgp1:"; + //for ( int j = 0; j < g1.size(); j++ ) { + // Log::info() << rgp1[pos * g1.size() + j] << " "; + //}; + //Log::info() << std::endl << "analytic1:"; + //for ( int j = 0; j < g1.size(); j++ ) { + // Log::info() << rgp1_analytic[j] << " "; + //}; + //Log::info() << std::endl; + + ATLAS_TRACE_SCOPE( "translocal2" ) + EXPECT_NO_THROW( transLocal2.invtrans( nb_scalar, sp.data(), nb_vordiv, vor.data(), + div.data(), rgp2.data() ) ); + + double rms_gen2 = + compute_rms( g2.size(), rgp2.data() + pos * g2.size(), rgp2_analytic.data() ); + + //Log::info() << std::endl << "rgp2:"; + //for ( int j = 0; j < g2.size(); j++ ) { + // Log::info() << rgp2[pos * g2.size() + j] << " "; + //}; + //Log::info() << std::endl << "analytic2:"; + //for ( int j = 0; j < g2.size(); j++ ) { + // Log::info() << rgp2_analytic[j] << " "; + //}; + //Log::info() << std::endl; + rav1 += rms_gen1; + rav2 += rms_gen2; + if ( !( rms_gen1 < tolerance ) || !( rms_gen2 < tolerance ) ) { + Log::info() + << "Case " << icase << " ivar_in=" << ivar_in << " ivar_out=" << ivar_out + << " m=" << m << " n=" << n << " imag=" << imag << " k=" << k << std::endl; + ATLAS_DEBUG_VAR( rms_gen1 ); + ATLAS_DEBUG_VAR( rms_gen2 ); + ATLAS_DEBUG_VAR( tolerance ); + } + EXPECT( rms_gen1 < tolerance ); + EXPECT( rms_gen2 < tolerance ); + icase++; + auto end = std::chrono::system_clock::now(); // + std::chrono::duration elapsed_seconds = end - start; + std::time_t end_time = std::chrono::system_clock::to_time_t( end ); + std::string time_str = std::ctime( &end_time ); + //Log::info() << "case " << icase << ", elapsed time: " << elapsed_seconds.count() + // << "s. Now: " << time_str.substr( 0, time_str.length() - 1 ) << std::endl; + } + k++; + } + } + } + } + } + } + Log::info() << "Vordiv+scalar comparison with trans: all " << icase << " cases successfully passed!" << std::endl; + rav1 /= icase; + Log::info() << "average RMS error of transLocal1: " << rav1 << std::endl; + rav2 /= icase; + Log::info() << "average RMS error of transLocal2: " << rav2 << std::endl; +} +#endif +//----------------------------------------------------------------------------- +#if 1 CASE( "test_trans_unstructured" ) { Log::info() << "test_trans_unstructured" << std::endl; // test transgeneral by comparing with analytic solution on an unstructured grid From a4e8ed7280c663963c4dafd993b7a18366156fd0 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Wed, 29 Aug 2018 13:09:49 +0100 Subject: [PATCH 19/21] Too strict tolerance for test_trans_pole --- src/tests/trans/test_transgeneral.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests/trans/test_transgeneral.cc b/src/tests/trans/test_transgeneral.cc index 5af4fafcd..a12685406 100644 --- a/src/tests/trans/test_transgeneral.cc +++ b/src/tests/trans/test_transgeneral.cc @@ -875,7 +875,7 @@ CASE( "test_trans_pole" ) { nb_fld = nb_scalar; } else { - tolerance = 2.e-6; + tolerance = 2.e-5; nb_fld = nb_vordiv; } for ( int jfld = 0; jfld < nb_fld; jfld++ ) { // multiple fields From c77fbee5a59b0aeacf37c2166869652f6b679c1b Mon Sep 17 00:00:00 2001 From: Tiago Quintino Date: Fri, 24 Aug 2018 11:33:47 +0100 Subject: [PATCH 20/21] version bump --- VERSION.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION.cmake b/VERSION.cmake index 3465fe8e8..dc6a640b7 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.15.1" ) +set ( ${PROJECT_NAME}_VERSION_STR "0.15.2" ) From 13f013ca5bd2e020d6b09b569dddea61adb1ed61 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 31 Aug 2018 15:51:23 +0100 Subject: [PATCH 21/21] Update changelog --- CHANGELOG.md | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 353ade829..b3ee8a03e 100644 --- a/CHANGELOG.md +++ b/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