diff --git a/src/apps/atlas-meshgen.cc b/src/apps/atlas-meshgen.cc index 1608f077b..ccc349f92 100644 --- a/src/apps/atlas-meshgen.cc +++ b/src/apps/atlas-meshgen.cc @@ -125,6 +125,7 @@ class Meshgen2Gmsh : public AtlasTool { std::string key; long halo; bool edges; + bool cells; bool brick; bool stats; bool info; @@ -165,6 +166,7 @@ Meshgen2Gmsh::Meshgen2Gmsh(int argc, char** argv): AtlasTool(argc, argv) { add_option(new Separator("Advanced")); add_option(new SimpleOption("halo", "Halo size")); add_option(new SimpleOption("edges", "Build edge datastructure")); + add_option(new SimpleOption("cells", "Build cells datastructure")); add_option(new SimpleOption("brick", "Build brick dual mesh")); add_option(new SimpleOption("stats", "Write statistics file")); add_option(new SimpleOption("info", "Write Info")); @@ -203,6 +205,8 @@ int Meshgen2Gmsh::execute(const Args& args) { edges = false; args.get("edges", edges); + cells = false; + args.get("cells", cells); stats = false; args.get("stats", stats); info = false; @@ -280,6 +284,7 @@ int Meshgen2Gmsh::execute(const Args& args) { throw; } + if (grid.projection().units() == "degrees") { functionspace::NodeColumns nodes_fs(mesh, option::halo(halo)); } @@ -299,6 +304,10 @@ int Meshgen2Gmsh::execute(const Args& args) { } } + if (cells) { + functionspace::CellColumns cells_fs(mesh, option::halo(halo)); + } + if (stats) { build_statistics(mesh); } diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 3001357f7..9cf5e9ea8 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -728,6 +728,8 @@ util/detail/BlackMagic.h util/detail/Cache.h util/detail/Debug.h util/detail/KDTree.h +util/function/VortexRollup.h +util/function/VortexRollup.cc ) list( APPEND atlas_internals_srcs diff --git a/src/atlas/functionspace/CellColumns.cc b/src/atlas/functionspace/CellColumns.cc index 6d39e0483..760bc04b5 100644 --- a/src/atlas/functionspace/CellColumns.cc +++ b/src/atlas/functionspace/CellColumns.cc @@ -272,10 +272,10 @@ CellColumns::CellColumns(const Mesh& mesh, const eckit::Configuration& config): mesh::actions::build_cells_parallel_fields(mesh_); mesh::actions::build_periodic_boundaries(mesh_); - if (halo_.size() > 0) { - mesh::actions::build_halo(mesh_, halo_.size()); - nb_cells_ = get_nb_cells_from_metadata(); - } + + mesh::actions::build_halo(mesh_, halo_.size()); + nb_cells_ = get_nb_cells_from_metadata(); + if (!nb_cells_) { nb_cells_ = mesh.cells().size(); } @@ -575,6 +575,10 @@ Field CellColumns::global_index() const { return mesh_.cells().global_index(); } +Field CellColumns::ghost() const { + return mesh_.cells().field("ghost"); +} + //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ diff --git a/src/atlas/functionspace/CellColumns.h b/src/atlas/functionspace/CellColumns.h index 3b904fdd1..1f84c6ea5 100644 --- a/src/atlas/functionspace/CellColumns.h +++ b/src/atlas/functionspace/CellColumns.h @@ -93,6 +93,8 @@ class CellColumns : public functionspace::FunctionSpaceImpl { Field global_index() const override; + Field ghost() const override; + private: // methods idx_t config_size(const eckit::Configuration& config) const; array::DataType config_datatype(const eckit::Configuration&) const; diff --git a/src/atlas/functionspace/CubedSphereColumns.cc b/src/atlas/functionspace/CubedSphereColumns.cc index 649d28481..c2a880b91 100644 --- a/src/atlas/functionspace/CubedSphereColumns.cc +++ b/src/atlas/functionspace/CubedSphereColumns.cc @@ -65,7 +65,7 @@ class CubedSphereStructureCache : public util::Cachesize()); util::ObjectHandle value = Base::get_or_create(key(mesh_impl), creator); return value; } @@ -78,8 +78,8 @@ class CubedSphereStructureCache : public util::Cache(mesh), getGhost(mesh)); + static value_type* create(const Mesh& mesh, idx_t size) { + value_type* value = new value_type(getTij(mesh), getGhost(mesh), size); return value; } }; @@ -121,13 +121,8 @@ idx_t CubedSphereColumns::invalid_index() const { } template -idx_t CubedSphereColumns::nb_elems() const { - return cubedSphereColumnsHandle_.get()->nb_elems(); -} - -template -idx_t CubedSphereColumns::nb_owned_elems() const { - return cubedSphereColumnsHandle_.get()->nb_owned_elems(); +idx_t CubedSphereColumns::sizeOwned() const { + return cubedSphereColumnsHandle_.get()->sizeOwned(); } template @@ -160,11 +155,6 @@ Field CubedSphereColumns::tij() const { return cubedSphereColumnsHandle_.get()->tij(); } -template -Field CubedSphereColumns::ghost() const { - return cubedSphereColumnsHandle_.get()->ghost(); -} - // Explicit instantiation of template classes. template class CubedSphereColumns; template class CubedSphereColumns; diff --git a/src/atlas/functionspace/CubedSphereColumns.h b/src/atlas/functionspace/CubedSphereColumns.h index eb1058d92..916a6c88d 100644 --- a/src/atlas/functionspace/CubedSphereColumns.h +++ b/src/atlas/functionspace/CubedSphereColumns.h @@ -38,11 +38,8 @@ class CubedSphereColumns : public BaseFunctionSpace { /// Invalid index. idx_t invalid_index() const; - /// Get total number of elements. - idx_t nb_elems() const; - /// Get number of owned elements. - idx_t nb_owned_elems() const; + idx_t sizeOwned() const; /// i lower bound for tile t (including halo) idx_t i_begin(idx_t t) const; @@ -60,16 +57,12 @@ class CubedSphereColumns : public BaseFunctionSpace { /// Return tij field. Field tij() const; - /// Return ghost field. - Field ghost() const; - private: class For { public: For(const CubedSphereColumns& functionSpace, const util::Config& config = util::NoConfig()): functionSpace_{functionSpace}, - indexMax_{config.getBool("include_halo", false) ? functionSpace.nb_elems() - : functionSpace.nb_owned_elems()}, + indexMax_{config.getBool("include_halo", false) ? functionSpace.size() : functionSpace.sizeOwned()}, levels_{config.getInt("levels", functionSpace_.levels())}, tijView_(array::make_view(functionSpace_.tij())) {} @@ -78,7 +71,7 @@ class CubedSphereColumns : public BaseFunctionSpace { using EnableFunctor = typename std::enable_if>::value>::type*; - // Functor: void f( index, t, i, j, k) + // Functor: void f(index, t, i, j, k) template = nullptr> void operator()(const Functor& f) const { using namespace meshgenerator::detail::cubedsphere; @@ -94,7 +87,7 @@ class CubedSphereColumns : public BaseFunctionSpace { } } - // Functor: void f( index, t, i, j) + // Functor: void f(index, t, i, j) template = nullptr> void operator()(const Functor& f) const { using namespace meshgenerator::detail::cubedsphere; @@ -108,7 +101,7 @@ class CubedSphereColumns : public BaseFunctionSpace { } } - // Functor: void f( index, k) + // Functor: void f(index, k) template = nullptr> void operator()(const Functor& f) const { using namespace meshgenerator::detail::cubedsphere; @@ -121,7 +114,7 @@ class CubedSphereColumns : public BaseFunctionSpace { } } - // Functor: void f( index ) + // Functor: void f(index ) template = nullptr> void operator()(const Functor& f) const { using namespace meshgenerator::detail::cubedsphere; diff --git a/src/atlas/functionspace/detail/CubedSphereStructure.cc b/src/atlas/functionspace/detail/CubedSphereStructure.cc index 0f1898c94..612825851 100644 --- a/src/atlas/functionspace/detail/CubedSphereStructure.cc +++ b/src/atlas/functionspace/detail/CubedSphereStructure.cc @@ -32,7 +32,8 @@ CubedSphereStructure::BoundingBox::BoundingBox() { jEnd = std::numeric_limits::min(); } -CubedSphereStructure::CubedSphereStructure(const Field& tij, const Field& ghost): tij_(tij), ghost_(ghost) { +CubedSphereStructure::CubedSphereStructure(const Field& tij, const Field& ghost, idx_t size): + tij_(tij), ghost_(ghost), nElems_(size) { ATLAS_TRACE(); Log::debug() << "CubedSphereStructure bounds checking is set to " + std::to_string(checkBounds) << std::endl; @@ -40,8 +41,6 @@ CubedSphereStructure::CubedSphereStructure(const Field& tij, const Field& ghost) const auto tijView_ = array::make_view(tij_); const auto ghostView_ = array::make_view(ghost_); - nElems_ = tijView_.shape(0); - // loop over tij and find min and max ij bounds. for (idx_t index = 0; index < nElems_; ++index) { const size_t t = static_cast(tijView_(index, Coordinates::T)); @@ -76,11 +75,11 @@ CubedSphereStructure::CubedSphereStructure(const Field& tij, const Field& ghost) } } -idx_t CubedSphereStructure::nb_elems() const { +idx_t CubedSphereStructure::size() const { return nElems_; } -idx_t CubedSphereStructure::nb_owned_elems() const { +idx_t CubedSphereStructure::sizeOwned() const { return nOwnedElems_; } diff --git a/src/atlas/functionspace/detail/CubedSphereStructure.h b/src/atlas/functionspace/detail/CubedSphereStructure.h index 81803b466..aaa8354fc 100644 --- a/src/atlas/functionspace/detail/CubedSphereStructure.h +++ b/src/atlas/functionspace/detail/CubedSphereStructure.h @@ -27,16 +27,16 @@ namespace detail { class CubedSphereStructure : public util::Object { public: CubedSphereStructure() = default; - CubedSphereStructure(const Field& tij, const Field& ghost); + CubedSphereStructure(const Field& tij, const Field& ghost, idx_t size); /// Invalid index. static constexpr idx_t invalid_index() { return -1; } /// Number of elements. - idx_t nb_elems() const; + idx_t size() const; /// Number of owned elements. - idx_t nb_owned_elems() const; + idx_t sizeOwned() const; /// i lower bound for tile t (including halo) idx_t i_begin(idx_t) const; diff --git a/src/atlas/grid/detail/partitioner/EqualRegionsPartitioner.cc b/src/atlas/grid/detail/partitioner/EqualRegionsPartitioner.cc index 66feea6b0..92b551577 100644 --- a/src/atlas/grid/detail/partitioner/EqualRegionsPartitioner.cc +++ b/src/atlas/grid/detail/partitioner/EqualRegionsPartitioner.cc @@ -269,6 +269,10 @@ void cap_colats(int N, int n_collars, const double& c_polar, int n_regions[], do c_caps[n_collars + 1] = M_PI; } +// Disable optimisation for Cray (See ATLAS-327) +#ifdef _CRAYC +#pragma _CRI noopt +#endif void eq_caps(int N, std::vector& n_regions, std::vector& s_cap) { // // eq_regions uses the zonal equal area sphere partitioning algorithm to @@ -291,14 +295,14 @@ void eq_caps(int N, std::vector& n_regions, std::vector& s_cap) { // Given N, determine c_polar // the colatitude of the North polar spherical cap. // - double c_polar = polar_colat(N); + const double c_polar = polar_colat(N); // // Given N, determine the ideal angle for spherical collars. // Based on N, this ideal angle, and c_polar, // determine n_collars, the number of collars between the polar caps. // - int n_collars = num_collars(N, c_polar, ideal_collar_angle(N)); + const int n_collars = num_collars(N, c_polar, ideal_collar_angle(N)); // int n_regions_ns=n_collars+2; // @@ -337,6 +341,10 @@ void eq_caps(int N, std::vector& n_regions, std::vector& s_cap) { } // int n_regions_ew=maxval(n_regions(:)); } +// Reenable optimisation for Cray (See ATLAS-327) +#ifdef _CRAYC +#pragma _CRI opt +#endif void eq_regions(int N, double xmin[], double xmax[], double ymin[], double ymax[]) { // EQ_REGIONS Recursive zonal equal area (EQ) partition of sphere diff --git a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc b/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc index 3a52809e5..dee4a9ac4 100644 --- a/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc +++ b/src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc @@ -10,6 +10,8 @@ #include "atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.h" +#include +#include #include #include "eckit/config/Resource.h" @@ -19,6 +21,7 @@ #include "atlas/grid/Iterator.h" #include "atlas/mesh/Nodes.h" #include "atlas/parallel/mpi/mpi.h" +#include "atlas/parallel/omp/fill.h" #include "atlas/runtime/Exception.h" #include "atlas/runtime/Log.h" #include "atlas/util/CoordinateEnums.h" @@ -44,36 +47,63 @@ void MatchingMeshPartitionerLonLatPolygon::partition(const Grid& grid, int parti Log::debug() << "MatchingMeshPartitionerLonLatPolygon::partition" << std::endl; - // FIXME: THIS IS A HACK! the coordinates include North/South Pole (first/last - // partitions only) - bool includesNorthPole = (mpi_rank == 0); - bool includesSouthPole = (mpi_rank == mpi_size - 1); - const util::PolygonXY poly{prePartitionedMesh_.polygon(0)}; + + double west = poly.coordinatesMin().x(); + double east = poly.coordinatesMax().x(); + comm.allReduceInPlace(west, eckit::mpi::Operation::MIN); + comm.allReduceInPlace(east, eckit::mpi::Operation::MAX); + Projection projection = prePartitionedMesh_.projection(); + omp::fill(partitioning, partitioning + grid.size(), -1); - { - eckit::ProgressTimer timer("Partitioning", grid.size(), "point", double(10), atlas::Log::trace()); + auto compute = [&](double west) { size_t i = 0; for (PointLonLat P : grid.lonlat()) { - ++timer; - projection.lonlat2xy(P); - const bool atThePole = (includesNorthPole && P[LAT] >= poly.coordinatesMax()[LAT]) || - (includesSouthPole && P[LAT] < poly.coordinatesMin()[LAT]); - - partitioning[i++] = atThePole || poly.contains(P) ? mpi_rank : -1; + if (partitioning[i] < 0) { + projection.lonlat2xy(P); + P.normalise(west); + partitioning[i] = poly.contains(P) ? mpi_rank : -1; + } + ++i; } - } + // Synchronize partitioning + comm.allReduceInPlace(partitioning, grid.size(), eckit::mpi::Operation::MAX); - // Synchronize partitioning, do a sanity check - comm.allReduceInPlace(partitioning, grid.size(), eckit::mpi::Operation::MAX); - const int min = *std::min_element(partitioning, partitioning + grid.size()); + return *std::min_element(partitioning, partitioning + grid.size()); + }; + + int min = compute(east - 360.); + constexpr double eps = 1.e-10; + if (min < 0 && east - west > 360. + eps) { + min = compute(west - eps); + } if (min < 0) { - throw_Exception( - "Could not find partition for target node (source " - "mesh does not contain all target grid points)", - Here()); + size_t i = 0; + size_t max_failures = grid.size(); + std::vector failed_index; + std::vector failed_lonlat; + failed_index.reserve(max_failures); + failed_lonlat.reserve(max_failures); + for (PointLonLat P : grid.lonlat()) { + if (partitioning[i] < 0) { + failed_index.emplace_back(i); + failed_lonlat.emplace_back(P); + } + ++i; + } + size_t nb_failures = failed_index.size(); + std::stringstream err; + err << "Could not find partition of " << nb_failures + << " target grid points (source " + "mesh does not contain all target grid points)\n" + "Failed target grid points with global index:\n"; + for (size_t n = 0; n < nb_failures; ++n) { + err << " - " << std::setw(10) << std::left << failed_index[n] + 1 << " {lon,lat} : " << failed_lonlat[n] + << "\n"; + } + throw_Exception(err.str(), Here()); } } diff --git a/src/atlas/library/config.h b/src/atlas/library/config.h index 34e9c9949..08fa218b8 100644 --- a/src/atlas/library/config.h +++ b/src/atlas/library/config.h @@ -45,4 +45,12 @@ using idx_t = long; /// @typedef uidx_t /// Integer type for unique indices typedef gidx_t uidx_t; + +#if ATLAS_ECKIT_VERSION_INT >= 11900 // eckit >= 1.19 +#define ATLAS_ECKIT_HAVE_ECKIT_585 1 +#else +#define ATLAS_ECKIT_HAVE_ECKIT_585 0 +#endif + + } // namespace atlas diff --git a/src/atlas/linalg/Introspection.h b/src/atlas/linalg/Introspection.h index 355c2f8d1..a0e37afbb 100644 --- a/src/atlas/linalg/Introspection.h +++ b/src/atlas/linalg/Introspection.h @@ -12,14 +12,11 @@ #include -#include "eckit/linalg/LinearAlgebra.h" #include "eckit/linalg/Matrix.h" -#include "eckit/linalg/SparseMatrix.h" #include "eckit/linalg/Vector.h" +#include "eckit/linalg/types.h" -#include "atlas/array.h" -#include "atlas/parallel/omp/omp.h" -#include "atlas/runtime/Exception.h" +#include "atlas/library/config.h" namespace atlas { namespace linalg { diff --git a/src/atlas/linalg/dense/Backend.cc b/src/atlas/linalg/dense/Backend.cc index 6dda8f0c0..f72a43611 100644 --- a/src/atlas/linalg/dense/Backend.cc +++ b/src/atlas/linalg/dense/Backend.cc @@ -12,7 +12,13 @@ #include +#include "atlas/library/config.h" +#if ATLAS_ECKIT_HAVE_ECKIT_585 +#include "eckit/linalg/LinearAlgebraDense.h" +#else #include "eckit/linalg/LinearAlgebra.h" +#endif + #include "eckit/utils/Tokenizer.h" #include "atlas/library.h" @@ -66,7 +72,11 @@ struct backends { // // However to be identical in behaviour for TransLocal with atlas 0.26.0 and earlier before // any other codes can be adapted in short notice: +#if ATLAS_ECKIT_HAVE_ECKIT_585 + if (eckit::linalg::LinearAlgebraDense::hasBackend("mkl")) { +#else if (eckit::linalg::LinearAlgebra::hasBackend("mkl")) { +#endif current_backend_ = "mkl"; } else { @@ -110,11 +120,19 @@ bool Backend::available() const { std::string t = type(); if (t == backend::eckit_linalg::type()) { if (has("backend")) { +#if ATLAS_ECKIT_HAVE_ECKIT_585 + return eckit::linalg::LinearAlgebraDense::hasBackend(getString("backend")); +#else return eckit::linalg::LinearAlgebra::hasBackend(getString("backend")); +#endif } return true; } +#if ATLAS_ECKIT_HAVE_ECKIT_585 + return eckit::linalg::LinearAlgebraDense::hasBackend(t); +#else return eckit::linalg::LinearAlgebra::hasBackend(t); +#endif } } // namespace dense diff --git a/src/atlas/linalg/dense/MatrixMultiply.h b/src/atlas/linalg/dense/MatrixMultiply.h index 8194cd585..77f4167f0 100644 --- a/src/atlas/linalg/dense/MatrixMultiply.h +++ b/src/atlas/linalg/dense/MatrixMultiply.h @@ -16,6 +16,7 @@ #include "atlas/linalg/Indexing.h" #include "atlas/linalg/View.h" #include "atlas/linalg/dense/Backend.h" +#include "atlas/runtime/Exception.h" #include "atlas/util/Config.h" namespace atlas { diff --git a/src/atlas/linalg/dense/MatrixMultiply.tcc b/src/atlas/linalg/dense/MatrixMultiply.tcc index 23dffe372..a24096a56 100644 --- a/src/atlas/linalg/dense/MatrixMultiply.tcc +++ b/src/atlas/linalg/dense/MatrixMultiply.tcc @@ -12,14 +12,18 @@ #include "MatrixMultiply.h" -#include - #include "atlas/linalg/Indexing.h" #include "atlas/linalg/Introspection.h" #include "atlas/linalg/View.h" #include "atlas/linalg/dense/Backend.h" #include "atlas/runtime/Exception.h" +#if ATLAS_ECKIT_HAVE_ECKIT_585 +#include "eckit/linalg/LinearAlgebraDense.h" +#else +#include "eckit/linalg/LinearAlgebra.h" +#endif + namespace atlas { namespace linalg { @@ -46,16 +50,14 @@ void matrix_multiply( const Mat& A, const Mat& B, Mat& C, const eckit::Configura } else { if( type == "openmp" ) { - // with ECKIT-578 it will be guaranteed that openmp backend is always available. - if( eckit::linalg::LinearAlgebra::hasBackend("openmp") ) { - type = "openmp"; - } - else { - type = "generic"; - } + type = "generic"; dense::MatrixMultiplyHelper::apply( A, B, C, util::Config("backend",type) ); } +#if ATLAS_ECKIT_HAVE_ECKIT_585 + else if( eckit::linalg::LinearAlgebraDense::hasBackend(type) ) { +#else else if( eckit::linalg::LinearAlgebra::hasBackend(type) ) { +#endif dense::MatrixMultiplyHelper::apply( A, B, C, util::Config("backend",type) ); } else { diff --git a/src/atlas/linalg/dense/MatrixMultiply_EckitLinalg.cc b/src/atlas/linalg/dense/MatrixMultiply_EckitLinalg.cc index b972dbaeb..86885e5b5 100644 --- a/src/atlas/linalg/dense/MatrixMultiply_EckitLinalg.cc +++ b/src/atlas/linalg/dense/MatrixMultiply_EckitLinalg.cc @@ -20,7 +20,13 @@ #include "MatrixMultiply_EckitLinalg.h" +#include "atlas/library/config.h" +#if ATLAS_ECKIT_HAVE_ECKIT_585 +#include "eckit/linalg/LinearAlgebraDense.h" +#else #include "eckit/linalg/LinearAlgebra.h" +#endif + #include "eckit/linalg/Matrix.h" #include "eckit/linalg/Vector.h" @@ -31,6 +37,17 @@ namespace linalg { namespace dense { namespace { +#if ATLAS_ECKIT_HAVE_ECKIT_585 +const eckit::linalg::LinearAlgebraDense& eckit_linalg_backend(const Configuration& config) { + std::string backend = "default"; + config.get("backend", backend); + if (backend == "default") { + return eckit::linalg::LinearAlgebraDense::backend(); + } + ATLAS_ASSERT(eckit::linalg::LinearAlgebraDense::hasBackend(backend)); + return eckit::linalg::LinearAlgebraDense::getBackend(backend); +} +#else const eckit::linalg::LinearAlgebra& eckit_linalg_backend(const Configuration& config) { std::string backend = "default"; config.get("backend", backend); @@ -40,6 +57,7 @@ const eckit::linalg::LinearAlgebra& eckit_linalg_backend(const Configuration& co ATLAS_ASSERT(eckit::linalg::LinearAlgebra::hasBackend(backend)); return eckit::linalg::LinearAlgebra::getBackend(backend); } +#endif } // namespace diff --git a/src/atlas/linalg/sparse/Backend.cc b/src/atlas/linalg/sparse/Backend.cc index 1ce1e55ae..b9084884a 100644 --- a/src/atlas/linalg/sparse/Backend.cc +++ b/src/atlas/linalg/sparse/Backend.cc @@ -12,7 +12,12 @@ #include +#include "atlas/library/config.h" +#if ATLAS_ECKIT_HAVE_ECKIT_585 +#include "eckit/linalg/LinearAlgebraSparse.h" +#else #include "eckit/linalg/LinearAlgebra.h" +#endif #include "eckit/utils/Tokenizer.h" #include "atlas/library.h" @@ -96,11 +101,19 @@ bool Backend::available() const { } if (t == backend::eckit_linalg::type()) { if (has("backend")) { +#if ATLAS_ECKIT_HAVE_ECKIT_585 + return eckit::linalg::LinearAlgebraSparse::hasBackend(getString("backend")); +#else return eckit::linalg::LinearAlgebra::hasBackend(getString("backend")); +#endif } return true; } +#if ATLAS_ECKIT_HAVE_ECKIT_585 + return eckit::linalg::LinearAlgebraSparse::hasBackend(t); +#else return eckit::linalg::LinearAlgebra::hasBackend(t); +#endif } diff --git a/src/atlas/linalg/sparse/SparseMatrixMultiply.h b/src/atlas/linalg/sparse/SparseMatrixMultiply.h index 899b200d6..d1776975b 100644 --- a/src/atlas/linalg/sparse/SparseMatrixMultiply.h +++ b/src/atlas/linalg/sparse/SparseMatrixMultiply.h @@ -16,6 +16,7 @@ #include "atlas/linalg/Indexing.h" #include "atlas/linalg/View.h" #include "atlas/linalg/sparse/Backend.h" +#include "atlas/runtime/Exception.h" #include "atlas/util/Config.h" namespace atlas { diff --git a/src/atlas/linalg/sparse/SparseMatrixMultiply.tcc b/src/atlas/linalg/sparse/SparseMatrixMultiply.tcc index 0e1ca4506..b4e1ed9cd 100644 --- a/src/atlas/linalg/sparse/SparseMatrixMultiply.tcc +++ b/src/atlas/linalg/sparse/SparseMatrixMultiply.tcc @@ -12,14 +12,19 @@ #include "SparseMatrixMultiply.h" -#include - #include "atlas/linalg/Indexing.h" #include "atlas/linalg/Introspection.h" #include "atlas/linalg/View.h" #include "atlas/linalg/sparse/Backend.h" #include "atlas/runtime/Exception.h" +#if ATLAS_ECKIT_HAVE_ECKIT_585 +#include "eckit/linalg/LinearAlgebraSparse.h" +#else +#include "eckit/linalg/LinearAlgebra.h" +#endif + + namespace atlas { namespace linalg { @@ -75,7 +80,11 @@ void sparse_matrix_multiply( const Matrix& matrix, const SourceView& src, Target else if ( type == sparse::backend::eckit_linalg::type() ) { sparse::dispatch_sparse_matrix_multiply( matrix, src, tgt, indexing, config ); } +#if ATLAS_ECKIT_HAVE_ECKIT_585 + else if( eckit::linalg::LinearAlgebraSparse::hasBackend(type) ) { +#else else if( eckit::linalg::LinearAlgebra::hasBackend(type) ) { +#endif sparse::dispatch_sparse_matrix_multiply( matrix, src, tgt, indexing, util::Config("backend",type) ); } else { diff --git a/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.cc b/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.cc index 72c0b4ecb..72d55c1a9 100644 --- a/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.cc +++ b/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.cc @@ -20,7 +20,14 @@ #include "SparseMatrixMultiply_EckitLinalg.h" +#include "atlas/library/config.h" + +#if ATLAS_ECKIT_HAVE_ECKIT_585 +#include "eckit/linalg/LinearAlgebraSparse.h" +#else #include "eckit/linalg/LinearAlgebra.h" +#endif + #include "eckit/linalg/Matrix.h" #include "eckit/linalg/Vector.h" @@ -31,6 +38,19 @@ namespace linalg { namespace sparse { namespace { + +#if ATLAS_ECKIT_HAVE_ECKIT_585 +const eckit::linalg::LinearAlgebraSparse& eckit_linalg_backend(const Configuration& config) { + std::string backend = "default"; + config.get("backend", backend); + + if (backend == "default") { + return eckit::linalg::LinearAlgebraSparse::backend(); + } + ATLAS_ASSERT(eckit::linalg::LinearAlgebraSparse::hasBackend(backend)); + return eckit::linalg::LinearAlgebraSparse::getBackend(backend); +} +#else const eckit::linalg::LinearAlgebra& eckit_linalg_backend(const Configuration& config) { std::string backend = "default"; config.get("backend", backend); @@ -40,6 +60,7 @@ const eckit::linalg::LinearAlgebra& eckit_linalg_backend(const Configuration& co ATLAS_ASSERT(eckit::linalg::LinearAlgebra::hasBackend(backend)); return eckit::linalg::LinearAlgebra::getBackend(backend); } +#endif } // namespace diff --git a/src/atlas/mesh/actions/BuildEdges.cc b/src/atlas/mesh/actions/BuildEdges.cc index 198efe062..ac773bb01 100644 --- a/src/atlas/mesh/actions/BuildEdges.cc +++ b/src/atlas/mesh/actions/BuildEdges.cc @@ -289,89 +289,6 @@ class AccumulatePoleEdges { } }; -void accumulate_pole_edges(mesh::Nodes& nodes, std::vector& pole_edge_nodes, idx_t& nb_pole_edges) { - enum - { - NORTH = 0, - SOUTH = 1 - }; - - const auto xy = array::make_view(nodes.xy()); - const auto flags = array::make_view(nodes.flags()); - const auto part = array::make_view(nodes.partition()); - const idx_t nb_nodes = nodes.size(); - - double min[2], max[2]; - min[XX] = std::numeric_limits::max(); - min[YY] = std::numeric_limits::max(); - max[XX] = -std::numeric_limits::max(); - max[YY] = -std::numeric_limits::max(); - for (idx_t node = 0; node < nb_nodes; ++node) { - min[XX] = std::min(min[XX], xy(node, XX)); - min[YY] = std::min(min[YY], xy(node, YY)); - max[XX] = std::max(max[XX], xy(node, XX)); - max[YY] = std::max(max[YY], xy(node, YY)); - } - - ATLAS_TRACE_MPI(ALLREDUCE) { - mpi::comm().allReduceInPlace(min, 2, eckit::mpi::min()); - mpi::comm().allReduceInPlace(max, 2, eckit::mpi::max()); - } - - double tol = 1e-6; - - // Collect all nodes closest to poles - std::vector> pole_nodes(2); - for (idx_t node = 0; node < nb_nodes; ++node) { - if (std::abs(xy(node, YY) - max[YY]) < tol) { - pole_nodes[NORTH].insert(node); - } - else if (std::abs(xy(node, YY) - min[YY]) < tol) { - pole_nodes[SOUTH].insert(node); - } - } - - // Sanity check - { - for (idx_t NS = 0; NS < 2; ++NS) { - int npart = -1; - for (std::set::iterator it = pole_nodes[NS].begin(); it != pole_nodes[NS].end(); ++it) { - int node = *it; - if (npart == -1) { - npart = part(node); - } - else if (part(node) != npart) { - // Not implemented yet, when pole-lattitude is split. - std::stringstream msg; - msg << "Split pole-latitude is not supported yet... node " << node << "[p" << part(node) - << "] should belong to part " << npart; - throw_NotImplemented(msg.str(), Here()); - } - } - } - } - - // Create connections over the poles and store in pole_edge_nodes - nb_pole_edges = 0; - for (idx_t NS = 0; NS < 2; ++NS) { - for (std::set::iterator it = pole_nodes[NS].begin(); it != pole_nodes[NS].end(); ++it) { - int node = *it; - if (!Topology::check(flags(node), Topology::PERIODIC | Topology::GHOST)) { - int x2 = microdeg(xy(node, XX) + 180.); - for (std::set::iterator itr = pole_nodes[NS].begin(); itr != pole_nodes[NS].end(); ++itr) { - int other_node = *itr; - if (microdeg(xy(other_node, XX)) == x2) { - if (!Topology::check(flags(other_node), Topology::PERIODIC)) { - pole_edge_nodes.push_back(node); - pole_edge_nodes.push_back(other_node); - ++nb_pole_edges; - } - } - } - } - } - } -} struct ComputeUniquePoleEdgeIndex { // Already assumes that the edges cross the pole @@ -456,11 +373,6 @@ void build_edges(Mesh& mesh, const eckit::Configuration& config) { accumulate_facets_ordered_by_halo(mesh.cells(), mesh.nodes(), edge_nodes_data, edge_to_elem_data, nb_edges, nb_inner_edges, missing_value, edge_halo_offsets); - std::shared_ptr pole_edge_accumulator; - if (pole_edges) { - pole_edge_accumulator = std::make_shared(nodes); - } - std::vector sorted_edge_nodes_data; std::vector sorted_edge_to_elem_data; @@ -551,6 +463,8 @@ void build_edges(Mesh& mesh, const eckit::Configuration& config) { edge_to_elem_data.data() + edge_halo_offsets[halo] * 2); if (pole_edges) { + auto pole_edge_accumulator = std::make_shared(nodes); + idx_t nb_pole_edges; std::vector pole_edge_nodes; diff --git a/src/atlas/mesh/actions/BuildParallelFields.cc b/src/atlas/mesh/actions/BuildParallelFields.cc index 10d67a3aa..f096419a9 100644 --- a/src/atlas/mesh/actions/BuildParallelFields.cc +++ b/src/atlas/mesh/actions/BuildParallelFields.cc @@ -10,6 +10,7 @@ #include +#include #include #include @@ -981,8 +982,6 @@ Field& build_cells_remote_idx(mesh::Cells& cells, const mesh::Nodes& nodes) { idx_t mypart = static_cast(mpi::rank()); idx_t nparts = static_cast(mpi::size()); - UniqueLonLat compute_uid(nodes); - // This piece should be somewhere central ... could be NPROMA ? // ----------> std::vector proc(nparts); @@ -992,10 +991,29 @@ Field& build_cells_remote_idx(mesh::Cells& cells, const mesh::Nodes& nodes) { // <--------- auto ridx = array::make_indexview(cells.remote_index()); - auto part = array::make_view(cells.partition()); + const auto part = array::make_view(cells.partition()); + const auto gidx = array::make_view(cells.global_index()); + const auto flags = array::make_view(cells.flags()); const auto& element_nodes = cells.node_connectivity(); idx_t nb_cells = cells.size(); + const PeriodicTransform transform_periodic_east(-360.); + const PeriodicTransform transform_periodic_west(+360.); + const UniqueLonLat compute_uid_lonlat(nodes); + auto compute_uid = [&](idx_t jcell) { + constexpr int PERIODIC = util::Topology::PERIODIC; + constexpr int EAST = util::Topology::EAST; + constexpr int WEST = util::Topology::WEST; + const auto flags_view = util::Bitflags::view(flags(jcell)); + if (flags_view.check(PERIODIC | EAST)) { + return compute_uid_lonlat(element_nodes.row(jcell), transform_periodic_east); + } + if (flags_view.check(PERIODIC | WEST)) { + return compute_uid_lonlat(element_nodes.row(jcell), transform_periodic_west); + } + return compute_uid_lonlat(element_nodes.row(jcell)); + }; + idx_t varsize = 2; std::vector> send_needed(mpi::size()); @@ -1003,7 +1021,7 @@ Field& build_cells_remote_idx(mesh::Cells& cells, const mesh::Nodes& nodes) { int sendcnt = 0; std::map lookup; for (idx_t jcell = 0; jcell < nb_cells; ++jcell) { - uid_t uid = compute_uid(element_nodes.row(jcell)); + uid_t uid = compute_uid(jcell); if (idx_t(part(jcell)) == mypart) { lookup[uid] = jcell; @@ -1034,38 +1052,39 @@ Field& build_cells_remote_idx(mesh::Cells& cells, const mesh::Nodes& nodes) { for (idx_t jpart = 0; jpart < nparts; ++jpart) { const std::vector& recv_cell = recv_needed[proc[jpart]]; const idx_t nb_recv_cells = idx_t(recv_cell.size()) / varsize; - // array::ArrayView recv_node( make_view( Array::wrap(shape, - // recv_needed[ proc[jpart] ].data()) ), - // array::make_shape(recv_needed[ proc[jpart] ].size()/varsize,varsize) - // ); for (idx_t jcell = 0; jcell < nb_recv_cells; ++jcell) { uid_t uid = recv_cell[jcell * varsize + 0]; int icell = recv_cell[jcell * varsize + 1]; - if (lookup.count(uid)) { - send_found[proc[jpart]].push_back(icell); - send_found[proc[jpart]].push_back(lookup[uid]); - } - else { - std::stringstream msg; - msg << "[" << mpi::rank() << "] " - << "Node requested by rank [" << jpart << "] with uid [" << uid - << "] that should be owned is not found"; - throw_Exception(msg.str(), Here()); - } + send_found[proc[jpart]].push_back(icell); + send_found[proc[jpart]].push_back(lookup.count(uid) ? lookup[uid] : -1); } } ATLAS_TRACE_MPI(ALLTOALL) { mpi::comm().allToAll(send_found, recv_found); } + std::stringstream errstream; + size_t failed{0}; for (idx_t jpart = 0; jpart < nparts; ++jpart) { const std::vector& recv_cell = recv_found[proc[jpart]]; const idx_t nb_recv_cells = recv_cell.size() / 2; // array::ArrayView recv_node( recv_found[ proc[jpart] ].data(), // array::make_shape(recv_found[ proc[jpart] ].size()/2,2) ); for (idx_t jcell = 0; jcell < nb_recv_cells; ++jcell) { - ridx(recv_cell[jcell * 2 + 0]) = recv_cell[jcell * 2 + 1]; + idx_t icell = recv_cell[jcell * 2 + 0]; + idx_t ridx_icell = recv_cell[jcell * 2 + 1]; + if (ridx_icell >= 0) { + ridx(icell) = ridx_icell; + } + else { + ++failed; + errstream << "\n[" << mpi::rank() << "] " + << "Cell " << gidx(icell) << " not found on part [" << part(icell) << "]"; + } } } + if (failed) { + throw_AssertionFailed(errstream.str(), Here()); + } return cells.remote_index(); } diff --git a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc index b7dde8f4a..b6d3eea46 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc @@ -189,17 +189,16 @@ void CubedSphereDualMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, copyField(csCells.partition(), nodes.partition()); copyField(csCells.halo(), nodes.halo()); copyField(csCells.flags(), nodes.flags()); + copyField(csCells.field("ghost"), nodes.ghost()); copyField(csCells.field("xy"), nodes.xy()); copyField(csCells.field("lonlat"), nodes.lonlat()); copyField(csCells.field("tij"), nodes.field("tij")); - // Need to set a ghost field and decrement halo by one. - auto nodesGhost = array::make_view(nodes.ghost()); - auto nodesHalo = array::make_view(nodes.halo()); + // Need to decrement halo by one. + auto nodesHalo = array::make_view(nodes.halo()); for (idx_t idx = 0; idx < nodes.size(); ++idx) { - nodesGhost(idx) = nodesHalo(idx) > 0; - nodesHalo(idx) = std::max(0, nodesHalo(idx) - 1); + nodesHalo(idx) = std::max(0, nodesHalo(idx) - 1); } //-------------------------------------------------------------------------- @@ -278,10 +277,13 @@ void CubedSphereDualMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cells.add(Field("lonlat", array::make_datatype(), array::make_shape(cells.size(), 2))); + cells.add(Field("ghost", array::make_datatype(), array::make_shape(cells.size()))); + // Copy dual cells fields from nodes. copyField(csNodes.global_index(), cells.global_index()); copyField(csNodes.remote_index(), cells.remote_index()); copyField(csNodes.partition(), cells.partition()); + copyField(csNodes.ghost(), cells.field("ghost")); copyField(csNodes.halo(), cells.halo()); copyField(csNodes.flags(), cells.flags()); copyField(csNodes.field("tij"), cells.field("tij")); diff --git a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc index 5a080d19d..6dad6bbc7 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc @@ -830,10 +830,13 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons Field lonLatField = cells.add(Field("lonlat", make_datatype(), make_shape(cells.size(), 2))); lonLatField.set_variables(2); + Field ghostField = cells.add(Field("ghost", make_datatype(), make_shape(cells.size()))); + // Set field views. auto cellsGlobalIdx = array::make_view(cells.global_index()); auto cellsRemoteIdx = array::make_indexview(cells.remote_index()); auto cellsPart = array::make_view(cells.partition()); + auto cellsGhost = array::make_view(ghostField); auto cellsHalo = array::make_view(cells.halo()); auto cellsFlags = array::make_view(cells.flags()); auto cellsTij = array::make_view(tijField); @@ -894,6 +897,9 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons // Set halo. cellsHalo(cellLocalIdx) = localCell.halo; + // Set ghost. + cellsGhost(cellLocalIdx) = cellsHalo(cellLocalIdx) > 0; + // Set flags. Topology::reset(cellsFlags(cellLocalIdx)); switch (localCell.type) { @@ -921,16 +927,48 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons // Done. That was rather a lot of bookkeeping! // --------------------------------------------------------------------------- + set_metadata(mesh); + + return; +} + +// ----------------------------------------------------------------------------- + +void CubedSphereMeshGenerator::set_metadata(Mesh& mesh) const { + const auto nHalo = options.get("halo"); + + // Set basic halo metadata. mesh.metadata().set("halo", nHalo); mesh.metadata().set("halo_locked", true); - mesh.nodes().metadata().set("parallel", true); - mesh.nodes().metadata().set("nb_owned", nodeLocalIdxCount[thisPart]); - mesh.cells().metadata().set("parallel", true); - mesh.cells().metadata().set("nb_owned", cellLocalIdxCount[thisPart]); - return; + // Loop over nodes and count number of halo elements. + auto nNodes = std::vector(nHalo + 1, 0); + const auto nodeHalo = array::make_view(mesh.nodes().halo()); + for (idx_t i = 0; i < mesh.nodes().size(); ++i) { + ++nNodes[static_cast(nodeHalo(i))]; + } + std::partial_sum(nNodes.begin(), nNodes.end(), nNodes.begin()); + + // Set node halo metadata. + for (size_t i = 0; i < nNodes.size(); ++i) { + const auto str = "nb_nodes_including_halo[" + std::to_string(i) + "]"; + mesh.metadata().set(str, nNodes[i]); + } + // Loop over cells and count number of halo elements. + auto nCells = std::vector(nHalo + 1, 0); + const auto cellHalo = array::make_view(mesh.cells().halo()); + for (idx_t i = 0; i < mesh.cells().size(); ++i) { + ++nCells[static_cast(cellHalo(i))]; + } + std::partial_sum(nCells.begin(), nCells.end(), nCells.begin()); + + // Set cell halo metadata. + for (size_t i = 0; i < nCells.size(); ++i) { + const auto str = "nb_cells_including_halo[0][" + std::to_string(i) + "]"; + mesh.metadata().set(str, nCells[i]); + } } // ----------------------------------------------------------------------------- diff --git a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.h b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.h index 9bfa0e40d..a350ea847 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.h +++ b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.h @@ -56,6 +56,8 @@ class CubedSphereMeshGenerator : public MeshGenerator::Implementation { void generate_mesh(const CubedSphereGrid&, const grid::Distribution&, Mesh&) const; + void set_metadata(Mesh&) const; + private: util::Metadata options; }; diff --git a/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc b/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc index 63ee66136..3926db641 100644 --- a/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc @@ -369,6 +369,8 @@ void HealpixMeshGenerator::hash(eckit::Hash& h) const { } void HealpixMeshGenerator::generate(const Grid& grid, const grid::Distribution& distribution, Mesh& mesh) const { + ATLAS_TRACE(); + Log::debug() << "HealpixMeshGenerator generating mesh from " << grid.name() << std::endl; ATLAS_ASSERT(HealpixGrid(grid), "Grid could not be cast to a HealpixGrid"); ATLAS_ASSERT(!mesh.generated()); @@ -643,6 +645,7 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: auto remote_idx = array::make_indexview(nodes.remote_index()); auto part = array::make_view(nodes.partition()); auto ghost = array::make_view(nodes.ghost()); + auto halo = array::make_view(nodes.halo()); auto flags = array::make_view(nodes.flags()); // define cells and associated properties @@ -739,6 +742,7 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: part(inode) = parts_SB[iil]; ghost(inode) = is_ghost_SB[iil]; + halo(inode) = 0; if (ghost(inode)) { Topology::set(flags(inode), Topology::GHOST); diff --git a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc index 2d21bfa4e..a32d43a9e 100644 --- a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc @@ -218,6 +218,7 @@ void StructuredMeshGenerator::hash(eckit::Hash& h) const { void StructuredMeshGenerator::generate(const Grid& grid, const grid::Distribution& distribution, Mesh& mesh) const { ATLAS_TRACE(); + Log::debug() << "StructuredMeshGenerator generating mesh from " << grid.name() << std::endl; const StructuredGrid rg = StructuredGrid(grid); if (!rg) { diff --git a/src/atlas/projection/detail/SchmidtProjection.cc b/src/atlas/projection/detail/SchmidtProjection.cc index 593aee191..16a6cb7aa 100644 --- a/src/atlas/projection/detail/SchmidtProjection.cc +++ b/src/atlas/projection/detail/SchmidtProjection.cc @@ -7,7 +7,6 @@ * granted to it by virtue of its status as an intergovernmental organisation * nor does it submit to any jurisdiction. */ - #include #include "eckit/config/Parametrisation.h" @@ -33,6 +32,14 @@ namespace atlas { namespace projection { namespace detail { +PointLonLat rotation_north_pole(const Rotated& rotation) { + return rotation.northPole(); +} + +PointLonLat rotation_north_pole(const NotRotated& rotation) { + return PointLonLat{0., 90.}; +} + // constructor template SchmidtProjectionT::SchmidtProjectionT(const eckit::Parametrisation& params): @@ -40,6 +47,12 @@ SchmidtProjectionT::SchmidtProjectionT(const eckit::Parametrisation& p if (!params.get("stretching_factor", c_)) { throw_Exception("stretching_factor missing in Params", Here()); } + ATLAS_ASSERT(c_ != 0.); + + north0_ = {0.0, 0.0, 1.0}; + + atlas::util::UnitSphere::convertSphericalToCartesian(rotation_north_pole(rotation_), north1_); + north1_ = PointXYZ::normalize(north1_); } // constructor @@ -75,36 +88,34 @@ ProjectionImpl::Jacobian SchmidtProjectionT::jacobian(const PointLonLa lonlat2xy(xy); - PointXYZ xyz, north1, north0(0.0, 0.0, 1.0); - atlas::util::UnitSphere::convertSphericalToCartesian(rotation_.northPole(), north1); + PointXYZ xyz; atlas::util::UnitSphere::convertSphericalToCartesian(lonlat, xyz); - north1 = PointXYZ::normalize(north1); + + double zomc2 = 1.0 - 1.0 / (c_ * c_); + double zopc2 = 1.0 + 1.0 / (c_ * c_); + + double zcosy = std::cos(D2R(xy[1])); + double zsiny = std::sin(D2R(xy[1])); + double zcoslat = std::cos(D2R(lonlat.lat())); + + double zfactor = std::sqrt((zopc2 + zsiny * zomc2) * (zopc2 + zsiny * zomc2) / (zopc2 * zopc2 - zomc2 * zomc2)); // Base vectors in unrotated frame - auto u0 = PointXYZ::normalize(PointXYZ::cross(north0, xyz)); + auto u0 = PointXYZ::normalize(PointXYZ::cross(north0_, xyz)); auto v0 = PointXYZ::normalize(PointXYZ::cross(xyz, u0)); // Base vectors in rotated frame - auto u1 = PointXYZ::normalize(PointXYZ::cross(north1, xyz)); + auto u1 = PointXYZ::normalize(PointXYZ::cross(north1_, xyz)); auto v1 = PointXYZ::normalize(PointXYZ::cross(xyz, u1)); double u0u1 = PointXYZ::dot(u0, u1); double v0u1 = PointXYZ::dot(v0, u1); + double u0v1 = PointXYZ::dot(u0, v1); double v0v1 = PointXYZ::dot(v0, v1); Jacobian jac; - - double zomc2 = 1.0 - 1.0 / (c_ * c_); - double zopc2 = 1.0 + 1.0 / (c_ * c_); - - double zcosy = cos(D2R(xy[1])), zsiny = sin(D2R(xy[1])); - double zcoslat = cos(D2R(lonlat.lat())); - - double zfactor = sqrt((zopc2 + zsiny * zomc2) * (zopc2 + zsiny * zomc2) / (zopc2 * zopc2 - zomc2 * zomc2)); - - jac[0] = {zcoslat * u0u1 * zfactor / zcosy, v0u1 * zfactor / zcosy}; jac[1] = {zcoslat * u0v1 * zfactor, v0v1 * zfactor}; diff --git a/src/atlas/projection/detail/SchmidtProjection.h b/src/atlas/projection/detail/SchmidtProjection.h index 92200e430..24509e795 100644 --- a/src/atlas/projection/detail/SchmidtProjection.h +++ b/src/atlas/projection/detail/SchmidtProjection.h @@ -49,6 +49,8 @@ class SchmidtProjectionT final : public ProjectionImpl { private: double c_; // stretching factor Rotation rotation_; + PointXYZ north0_; + PointXYZ north1_; }; using SchmidtProjection = SchmidtProjectionT; diff --git a/src/atlas/projection/detail/VariableResolutionProjection.cc b/src/atlas/projection/detail/VariableResolutionProjection.cc index 6d15313f9..64317b4a0 100644 --- a/src/atlas/projection/detail/VariableResolutionProjection.cc +++ b/src/atlas/projection/detail/VariableResolutionProjection.cc @@ -45,6 +45,12 @@ * */ +#ifdef __NVCOMPILER +#define PREVENT_OPT volatile +#else +#define PREVENT_OPT +#endif + namespace atlas { namespace projection { namespace detail { @@ -59,10 +65,11 @@ static double new_ratio(int n_stretched, double var_ratio) { constexpr float epstest = std::numeric_limits::epsilon(); ///< number of variable (stretched) grid points in one side - int var_ints = (n_stretched + epstest) / 2.; - double var_ints_f = n_stretched / 2.; - double logr = std::log(var_ratio); - double log_ratio = (var_ints_f - 0.5) * logr; + PREVENT_OPT int var_ints = (n_stretched + epstest) / 2.; + double var_ints_f = n_stretched / 2.; + double logr = std::log(var_ratio); + double log_ratio = (var_ints_f - 0.5) * logr; + return std::exp(log_ratio / var_ints); }; @@ -85,7 +92,6 @@ typename VariableResolutionProjectionT::Spec VariableResolutionProject proj_st.set("endy", endy_); ///< original domain endy proj_st.set("rim_widthx", rim_widthx_); ///< xsize of the rim proj_st.set("rim_widthy", rim_widthy_); ///< ysize of the rim - proj_st.set("north_pole", north_pole_); ///< north_pole rotation_.spec(proj_st); return proj_st; } @@ -95,20 +101,19 @@ typename VariableResolutionProjectionT::Spec VariableResolutionProject template VariableResolutionProjectionT::VariableResolutionProjectionT(const eckit::Parametrisation& proj_st): ProjectionImpl(), rotation_(proj_st) { - proj_st.get("delta_low", delta_outer = 0.); ///< resolution of the external regular grid (rim) - proj_st.get("delta_hi", delta_inner = 0.); ///< resolution of the regional model (regular grid) - proj_st.get("var_ratio", var_ratio_ = 0.); ///< power used for the stretching - proj_st.get("x_reg_start", x_reg_start_ = 0.); ///< xstart of the internal regional grid - proj_st.get("y_reg_start", y_reg_start_ = 0.); ///< ystart of the internal regional grid - proj_st.get("x_reg_end", x_reg_end_ = 0.); ///< xend of the regular part of stretched internal grid - proj_st.get("y_reg_end", y_reg_end_ = 0.); ///< yend of the regular part of stretched internal grid - proj_st.get("startx", startx_ = 0.); ///< original domain startx - proj_st.get("endx", endx_ = 0.); ///< original domain endx - proj_st.get("starty", starty_ = 0.); ///< original domain starty - proj_st.get("endy", endy_ = 0.); ///< original domain endy - proj_st.get("rim_widthx", rim_widthx_); ///< xsize of the rim - proj_st.get("rim_widthy", rim_widthy_); ///< ysize of the rim - proj_st.get("north_pole", north_pole_ = {0., 90.}); ///< north_pole + proj_st.get("delta_low", delta_outer = 0.); ///< resolution of the external regular grid (rim) + proj_st.get("delta_hi", delta_inner = 0.); ///< resolution of the regional model (regular grid) + proj_st.get("var_ratio", var_ratio_ = 0.); ///< power used for the stretching + proj_st.get("x_reg_start", x_reg_start_ = 0.); ///< xstart of the internal regional grid + proj_st.get("y_reg_start", y_reg_start_ = 0.); ///< ystart of the internal regional grid + proj_st.get("x_reg_end", x_reg_end_ = 0.); ///< xend of the regular part of stretched internal grid + proj_st.get("y_reg_end", y_reg_end_ = 0.); ///< yend of the regular part of stretched internal grid + proj_st.get("startx", startx_ = 0.); ///< original domain startx + proj_st.get("endx", endx_ = 0.); ///< original domain endx + proj_st.get("starty", starty_ = 0.); ///< original domain starty + proj_st.get("endy", endy_ = 0.); ///< original domain endy + proj_st.get("rim_widthx", rim_widthx_); ///< xsize of the rim + proj_st.get("rim_widthy", rim_widthy_); ///< ysize of the rim constexpr float epsilon = std::numeric_limits::epsilon(); ///< value used to check if the values are equal constexpr float epstest = diff --git a/src/atlas/projection/detail/VariableResolutionProjection.h b/src/atlas/projection/detail/VariableResolutionProjection.h index fdbbf26e3..f9f3076a7 100644 --- a/src/atlas/projection/detail/VariableResolutionProjection.h +++ b/src/atlas/projection/detail/VariableResolutionProjection.h @@ -76,6 +76,7 @@ class VariableResolutionProjectionT final : public ProjectionImpl { std::vector north_pole_; ///< north_pole lon,lat + //< variables derived from the configuration used for the projection double deltax_all, deltay_all; double add_xf_, add_yf_; diff --git a/src/atlas/util/Point.h b/src/atlas/util/Point.h index 2f8443463..91994a221 100644 --- a/src/atlas/util/Point.h +++ b/src/atlas/util/Point.h @@ -95,6 +95,20 @@ class PointXYZ : public eckit::geometry::Point3 { x_[1] = y; x_[2] = z; } + + PointXYZ& operator*=(double m) { + x_[0] *= m; + x_[1] *= m; + x_[2] *= m; + return *this; + } + + PointXYZ& operator/=(double m) { + x_[0] /= m; + x_[1] /= m; + x_[2] /= m; + return *this; + } }; /// @brief Point in longitude-latitude coordinate system. diff --git a/src/atlas/util/Polygon.cc b/src/atlas/util/Polygon.cc index 2bfa0b84f..40a88cb63 100644 --- a/src/atlas/util/Polygon.cc +++ b/src/atlas/util/Polygon.cc @@ -152,7 +152,7 @@ PolygonCoordinates::PolygonCoordinates(const Polygon& poly, const atlas::Field& if ((coordinates_.size() >= 2) && removeAlignedPoints) { const Point2& B = coordinates_.back(); const Point2& C = coordinates_[coordinates_.size() - 2]; - if (eckit::types::is_approximately_equal(0., cross_product_analog(A, B, C))) { + if (eckit::types::is_approximately_equal(0., cross_product_analog(A, B, C), 1.e-10)) { coordinates_.back() = A; ++nb_removed_points_due_to_alignment; continue; @@ -208,7 +208,7 @@ PolygonCoordinates::PolygonCoordinates(const PointContainer& points, bool remove if ((coordinates_.size() >= 2) && removeAlignedPoints) { const Point2& B = coordinates_.back(); const Point2& C = coordinates_[coordinates_.size() - 2]; - if (eckit::types::is_approximately_equal(0., cross_product_analog(A, B, C))) { + if (eckit::types::is_approximately_equal(0., cross_product_analog(A, B, C), 1.e-10)) { coordinates_.back() = A; ++nb_removed_points_due_to_alignment; continue; @@ -257,6 +257,11 @@ void PolygonCoordinates::print(std::ostream& out) const { out << "]"; } +std::ostream& operator<<(std::ostream& out, const PolygonCoordinates& pc) { + pc.print(out); + return out; +} + Polygon::edge_set_t ExplicitPartitionPolygon::compute_edges(idx_t points_size) { util::Polygon::edge_set_t edges; auto add_edge = [&](idx_t p1, idx_t p2) { diff --git a/src/atlas/util/Polygon.h b/src/atlas/util/Polygon.h index 8ff6f1134..55d7a4043 100644 --- a/src/atlas/util/Polygon.h +++ b/src/atlas/util/Polygon.h @@ -184,10 +184,17 @@ class PolygonCoordinates { const Point2& coordinatesMin() const; const Point2& centroid() const; + template + const Point2& operator[](Index i) const { + return coordinates_[i]; + } + idx_t size() const { return coordinates_.size(); } void print(std::ostream&) const; + friend std::ostream& operator<<(std::ostream& out, const PolygonCoordinates& pc); + protected: // -- Members diff --git a/src/atlas/util/PolygonXY.cc b/src/atlas/util/PolygonXY.cc index d79e4b044..190188904 100644 --- a/src/atlas/util/PolygonXY.cc +++ b/src/atlas/util/PolygonXY.cc @@ -28,7 +28,7 @@ namespace util { namespace { double cross_product_analog(const Point2& A, const Point2& B, const Point2& C) { - return (A[LON] - C[LON]) * (B[LAT] - C[LAT]) - (A[LAT] - C[LAT]) * (B[LON] - C[LON]); + return (A.x() - C.x()) * (B.y() - C.y()) - (A.y() - C.y()) * (B.x() - C.x()); } @@ -90,7 +90,7 @@ PolygonXY::PolygonXY(const PartitionPolygon& partition_polygon): PolygonCoordina centroid_ = compute_centroid(coordinates_); inner_radius_squared_ = compute_inner_radius_squared(coordinates_, centroid_); - ATLAS_ASSERT(contains(centroid_)); + ATLAS_ASSERT(PolygonXY::contains(centroid_)); } } @@ -102,14 +102,14 @@ bool PolygonXY::contains(const Point2& P) const { }; // check first bounding box - if (coordinatesMax_[LAT] < P[LAT] || P[LAT] < coordinatesMin_[LAT] || coordinatesMax_[LON] < P[LON] || - P[LON] < coordinatesMin_[LON]) { + if (coordinatesMax_.y() < P.y() || P.y() < coordinatesMin_.y() || coordinatesMax_.x() < P.x() || + P.x() < coordinatesMin_.x()) { return false; } if (inner_radius_squared_ == 0) { // check inner bounding box - if (inner_coordinatesMin_[LON] <= P[LON] && P[LON] <= inner_coordinatesMax_[LON] && - inner_coordinatesMin_[LAT] <= P[LAT] && P[LAT] <= inner_coordinatesMax_[LAT]) { + if (inner_coordinatesMin_.x() <= P.x() && P.x() <= inner_coordinatesMax_.x() && + inner_coordinatesMin_.y() <= P.y() && P.y() <= inner_coordinatesMax_.y()) { return true; } } @@ -132,8 +132,8 @@ bool PolygonXY::contains(const Point2& P) const { // intersecting either: // - "up" on upward crossing & P left of edge, or // - "down" on downward crossing & P right of edge - const bool APB = (A[LAT] <= P[LAT] && P[LAT] <= B[LAT]); - const bool BPA = (B[LAT] <= P[LAT] && P[LAT] <= A[LAT]); + const bool APB = (A.y() <= P.y() && P.y() <= B.y()); + const bool BPA = (B.y() <= P.y() && P.y() <= A.y()); if (APB != BPA) { const double side = cross_product_analog(P, A, B); diff --git a/src/atlas/util/PolygonXY.h b/src/atlas/util/PolygonXY.h index ca2489acd..403d913cc 100644 --- a/src/atlas/util/PolygonXY.h +++ b/src/atlas/util/PolygonXY.h @@ -37,10 +37,10 @@ class PolygonXY : public PolygonCoordinates { bool contains(const Point2& Pxy) const override; private: - PointLonLat centroid_; + Point2 centroid_; double inner_radius_squared_{0}; - PointLonLat inner_coordinatesMin_; - PointLonLat inner_coordinatesMax_; + Point2 inner_coordinatesMin_; + Point2 inner_coordinatesMax_; }; diff --git a/src/atlas/util/function/VortexRollup.cc b/src/atlas/util/function/VortexRollup.cc new file mode 100644 index 000000000..a4aa591ac --- /dev/null +++ b/src/atlas/util/function/VortexRollup.cc @@ -0,0 +1,43 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include "atlas/util/Constants.h" +#include "atlas/util/Earth.h" + +#include "atlas/util/function/VortexRollup.h" + +namespace atlas { + +namespace util { + +namespace function { + +double vortex_rollup(double lon, double lat, double t) { + lon *= Constants::degreesToRadians(); + lat *= Constants::degreesToRadians(); + + auto sqr = [](const double x) { return x * x; }; + auto sech = [](const double x) { return 1. / std::cosh(x); }; + constexpr double two_pi = 2. * M_PI; + const double lambda_prime = std::atan2(-std::cos(lon - two_pi * t), std::tan(lat)); + const double rho = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - two_pi * t))); + double omega = 0.; + double a = Earth::radius(); + if (rho != 0.) { + omega = 0.5 * 3 * std::sqrt(3) * a * two_pi * sqr(sech(rho)) * std::tanh(rho) / rho; + } + return -std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t)); +} + +} // namespace function + +} // namespace util + +} // namespace atlas diff --git a/src/atlas/util/function/VortexRollup.h b/src/atlas/util/function/VortexRollup.h new file mode 100644 index 000000000..46f6138f4 --- /dev/null +++ b/src/atlas/util/function/VortexRollup.h @@ -0,0 +1,39 @@ +/* + * (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. + */ + + +#pragma once + +namespace atlas { + +namespace util { + +namespace function { + +/// \brief An analytic function that provides a vortex rollup on a 2D sphere +/// +/// \detailed The formula is found in +/// "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere" +/// by Peter Bosler, James Kent, Robert Krasny, Christiane Jablonowski, JCP 2015 +/// as the tracer density in Test Case 1. +/// The longitude (lon) and latitude (lat) are assumed to be in degrees, +/// The time parameter associated with the vortex rollup is set by (t). +/// +/// The time it takes for the counter-rotating vortices along +/// the equator to return to its original position takes +/// time t = 1.0; +/// +double vortex_rollup(double lon, double lat, double t); + +} // namespace function + +} // namespace util + +} // namespace atlas diff --git a/src/sandbox/interpolation/atlas-parallel-interpolation.cc b/src/sandbox/interpolation/atlas-parallel-interpolation.cc index 7b3aa4527..5a96dcbd2 100644 --- a/src/sandbox/interpolation/atlas-parallel-interpolation.cc +++ b/src/sandbox/interpolation/atlas-parallel-interpolation.cc @@ -10,41 +10,21 @@ #include #include + +#include "eckit/config/Resource.h" +#include "eckit/log/Plural.h" + +#include "PartitionedMesh.h" #include "atlas/field.h" #include "atlas/functionspace.h" #include "atlas/interpolation.h" +#include "atlas/linalg/sparse/Backend.h" #include "atlas/runtime/AtlasTool.h" #include "atlas/runtime/Log.h" -#include "eckit/config/Resource.h" -#include "eckit/linalg/LinearAlgebra.h" -#include "eckit/log/Plural.h" - -#include "PartitionedMesh.h" +#include "atlas/util/function/VortexRollup.h" using namespace atlas; -auto vortex_rollup = [](double lon, double lat, double t) { - // lon and lat in radians! - - // Formula found in "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere" - // by Peter Bosler, James Kent, Robert Krasny, CHristiane Jablonowski, JCP 2015 - - auto sqr = [](const double x) { return x * x; }; - auto sech = [](const double x) { return 1. / std::cosh(x); }; - const double T = 1.; - const double Omega = 2. * M_PI / T; - t *= T; - const double lambda_prime = std::atan2(-std::cos(lon - Omega * t), std::tan(lat)); - const double rho = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - Omega * t))); - double omega = 0.; - double a = util::Earth::radius(); - if (rho != 0.) { - omega = 0.5 * 3 * std::sqrt(3) * a * Omega * sqr(sech(rho)) * std::tanh(rho) / rho; - } - double q = 1. - std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t)); - return q; -}; - class AtlasParallelInterpolation : public AtlasTool { int execute(const AtlasTool::Args& args) override; std::string briefDescription() override { return "Demonstration of parallel interpolation"; } @@ -122,7 +102,7 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { if (args.get("backend", option)) { - eckit::linalg::LinearAlgebra::backend(option); + linalg::sparse::current_backend(option); } // Generate and partition source & target mesh @@ -243,7 +223,7 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { // src_scalar_1( j ) = -std::tanh( y / 10. * std::cos( 50. / std::sqrt( x * x + y * y ) ) - // x / 10. * std::sin( 50. / std::sqrt( x * x + y * y ) ) ); - src_scalar_1(j) = vortex_rollup(lon, lat, 1.); + src_scalar_1(j) = util::function::vortex_rollup(lonlat(j, 0), lonlat(j, 1), 1.); } } diff --git a/src/sandbox/interpolation/atlas-parallel-structured-interpolation.cc b/src/sandbox/interpolation/atlas-parallel-structured-interpolation.cc index da06f7097..2be8b3b07 100644 --- a/src/sandbox/interpolation/atlas-parallel-structured-interpolation.cc +++ b/src/sandbox/interpolation/atlas-parallel-structured-interpolation.cc @@ -19,6 +19,7 @@ #include "atlas/runtime/AtlasTool.h" #include "atlas/runtime/Log.h" #include "atlas/util/CoordinateEnums.h" +#include "atlas/util/function/VortexRollup.h" using namespace atlas; @@ -87,32 +88,6 @@ static Config processed_config(const eckit::Configuration& _config) { return config; } - -double vortex_rollup(double lon, double lat, double t) { - // lon and lat in degrees! - - // Formula found in "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere" - // by Peter Bosler, James Kent, Robert Krasny, Christiane Jablonowski, JCP 2015 - - lon *= M_PI / 180.; - lat *= M_PI / 180.; - - auto sqr = [](const double x) { return x * x; }; - auto sech = [](const double x) { return 1. / std::cosh(x); }; - const double T = 1.; - const double Omega = 2. * M_PI / T; - t *= T; - const double lambda_prime = std::atan2(-std::cos(lon - Omega * t), std::tan(lat)); - const double rho = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - Omega * t))); - double omega = 0.; - double a = util::Earth::radius(); - if (rho != 0.) { - omega = 0.5 * 3 * std::sqrt(3) * a * Omega * sqr(sech(rho)) * std::tanh(rho) / rho; - } - double q = 1. - std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t)); - return q; -}; - int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { ATLAS_TRACE("AtlasParallelInterpolation::execute"); auto source_gridname = args.getString("source", "O32"); @@ -180,7 +155,7 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { else { idx_t k = args.getInt("vortex-rollup", 0); atlas_omp_parallel_for(idx_t n = 0; n < size; ++n) { - source(n) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); + source(n) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); } } } diff --git a/src/tests/AtlasTestEnvironment.h b/src/tests/AtlasTestEnvironment.h index 528810c46..d656f9feb 100644 --- a/src/tests/AtlasTestEnvironment.h +++ b/src/tests/AtlasTestEnvironment.h @@ -106,8 +106,8 @@ Test& current_test() { #undef CASE #define CASE(description) \ void UNIQUE_NAME2(test_, __LINE__)(std::string&, int&, int); \ - static eckit::testing::TestRegister UNIQUE_NAME2(test_registration_, __LINE__)(description, \ - &UNIQUE_NAME2(test_, __LINE__)); \ + static const eckit::testing::TestRegister UNIQUE_NAME2(test_registration_, __LINE__)( \ + description, &UNIQUE_NAME2(test_, __LINE__)); \ void UNIQUE_NAME2(traced_test_, __LINE__)(std::string & _test_subsection, int& _num_subsections, int _subsection); \ void UNIQUE_NAME2(test_, __LINE__)(std::string & _test_subsection, int& _num_subsections, int _subsection) { \ Test UNIQUE_NAME2(testobj_, __LINE__)(description, Here()); \ diff --git a/src/tests/acceptance_tests/atest_mgrids.cc b/src/tests/acceptance_tests/atest_mgrids.cc index b4a452704..e5022f6ad 100644 --- a/src/tests/acceptance_tests/atest_mgrids.cc +++ b/src/tests/acceptance_tests/atest_mgrids.cc @@ -29,11 +29,10 @@ #include "atlas/util/Config.h" #include "atlas/util/Constants.h" #include "atlas/util/CoordinateEnums.h" +#include "atlas/util/function/VortexRollup.h" using namespace atlas; -double vortex_rollup(double lon, double lat, double t, double mean); - //------------------------------------------------------------------------------ class Program : public AtlasTool { @@ -100,7 +99,7 @@ int Program::execute(const Args& args) { double meanA = 1.; for (idx_t n = 0; n < fsA.size(); ++n) { - A(n) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 1., meanA); + A(n) = meanA + util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 1.); } fieldA.set_dirty(); fieldA.haloExchange(); @@ -204,33 +203,6 @@ int Program::execute(const Args& args) { //------------------------------------------------------------------------------ -double vortex_rollup(double lon, double lat, double t, double mean) { - // lon and lat in degrees! - - // Formula found in "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere" - // by Peter Bosler, James Kent, Robert Krasny, CHristiane Jablonowski, JCP 2015 - - lon *= util::Constants::degreesToRadians(); - lat *= util::Constants::degreesToRadians(); - - auto sqr = [](const double x) { return x * x; }; - auto sech = [](const double x) { return 1. / std::cosh(x); }; - const double T = 1.; - const double Omega = 2. * M_PI / T; - t *= T; - const double lambda_prime = std::atan2(-std::cos(lon - Omega * t), std::tan(lat)); - const double rho = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - Omega * t))); - double omega = 0.; - double a = util::Earth::radius(); - if (rho != 0.) { - omega = 0.5 * 3 * std::sqrt(3) * a * Omega * sqr(sech(rho)) * std::tanh(rho) / rho; - } - double q = mean - std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t)); - return q; -}; - -//------------------------------------------------------------------------------ - int main(int argc, char** argv) { Program tool(argc, argv); return tool.start(); diff --git a/src/tests/field/test_field_missingvalue.cc b/src/tests/field/test_field_missingvalue.cc index 99f4d6c2e..b55e82617 100644 --- a/src/tests/field/test_field_missingvalue.cc +++ b/src/tests/field/test_field_missingvalue.cc @@ -123,13 +123,20 @@ CASE("MissingValue (DataType specialisations)") { config.set("missing_value", n); config.set("missing_value_epsilon", eps); + Log::info().indent(); for (std::string type : {"nan", "equals", "approximately-equals"}) { + Log::info() << "type " << type << std::endl; + Log::info().indent(); auto mv = MissingValue(type + "-real64", config); EXPECT(bool(mv)); EXPECT(mv(type == "nan" ? nan : n)); - EXPECT(mv(n) != mv(nan)); + if (type == "nan") { + EXPECT(mv(n) != mv(nan)); + } EXPECT(mv(n + 1) == false); + Log::info().unindent(); } + Log::info().unindent(); } @@ -142,13 +149,20 @@ CASE("MissingValue (DataType specialisations)") { config.set("missing_value", n); config.set("missing_value_epsilon", eps); + Log::info().indent(); for (std::string type : {"nan", "equals", "approximately-equals"}) { + Log::info() << "type " << type << std::endl; + Log::info().indent(); auto mv = MissingValue(type + "-real32", config); EXPECT(bool(mv)); EXPECT(mv(type == "nan" ? nan : n)); - EXPECT(mv(n) != mv(nan)); + if (type == "nan") { + EXPECT(mv(n) != mv(nan)); + } EXPECT(mv(n + 1) == false); + Log::info().unindent(); } + Log::info().unindent(); } @@ -181,7 +195,7 @@ CASE("MissingValue (DataType specialisations)") { CASE("MissingValue from Field (basic)") { - std::vector values{1., nan, missingValue, missingValue, missingValue + missingValueEps / 2., 6., 7.}; + std::vector values{1., missingValue, missingValue, missingValue + missingValueEps / 2., 6., 7.}; Field field("field", values.data(), array::make_shape(values.size(), 1)); field.metadata().set("missing_value_type", "not defined"); @@ -192,12 +206,15 @@ CASE("MissingValue from Field (basic)") { SECTION("nan") { + std::vector values_with_nan = values; + values_with_nan.insert(values_with_nan.begin() + 1, nan); + // missing value type from user - EXPECT(std::count_if(values.begin(), values.end(), MissingValue("nan", field)) == 1); + EXPECT(std::count_if(values_with_nan.begin(), values_with_nan.end(), MissingValue("nan", field)) == 1); // missing value type from field field.metadata().set("missing_value_type", "nan"); - EXPECT(std::count_if(values.begin(), values.end(), MissingValue(field)) == 1); + EXPECT(std::count_if(values_with_nan.begin(), values_with_nan.end(), MissingValue(field)) == 1); } diff --git a/src/tests/functionspace/test_cubedsphere_functionspace.cc b/src/tests/functionspace/test_cubedsphere_functionspace.cc index 070c43d67..3ce49a710 100644 --- a/src/tests/functionspace/test_cubedsphere_functionspace.cc +++ b/src/tests/functionspace/test_cubedsphere_functionspace.cc @@ -18,6 +18,7 @@ #include "atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.h" #include "atlas/option.h" #include "atlas/util/CoordinateEnums.h" +#include "atlas/util/function/VortexRollup.h" #include "tests/AtlasTestEnvironment.h" @@ -27,11 +28,6 @@ namespace test { // Allow small differences in the last few bits of a double aprroximately equal to 1 constexpr double epsilon = std::numeric_limits::epsilon() * 16; - -double testFunction(double lon, double lat) { - return std::sin(3 * lon * M_PI / 180) * std::sin(2 * lat * M_PI / 180); -} - template void testFunctionSpace(const functionspace::CubedSphereColumns& functionspace) { // Make field. @@ -63,7 +59,7 @@ void testFunctionSpace(const functionspace::CubedSphereColumns(functionSpaceA.lonlat()); + const auto lonLatFieldB = array::make_view(functionSpaceB.lonlat()); + const auto ghostFieldA = array::make_view(functionSpaceA.ghost()); + const auto ghostFieldB = array::make_view(functionSpaceB.ghost()); + + // Loop over functionspaces. + for (idx_t i = 0; i < functionSpaceA.size(); ++i) { + EXPECT_EQ(lonLatFieldA(i, LON), lonLatFieldB(i, LON)); + EXPECT_EQ(lonLatFieldA(i, LAT), lonLatFieldB(i, LAT)); + EXPECT_EQ(ghostFieldA(i), ghostFieldB(i)); + } + }; + + // Check that primal cells and dual nodes are equivalent. + compareFields(primalCells, dualNodes); + + // Check that dual cells and primal nodes are equivalent. + compareFields(dualCells, primalNodes); +} + +CASE("Variable halo size functionspaces") { + // Create a mesh with a large halo, and a few functionspaces with different + // (smaller) halo sizes. These should create fields with a smaller memory + // footprint. + + // Set grid. + const auto grid = Grid("CS-LFR-C-12"); + + // Set mesh config. + const auto meshConfig = util::Config("partitioner", "equal_regions") | util::Config("halo", 3); + + // Set mesh. + const auto mesh = MeshGenerator("cubedsphere", meshConfig).generate(grid); + + // Set functionspaces. + const auto nodeColumns0 = functionspace::CubedSphereNodeColumns(mesh, util::Config("halo", 0)); + const auto nodeColumns1 = functionspace::CubedSphereNodeColumns(mesh, util::Config("halo", 1)); + const auto nodeColumns2 = functionspace::CubedSphereNodeColumns(mesh, util::Config("halo", 2)); + + const auto cellColumns0 = functionspace::CubedSphereCellColumns(mesh, util::Config("halo", 0)); + const auto cellColumns1 = functionspace::CubedSphereCellColumns(mesh, util::Config("halo", 1)); + const auto cellColumns2 = functionspace::CubedSphereCellColumns(mesh, util::Config("halo", 2)); + + // Check functionspace sizes. + EXPECT(nodeColumns0.size() < nodeColumns1.size()); + EXPECT(nodeColumns1.size() < nodeColumns2.size()); + EXPECT(nodeColumns2.size() < mesh.nodes().size()); + EXPECT(cellColumns0.size() < cellColumns1.size()); + EXPECT(cellColumns1.size() < cellColumns2.size()); + EXPECT(cellColumns2.size() < mesh.cells().size()); + + // Make sure size of owned cell data matches grid. + auto checkSize = [&](idx_t sizeOwned) { + mpi::comm().allReduceInPlace(sizeOwned, eckit::mpi::Operation::SUM); + EXPECT_EQ(sizeOwned, grid.size()); + }; + + checkSize(cellColumns0.sizeOwned()); + checkSize(cellColumns1.sizeOwned()); + checkSize(cellColumns2.sizeOwned()); +} } // namespace test } // namespace atlas diff --git a/src/tests/functionspace/test_stencil.cc b/src/tests/functionspace/test_stencil.cc index 63e76bff0..6ea6ada9f 100644 --- a/src/tests/functionspace/test_stencil.cc +++ b/src/tests/functionspace/test_stencil.cc @@ -9,7 +9,6 @@ */ #include -#include "eckit/linalg/LinearAlgebra.h" #include "eckit/linalg/Vector.h" #include "eckit/types/Types.h" diff --git a/src/tests/interpolation/test_interpolation_cubic_prototype.cc b/src/tests/interpolation/test_interpolation_cubic_prototype.cc index b5d1f94a5..ff1a3bdf7 100644 --- a/src/tests/interpolation/test_interpolation_cubic_prototype.cc +++ b/src/tests/interpolation/test_interpolation_cubic_prototype.cc @@ -11,10 +11,6 @@ #include #include -#include "eckit/linalg/LinearAlgebra.h" -#include "eckit/linalg/Vector.h" -#include "eckit/types/Types.h" - #include "atlas/array.h" #include "atlas/field/Field.h" #include "atlas/functionspace/NodeColumns.h" @@ -30,6 +26,16 @@ #include "atlas/util/CoordinateEnums.h" #include "atlas/util/MicroDeg.h" +#if ATLAS_ECKIT_HAVE_ECKIT_585 +#include "eckit/linalg/LinearAlgebraSparse.h" +#else +#include "eckit/linalg/LinearAlgebra.h" +#endif +#include "eckit/linalg/SparseMatrix.h" +#include "eckit/linalg/Vector.h" +#include "eckit/types/Types.h" + + #include "CubicInterpolationPrototype.h" #include "tests/AtlasTestEnvironment.h" @@ -233,7 +239,11 @@ CASE("test horizontal cubic interpolation triplets") { std::vector tgt(departure_points.size()); eckit::linalg::Vector v_src(const_cast(f.data()), f.size()); eckit::linalg::Vector v_tgt(tgt.data(), tgt.size()); +#if ATLAS_ECKIT_HAVE_ECKIT_585 + eckit::linalg::LinearAlgebraSparse::backend().spmv(matrix, v_src, v_tgt); +#else eckit::linalg::LinearAlgebra::backend().spmv(matrix, v_src, v_tgt); +#endif Log::info() << "output = " << tgt << std::endl; } diff --git a/src/tests/interpolation/test_interpolation_structured2D.cc b/src/tests/interpolation/test_interpolation_structured2D.cc index 266122c95..be4c13ad0 100644 --- a/src/tests/interpolation/test_interpolation_structured2D.cc +++ b/src/tests/interpolation/test_interpolation_structured2D.cc @@ -21,6 +21,7 @@ #include "atlas/meshgenerator.h" #include "atlas/output/Gmsh.h" #include "atlas/util/CoordinateEnums.h" +#include "atlas/util/function/VortexRollup.h" #include "tests/AtlasTestEnvironment.h" @@ -116,31 +117,6 @@ FunctionSpace output_functionspace(const Grid& grid) { return NodeColumns{output_mesh}; } -double vortex_rollup(double lon, double lat, double t) { - // lon and lat in degrees! - - // Formula found in "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere" - // by Peter Bosler, James Kent, Robert Krasny, CHristiane Jablonowski, JCP 2015 - - lon *= M_PI / 180.; - lat *= M_PI / 180.; - - auto sqr = [](const double x) { return x * x; }; - auto sech = [](const double x) { return 1. / std::cosh(x); }; - const double T = 1.; - const double Omega = 2. * M_PI / T; - t *= T; - const double lambda_prime = std::atan2(-std::cos(lon - Omega * t), std::tan(lat)); - const double rho = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - Omega * t))); - double omega = 0.; - double a = util::Earth::radius(); - if (rho != 0.) { - omega = 0.5 * 3 * std::sqrt(3) * a * Omega * sqr(sech(rho)) * std::tanh(rho) / rho; - } - double q = 1. - std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t)); - return q; -}; - CASE("which scheme?") { Log::info() << scheme().getString("type") << std::endl; } @@ -162,7 +138,7 @@ CASE("test_interpolation_structured using functionspace API") { auto lonlat = array::make_view(input_fs.xy()); auto source = array::make_view(field_source); for (idx_t n = 0; n < input_fs.size(); ++n) { - source(n) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 1.); + source(n) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 1.); } EXPECT(field_source.dirty()); @@ -202,7 +178,7 @@ CASE("test_interpolation_structured using grid API") { idx_t n{0}; for (auto p : input_grid.lonlat()) { - src_data[n++] = vortex_rollup(p.lon(), p.lat(), 1.); + src_data[n++] = util::function::vortex_rollup(p.lon(), p.lat(), 1.); } // Wrap memory in atlas Fields and interpolate @@ -259,7 +235,7 @@ CASE("test_interpolation_structured using fs API multiple levels") { auto source = array::make_view(field_source); for (idx_t n = 0; n < input_fs.size(); ++n) { for (idx_t k = 0; k < 3; ++k) { - source(n, k) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); + source(n, k) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); } } interpolation.execute(field_source, field_target); @@ -281,7 +257,7 @@ struct AdjointTolerance { template <> const double AdjointTolerance::value = 2.e-14; template <> -const float AdjointTolerance::value = 4.e-6; +const float AdjointTolerance::value = 2.e-5; template @@ -307,7 +283,7 @@ void test_interpolation_structured_using_fs_API_for_fieldset() { auto source = array::make_view(field_source); for (idx_t n = 0; n < input_fs.size(); ++n) { for (idx_t k = 0; k < 3; ++k) { - source(n, k) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); + source(n, k) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); } } } @@ -483,7 +459,7 @@ CASE("ATLAS-315: Target grid with domain West of 0 degrees Lon") { auto source = array::make_view(field_src); constexpr double k = 1; for (idx_t n = 0; n < source.size(); ++n) { - source(n) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); + source(n) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); } interpolation.execute(field_src, field_tgt); diff --git a/src/tests/interpolation/test_interpolation_structured2D_to_unstructured.cc b/src/tests/interpolation/test_interpolation_structured2D_to_unstructured.cc index 25d0704dc..65ac5e42b 100644 --- a/src/tests/interpolation/test_interpolation_structured2D_to_unstructured.cc +++ b/src/tests/interpolation/test_interpolation_structured2D_to_unstructured.cc @@ -21,6 +21,7 @@ #include "atlas/meshgenerator.h" #include "atlas/output/Gmsh.h" #include "atlas/util/CoordinateEnums.h" +#include "atlas/util/function/VortexRollup.h" #include "tests/AtlasTestEnvironment.h" @@ -37,32 +38,6 @@ std::string input_gridname(const std::string& default_grid) { return eckit::Resource("--input-grid", default_grid); } -double vortex_rollup(double lon, double lat, double t) { - // lon and lat in degrees! - - // Formula found in "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere" - // by Peter Bosler, James Kent, Robert Krasny, CHristiane Jablonowski, JCP 2015 - - lon *= M_PI / 180.; - lat *= M_PI / 180.; - - auto sqr = [](const double x) { return x * x; }; - auto sech = [](const double x) { return 1. / std::cosh(x); }; - const double T = 1.; - const double Omega = 2. * M_PI / T; - t *= T; - const double lambda_prime = std::atan2(-std::cos(lon - Omega * t), std::tan(lat)); - const double rho = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - Omega * t))); - double omega = 0.; - double a = util::Earth::radius(); - if (rho != 0.) { - omega = 0.5 * 3 * std::sqrt(3) * a * Omega * sqr(sech(rho)) * std::tanh(rho) / rho; - } - double q = 1. - std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t)); - return q; -}; - - FunctionSpace output_functionspace_match() { std::vector points; if (mpi::size() == 2) { @@ -124,7 +99,7 @@ FieldSet create_source_fields(StructuredColumns& fs, idx_t nb_fields, idx_t nb_l auto source = array::make_view(field_source); for (idx_t n = 0; n < fs.size(); ++n) { for (idx_t k = 0; k < nb_levels; ++k) { - source(n, k) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); + source(n, k) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2); } }; } diff --git a/src/tests/linalg/test_linalg_dense.cc b/src/tests/linalg/test_linalg_dense.cc index 687880798..a085fe5a9 100644 --- a/src/tests/linalg/test_linalg_dense.cc +++ b/src/tests/linalg/test_linalg_dense.cc @@ -11,6 +11,13 @@ #include #include +#include "atlas/library/config.h" +#if ATLAS_ECKIT_HAVE_ECKIT_585 +#include "eckit/linalg/LinearAlgebraDense.h" +#else +#include "eckit/linalg/LinearAlgebra.h" +#endif + #include "eckit/linalg/Matrix.h" #include "eckit/linalg/Vector.h" @@ -67,7 +74,11 @@ void expect_equal(const T1& v, const T2& r) { CASE("test configuration via resource") { if (atlas::Library::instance().linalgDenseBackend().empty()) { +#if ATLAS_ECKIT_HAVE_ECKIT_585 + if (eckit::linalg::LinearAlgebraDense::hasBackend("mkl")) { +#else if (eckit::linalg::LinearAlgebra::hasBackend("mkl")) { +#endif EXPECT_EQ(dense::current_backend().type(), "mkl"); } else { diff --git a/src/tests/linalg/test_linalg_sparse.cc b/src/tests/linalg/test_linalg_sparse.cc index 181341c5d..ee1b10302 100644 --- a/src/tests/linalg/test_linalg_sparse.cc +++ b/src/tests/linalg/test_linalg_sparse.cc @@ -14,6 +14,7 @@ #include "eckit/linalg/Matrix.h" #include "eckit/linalg/Vector.h" +#include "atlas/array.h" #include "atlas/linalg/sparse.h" #include "tests/AtlasTestEnvironment.h" diff --git a/src/tests/mesh/test_cubedsphere_meshgen.cc b/src/tests/mesh/test_cubedsphere_meshgen.cc index 533733436..fef7770ad 100644 --- a/src/tests/mesh/test_cubedsphere_meshgen.cc +++ b/src/tests/mesh/test_cubedsphere_meshgen.cc @@ -22,6 +22,7 @@ #include "atlas/option.h" #include "atlas/output/Gmsh.h" #include "atlas/util/CoordinateEnums.h" +#include "atlas/util/function/VortexRollup.h" #include "tests/AtlasTestEnvironment.h" #include "eckit/mpi/Operation.h" @@ -119,33 +120,6 @@ CASE("cubedsphere_mesh_jacobian_test") { } } - -double testFunction(double lon, double lat, double t = 1.) { - // lon and lat in degrees! - - // Formula found in "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere" - // by Peter Bosler, James Kent, Robert Krasny, CHristiane Jablonowski, JCP 2015 - - lon *= M_PI / 180.; - lat *= M_PI / 180.; - - auto sqr = [](const double x) { return x * x; }; - auto sech = [](const double x) { return 1. / std::cosh(x); }; - const double T = 1.; - const double Omega = 2. * M_PI / T; - t *= T; - const double lambda_prime = std::atan2(-std::cos(lon - Omega * t), std::tan(lat)); - const double rho = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - Omega * t))); - double omega = 0.; - double a = util::Earth::radius(); - if (rho != 0.) { - omega = 0.5 * 3 * std::sqrt(3) * a * Omega * sqr(sech(rho)) * std::tanh(rho) / rho; - } - double q = 1. - std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t)); - return q; -}; - - void testHaloExchange(const std::string& gridStr, const std::string& partitionerStr, idx_t halo, bool output = true) { // Set grid. const auto grid = Grid(gridStr); @@ -185,7 +159,7 @@ void testHaloExchange(const std::string& gridStr, const std::string& partitioner break; } - testView1(i) = testFunction(lonLatView(i, LON), lonLatView(i, LAT)); + testView1(i) = 1.0 + util::function::vortex_rollup(lonLatView(i, LON), lonLatView(i, LAT), 1.0); ++testFuncCallCount; } @@ -199,7 +173,7 @@ void testHaloExchange(const std::string& gridStr, const std::string& partitioner // Check all values after halo exchange. double maxError = 0; for (idx_t i = 0; i < nodeColumns.size(); ++i) { - const double testVal = testFunction(lonLatView(i, LON), lonLatView(i, LAT)); + const double testVal = 1.0 + util::function::vortex_rollup(lonLatView(i, LON), lonLatView(i, LAT), 1.0); maxError = std::max(maxError, std::abs(testView1(i) - testVal)); } @@ -230,7 +204,7 @@ void testHaloExchange(const std::string& gridStr, const std::string& partitioner break; } - testView2(i) = testFunction(lonLatView(i, LON), lonLatView(i, LAT)); + testView2(i) = 1.0 + util::function::vortex_rollup(lonLatView(i, LON), lonLatView(i, LAT), 1.0); ++testFuncCallCount; } @@ -245,7 +219,7 @@ void testHaloExchange(const std::string& gridStr, const std::string& partitioner maxError = 0; for (idx_t i = 0; i < cellColumns.size(); ++i) { // Test field and test function should be the same. - const double testVal = testFunction(lonLatView(i, LON), lonLatView(i, LAT)); + const double testVal = 1.0 + util::function::vortex_rollup(lonLatView(i, LON), lonLatView(i, LAT), 1.0); maxError = std::max(maxError, std::abs(testView2(i) - testVal)); } @@ -357,8 +331,9 @@ CASE("cubedsphere_dual_mesh_test") { for (idx_t i = 0; i < sourceView.shape(0); ++i) { for (idx_t j = 0; j < sourceView.shape(1); ++j) { - sourceView(i, j) = testFunction(sourceLonLatView(i, LON), sourceLonLatView(i, LAT), - static_cast(j) / (nLevels - 1)); + sourceView(i, j) = + 1.0 + util::function::vortex_rollup(sourceLonLatView(i, LON), sourceLonLatView(i, LAT), + static_cast(j) / (nLevels - 1)); } } @@ -384,8 +359,9 @@ CASE("cubedsphere_dual_mesh_test") { double maxError = 0.; for (idx_t i = 0; i < targetView.shape(0); ++i) { for (idx_t j = 0; j < targetView.shape(1); ++j) { - double referenceVal = testFunction(targetLonLatView(i, LON), targetLonLatView(i, LAT), - static_cast(j) / (nLevels - 1)); + double referenceVal = + 1.0 + util::function::vortex_rollup(targetLonLatView(i, LON), targetLonLatView(i, LAT), + static_cast(j) / (nLevels - 1)); double relativeError = std::abs((targetView(i, j) - referenceVal) / referenceVal); maxError = std::max(maxError, relativeError); @@ -438,8 +414,8 @@ CASE("cubedsphere_dual_mesh_test") { } for (idx_t j = 0; j < fieldView.shape(1); ++j) { - fieldView(i, j) = - testFunction(lonLatView(i, LON), lonLatView(i, LAT), static_cast(j) / (nLevels - 1)); + fieldView(i, j) = 1.0 + util::function::vortex_rollup(lonLatView(i, LON), lonLatView(i, LAT), + static_cast(j) / (nLevels - 1)); } } @@ -449,8 +425,9 @@ CASE("cubedsphere_dual_mesh_test") { // Check test field on *all* cells. for (idx_t i = 0; i < fieldView.shape(0); ++i) { for (idx_t j = 0; j < fieldView.shape(1); ++j) { - EXPECT_APPROX_EQ(fieldView(i, j), testFunction(lonLatView(i, LON), lonLatView(i, LAT), - static_cast(j) / (nLevels - 1))); + EXPECT_APPROX_EQ(fieldView(i, j), + 1.0 + util::function::vortex_rollup(lonLatView(i, LON), lonLatView(i, LAT), + static_cast(j) / (nLevels - 1))); } } diff --git a/src/tests/mesh/test_healpixmeshgen.cc b/src/tests/mesh/test_healpixmeshgen.cc index 695593880..c83be6aaa 100644 --- a/src/tests/mesh/test_healpixmeshgen.cc +++ b/src/tests/mesh/test_healpixmeshgen.cc @@ -14,6 +14,7 @@ #include "atlas/meshgenerator.h" #include "atlas/meshgenerator/detail/HealpixMeshGenerator.h" #include "atlas/output/Gmsh.h" +#include "atlas/util/PolygonXY.h" #include "tests/AtlasTestEnvironment.h" @@ -131,6 +132,33 @@ CASE("construction by integer") { EXPECT(HealpixGrid(3) == Grid("H3")); } +//----------------------------------------------------------------------------- + +CASE("matching mesh partitioner") { + auto grid = Grid{"H8"}; + auto mesh = HealpixMeshGenerator{}.generate(grid); + auto match = MatchingPartitioner{mesh}; + auto& polygons = mesh.polygons(); + + static bool do_once = [&]() { + for (idx_t i = 0; i < polygons.size(); ++i) { + auto poly = util::PolygonXY{polygons[i]}; + Log::info() << "polygon[" << i << "]:\n"; + for (idx_t j = 0; j < poly.size(); ++j) { + Log::info() << " " << std::setw(5) << std::left << j << std::setprecision(16) << poly[j] << std::endl; + } + } + return true; + }(); + + + SECTION("H8 -> O64") { EXPECT_NO_THROW(match.partition(Grid{"O64"})); } + SECTION("H8 -> L32x17") { EXPECT_NO_THROW(match.partition(Grid{"L32x17"})); } + SECTION("H8 -> S32x17") { EXPECT_NO_THROW(match.partition(Grid{"S32x17"})); } + SECTION("H8 -> F32") { EXPECT_NO_THROW(match.partition(Grid{"F32"})); } + SECTION("H8 -> L64x33 (west=-180)") { EXPECT_NO_THROW(match.partition(Grid{"L64x33", GlobalDomain(-180.)})); } +} + } // namespace test } // namespace atlas diff --git a/src/tests/mesh/test_mesh_build_edges.cc b/src/tests/mesh/test_mesh_build_edges.cc index 2055af1ca..97da39bae 100644 --- a/src/tests/mesh/test_mesh_build_edges.cc +++ b/src/tests/mesh/test_mesh_build_edges.cc @@ -247,7 +247,7 @@ CASE("test_build_edges") { // Accumulate facets of cells ( edges in 2D ) mesh::actions::build_edges(mesh, option::pole_edges(false)); - idx_t edge_nodes_check[] = { + std::vector edge_nodes_check{ 0, 21, 21, 22, 22, 1, 1, 0, 22, 23, 23, 2, 2, 1, 3, 25, 25, 26, 26, 4, 4, 3, 26, 27, 27, 5, 5, 4, 27, 28, 28, 6, 6, 5, 28, 29, 29, 7, 7, 6, 8, 31, 31, 32, 32, 9, 9, 8, 32, 33, 33, 10, 10, 9, 33, 34, 34, 11, 11, 10, 34, 35, 35, 12, 12, 11, 13, 37, 37, 38, 38, 14, 14, 13, 38, 39, 39, 15, 15, 14, 39, @@ -266,206 +266,208 @@ CASE("test_build_edges") { const mesh::HybridElements::Connectivity& edge_node_connectivity = mesh.edges().node_connectivity(); EXPECT(mesh.projection().units() == "degrees"); const util::UniqueLonLat compute_uid(mesh); + EXPECT_EQ(mesh.edges().size(), edge_nodes_check.size() / 2); for (idx_t jedge = 0; jedge < mesh.edges().size(); ++jedge) { if (compute_uid(edge_nodes_check[2 * jedge + 0]) < compute_uid(edge_nodes_check[2 * jedge + 1])) { - EXPECT(edge_nodes_check[2 * jedge + 0] == edge_node_connectivity(jedge, 0)); - EXPECT(edge_nodes_check[2 * jedge + 1] == edge_node_connectivity(jedge, 1)); + EXPECT_EQ(edge_nodes_check[2 * jedge + 0], edge_node_connectivity(jedge, 0)); + EXPECT_EQ(edge_nodes_check[2 * jedge + 1], edge_node_connectivity(jedge, 1)); } else { - EXPECT(edge_nodes_check[2 * jedge + 0] == edge_node_connectivity(jedge, 1)); - EXPECT(edge_nodes_check[2 * jedge + 1] == edge_node_connectivity(jedge, 0)); + EXPECT_EQ(edge_nodes_check[2 * jedge + 0], edge_node_connectivity(jedge, 1)); + EXPECT_EQ(edge_nodes_check[2 * jedge + 1], edge_node_connectivity(jedge, 0)); } } } - idx_t edge_to_cell_check[] = {0, missing_value, - 0, 16, - 0, 1, - 0, missing_value, - 1, 17, - 1, 56, - 1, missing_value, - 2, 58, - 2, 20, - 2, 3, - 2, missing_value, - 3, 21, - 3, 4, - 3, missing_value, - 4, 22, - 4, 5, - 4, missing_value, - 5, 23, - 5, 59, - 5, missing_value, - 6, 61, - 6, 26, - 6, 7, - 6, missing_value, - 7, 27, - 7, 8, - 7, missing_value, - 8, 28, - 8, 9, - 8, missing_value, - 9, 29, - 9, 62, - 9, missing_value, - 10, 64, - 10, 32, - 10, 11, - 10, missing_value, - 11, 33, - 11, 12, - 11, missing_value, - 12, 34, - 12, 13, - 12, missing_value, - 13, 35, - 13, 65, - 13, missing_value, - 14, 67, - 14, 38, - 14, 15, - 14, missing_value, - 15, 39, - 15, missing_value, - 15, missing_value, - 16, missing_value, - 16, 40, - 16, 17, - 17, 41, - 17, 18, - 18, 68, - 18, 19, - 18, 56, - 19, 70, - 19, 20, - 19, 58, - 20, 42, - 20, 21, - 21, 43, - 21, 22, - 22, 44, - 22, 23, - 23, 45, - 23, 24, - 24, 71, - 24, 25, - 24, 59, - 25, 73, - 25, 26, - 25, 61, - 26, 46, - 26, 27, - 27, 47, - 27, 28, - 28, 48, - 28, 29, - 29, 49, - 29, 30, - 30, 74, - 30, 31, - 30, 62, - 31, 76, - 31, 32, - 31, 64, - 32, 50, - 32, 33, - 33, 51, - 33, 34, - 34, 52, - 34, 35, - 35, 53, - 35, 36, - 36, 77, - 36, 37, - 36, 65, - 37, 79, - 37, 38, - 37, 67, - 38, 54, - 38, 39, - 39, 55, - 39, missing_value, - 40, missing_value, - 40, missing_value, - 40, 41, - 41, missing_value, - 41, 68, - 42, 70, - 42, missing_value, - 42, 43, - 43, missing_value, - 43, 44, - 44, missing_value, - 44, 45, - 45, missing_value, - 45, 71, - 46, 73, - 46, missing_value, - 46, 47, - 47, missing_value, - 47, 48, - 48, missing_value, - 48, 49, - 49, missing_value, - 49, 74, - 50, 76, - 50, missing_value, - 50, 51, - 51, missing_value, - 51, 52, - 52, missing_value, - 52, 53, - 53, missing_value, - 53, 77, - 54, 79, - 54, missing_value, - 54, 55, - 55, missing_value, - 55, missing_value, - 56, 57, - 57, 58, - 57, missing_value, - 59, 60, - 60, 61, - 60, missing_value, - 62, 63, - 63, 64, - 63, missing_value, - 65, 66, - 66, 67, - 66, missing_value, - 68, 69, - 69, missing_value, - 69, 70, - 71, 72, - 72, missing_value, - 72, 73, - 74, 75, - 75, missing_value, - 75, 76, - 77, 78, - 78, missing_value, - 78, 79}; + std::vector edge_to_cell_check{0, missing_value, + 0, 16, + 0, 1, + 0, missing_value, + 1, 17, + 1, 56, + 1, missing_value, + 2, 58, + 2, 20, + 2, 3, + 2, missing_value, + 3, 21, + 3, 4, + 3, missing_value, + 4, 22, + 4, 5, + 4, missing_value, + 5, 23, + 5, 59, + 5, missing_value, + 6, 61, + 6, 26, + 6, 7, + 6, missing_value, + 7, 27, + 7, 8, + 7, missing_value, + 8, 28, + 8, 9, + 8, missing_value, + 9, 29, + 9, 62, + 9, missing_value, + 10, 64, + 10, 32, + 10, 11, + 10, missing_value, + 11, 33, + 11, 12, + 11, missing_value, + 12, 34, + 12, 13, + 12, missing_value, + 13, 35, + 13, 65, + 13, missing_value, + 14, 67, + 14, 38, + 14, 15, + 14, missing_value, + 15, 39, + 15, missing_value, + 15, missing_value, + 16, missing_value, + 16, 40, + 16, 17, + 17, 41, + 17, 18, + 18, 68, + 18, 19, + 18, 56, + 19, 70, + 19, 20, + 19, 58, + 20, 42, + 20, 21, + 21, 43, + 21, 22, + 22, 44, + 22, 23, + 23, 45, + 23, 24, + 24, 71, + 24, 25, + 24, 59, + 25, 73, + 25, 26, + 25, 61, + 26, 46, + 26, 27, + 27, 47, + 27, 28, + 28, 48, + 28, 29, + 29, 49, + 29, 30, + 30, 74, + 30, 31, + 30, 62, + 31, 76, + 31, 32, + 31, 64, + 32, 50, + 32, 33, + 33, 51, + 33, 34, + 34, 52, + 34, 35, + 35, 53, + 35, 36, + 36, 77, + 36, 37, + 36, 65, + 37, 79, + 37, 38, + 37, 67, + 38, 54, + 38, 39, + 39, 55, + 39, missing_value, + 40, missing_value, + 40, missing_value, + 40, 41, + 41, missing_value, + 41, 68, + 42, 70, + 42, missing_value, + 42, 43, + 43, missing_value, + 43, 44, + 44, missing_value, + 44, 45, + 45, missing_value, + 45, 71, + 46, 73, + 46, missing_value, + 46, 47, + 47, missing_value, + 47, 48, + 48, missing_value, + 48, 49, + 49, missing_value, + 49, 74, + 50, 76, + 50, missing_value, + 50, 51, + 51, missing_value, + 51, 52, + 52, missing_value, + 52, 53, + 53, missing_value, + 53, 77, + 54, 79, + 54, missing_value, + 54, 55, + 55, missing_value, + 55, missing_value, + 56, 57, + 57, 58, + 57, missing_value, + 59, 60, + 60, 61, + 60, missing_value, + 62, 63, + 63, 64, + 63, missing_value, + 65, 66, + 66, 67, + 66, missing_value, + 68, 69, + 69, missing_value, + 69, 70, + 71, 72, + 72, missing_value, + 72, 73, + 74, 75, + 75, missing_value, + 75, 76, + 77, 78, + 78, missing_value, + 78, 79}; { const mesh::HybridElements::Connectivity& cell_node_connectivity = mesh.cells().node_connectivity(); const mesh::HybridElements::Connectivity& edge_cell_connectivity = mesh.edges().cell_connectivity(); const util::UniqueLonLat compute_uid(mesh); + EXPECT_EQ(mesh.edges().size(), edge_to_cell_check.size() / 2); for (idx_t jedge = 0; jedge < mesh.edges().size(); ++jedge) { idx_t e1 = edge_to_cell_check[2 * jedge + 0]; idx_t e2 = edge_to_cell_check[2 * jedge + 1]; if (e2 == edge_cell_connectivity.missing_value() || compute_uid(cell_node_connectivity.row(e1)) < compute_uid(cell_node_connectivity.row(e2))) { - EXPECT(edge_to_cell_check[2 * jedge + 0] == edge_cell_connectivity(jedge, 0)); - EXPECT(edge_to_cell_check[2 * jedge + 1] == edge_cell_connectivity(jedge, 1)); + EXPECT_EQ(edge_to_cell_check[2 * jedge + 0], edge_cell_connectivity(jedge, 0)); + EXPECT_EQ(edge_to_cell_check[2 * jedge + 1], edge_cell_connectivity(jedge, 1)); } else { std::cout << "jedge " << jedge << std::endl; - EXPECT(edge_to_cell_check[2 * jedge + 0] == edge_cell_connectivity(jedge, 1)); - EXPECT(edge_to_cell_check[2 * jedge + 1] == edge_cell_connectivity(jedge, 0)); + EXPECT_EQ(edge_to_cell_check[2 * jedge + 0], edge_cell_connectivity(jedge, 1)); + EXPECT_EQ(edge_to_cell_check[2 * jedge + 1], edge_cell_connectivity(jedge, 0)); } } } @@ -473,9 +475,9 @@ CASE("test_build_edges") { { const MultiBlockConnectivity& elem_edge_connectivity = mesh.cells().edge_connectivity(); for (idx_t jelem = 0; jelem < mesh.cells().size(); ++jelem) { - std::cout << jelem << " : "; + std::cout << std::setw(3) << jelem << " : "; for (idx_t jedge = 0; jedge < elem_edge_connectivity.cols(jelem); ++jedge) { - std::cout << elem_edge_connectivity(jelem, jedge) << " "; + std::cout << std::setw(3) << elem_edge_connectivity(jelem, jedge) << " "; } std::cout << std::endl; } @@ -494,13 +496,13 @@ CASE("test_build_edges_triangles_only") { const MultiBlockConnectivity& elem_edge_connectivity = mesh.cells().edge_connectivity(); const MultiBlockConnectivity& elem_node_connectivity = mesh.cells().node_connectivity(); for (idx_t jelem = 0; jelem < mesh.cells().size(); ++jelem) { - std::cout << jelem << " : edges ( "; + std::cout << "elem" << std::setw(3) << jelem << " : edges ( "; for (idx_t jedge = 0; jedge < elem_edge_connectivity.cols(jelem); ++jedge) { - std::cout << elem_edge_connectivity(jelem, jedge) << " "; + std::cout << std::setw(3) << elem_edge_connectivity(jelem, jedge) << " "; } - std::cout << ") | nodes ( "; + std::cout << ") | nodes ( "; for (idx_t jnode = 0; jnode < elem_node_connectivity.cols(jelem); ++jnode) { - std::cout << elem_node_connectivity(jelem, jnode) << " "; + std::cout << std::setw(3) << elem_node_connectivity(jelem, jnode) << " "; } std::cout << ")" << std::endl; } diff --git a/src/tests/projection/test_projection_variable_resolution.cc b/src/tests/projection/test_projection_variable_resolution.cc index e1ab5e4e2..f08147342 100644 --- a/src/tests/projection/test_projection_variable_resolution.cc +++ b/src/tests/projection/test_projection_variable_resolution.cc @@ -64,7 +64,7 @@ const std::vector lat_LAM_reg = {-8.264495551724226, -7.613347275862 9.31650789655281919, 9.967656172414931959, 10.61880444827704473, 11.269952724139157, 11.92110100000127}; -const int nx = 34; // How do you find this? +const int nx = 34; // How do you find this? Now in a test using x/yrange_outer and delta_inner const int ny = 32; const double xrange_outer[2] = {348.13120344827576, 369.6190965517242}; const double yrange_outer[2] = {-8.264495551724226, 11.92110100000127}; @@ -75,6 +75,10 @@ const double rim_width = 4.; const double delta_outer = 1.; const double delta_inner = 0.6511482758621128; +///< correction used +constexpr float epstest = std::numeric_limits::epsilon(); + + auto make_var_ratio_projection = [](double var_ratio) { Config conf; conf.set("type", "variable_resolution"); @@ -117,6 +121,124 @@ auto make_var_ratio_projection_rot = [](double var_ratio, std::vector no auto not_equal = [](double a, double b) { return std::abs(b - a) > 1.e-5; }; +static double new_ratio(int n_stretched, double var_ratio) { + /** + * compute ratio, + * change stretching factor so that high and low grids + * retain original sizes + */ + + + ///< number of variable (stretched) grid points in one side + double var_ints_f = n_stretched * 1.; + double logr = std::log(var_ratio); + double log_ratio = (var_ints_f - 0.5) * logr; + + return std::exp(log_ratio / n_stretched); +}; + +/** + * n_stretched and n_rim. Outside the regular grid, only one side, + * number of grid points in the stretched grid and in the rim region +*/ +auto create_stretched_grid = [](const int& n_points_, const int& n_stretched_, const double& var_ratio_, + const int& n_rim, double lamphi, const bool& L_long) { + double new_ratio_{}; + double range_outer[2]{}; + + auto normalised = [L_long](double p) { + if (L_long) { + p = (p < 180) ? p + 360.0 : p; + } + return p; + }; + + lamphi = normalised(lamphi); + + if (L_long) { + range_outer[0] = xrange_outer[0]; + range_outer[1] = xrange_outer[1]; + } + else { + range_outer[0] = yrange_outer[0]; + range_outer[1] = yrange_outer[1]; + } + + if (var_ratio_ > 1.) { + //< compute number of points in the original regular grid + // THIS IS TO NORMALIZE + int n_idx_ = ((lamphi + epstest - range_outer[0]) / delta_inner) + 1; + // first point in the internal regular grid + int nstart_int = n_rim + n_stretched_ + 1; + // last point in the internal regular grid + int nend_int = n_points_ - (n_rim + n_stretched_); + + + double delta = delta_inner; + double delta_last = delta_inner; + double delta_add{}; + double delta_tot{}; + + new_ratio_ = new_ratio(n_stretched_, var_ratio_); + + //< stretched area + //< n_idx_ start from 1 + if (((n_idx_ < nstart_int) && (n_idx_ > n_rim)) || ((n_idx_ > nend_int) && (n_idx_ < n_points_ - n_rim + 1))) { + //< number of stretched points for the loop + int n_st{}; + if (n_idx_ < nstart_int) { + n_st = nstart_int - n_idx_; + } + else { + n_st = n_idx_ - nend_int; + } + + delta_tot = 0.; + for (int ix = 0; ix < n_st; ix += 1) { + delta_last = delta * new_ratio_; + delta_add = delta_last - delta_inner; + delta = delta_last; + delta_tot += delta_add; + } + + if (n_idx_ < nstart_int) { + lamphi -= delta_tot; + } + else { + lamphi += delta_tot; + } + } + + //< rim region + if (((n_idx_ < n_rim + 1)) || (n_idx_ > n_points_ - n_rim)) { + delta_tot = 0.; + //compute total stretched + for (int ix = 0; ix < n_stretched_; ix += 1) { + delta_last = delta * new_ratio_; + delta_add = delta_last - delta_inner; + delta = delta_last; + delta_tot += delta_add; + } + + + double drim_size = ((rim_width / 2.) / n_rim) - delta_inner; + int ndrim{}; + + if (n_idx_ < nstart_int) { + ndrim = nstart_int - n_stretched_ - n_idx_; + lamphi = lamphi - delta_tot - (ndrim * drim_size); + } + else { + ndrim = n_idx_ - (n_points_ - n_rim); + lamphi = lamphi + delta_tot + ndrim * drim_size; + } + } + } + + //< return the new point in the stretched grid + return normalised(lamphi); +}; + }; // namespace @@ -156,9 +278,6 @@ CASE("Understanding of the above data") { int nx_outer = 4; double var_ratio = 1.13; - constexpr double epstest = - std::numeric_limits::epsilon(); ///< correction used to change from double to integer - int inner_i_begin = -1; int inner_i_end = -1; for (int i = 0; i < nx; ++i) { @@ -206,12 +325,70 @@ CASE("Understanding of the above data") { EXPECT_APPROX_EQ(new_ratio, 1.10722, 1.e-5); } +//< create stretched grid from regular grid + +CASE("var_ratio_create = 1.13") { + ///< definition of grid + constexpr float epstest = std::numeric_limits::epsilon(); + int nx_reg = ((xrange_outer[1] - xrange_outer[0] + epstest) / delta_inner) + 1; + int ny_reg = ((yrange_outer[1] - yrange_outer[0] + epstest) / delta_inner) + 1; + + + auto grid_st = + RegularGrid{grid::LinearSpacing{xrange_outer[0], xrange_outer[1], nx_reg}, + grid::LinearSpacing{yrange_outer[0], yrange_outer[1], ny_reg}, make_var_ratio_projection(1.13)}; + + + ///< check over regular grid points stretched using new atlas object and check using look-up table + for (idx_t j = 0; j < grid_st.ny(); ++j) { + for (idx_t i = 0; i < grid_st.nx(); ++i) { + auto ll = grid_st.lonlat(i, j); + EXPECT_APPROX_EQ(ll.lon(), lon_LAM_str[i], 1.e-10); + EXPECT_APPROX_EQ(ll.lat(), lat_LAM_str[j], 1.e-10); + } + } + + + ///< check over regular grid points stretched using new atlas object and check using function above + for (idx_t j = 0; j < grid_st.ny(); ++j) { + for (idx_t i = 0; i < grid_st.nx(); ++i) { + auto ll_st_lon = create_stretched_grid(nx_reg, 3, 1.13, 2, lon_LAM_reg[i], true); + auto ll_st_lat = create_stretched_grid(ny_reg, 3, 1.13, 2, lat_LAM_reg[j], false); + EXPECT_APPROX_EQ(ll_st_lon, lon_LAM_str[i], 1.e-10); + EXPECT_APPROX_EQ(ll_st_lat, lat_LAM_str[j], 1.e-10); + } + } + + ///< check internal regular grid + idx_t ymid = grid_st.ny() / 2; + idx_t xmid = grid_st.nx() / 2; + auto expect_equal_dlon = [&](int i, double dlon) { + EXPECT_APPROX_EQ(grid_st.lonlat(i + 1, ymid).lon() - grid_st.lonlat(i, ymid).lon(), dlon, 1.e-10); + }; + auto expect_equal_dlat = [&](int j, double dlat) { + EXPECT_APPROX_EQ(grid_st.lonlat(xmid, j + 1).lat() - grid_st.lonlat(xmid, j).lat(), dlat, 1.e-10); + }; + expect_equal_dlon(0, delta_outer); + expect_equal_dlon(xmid, delta_inner); + expect_equal_dlat(0, delta_outer); + expect_equal_dlat(ymid, delta_inner); + + // Check that the spacing in xy coordinates matches "delta_inner" + for (int i = 0; i < grid_st.nx() - 1; ++i) { + EXPECT_APPROX_EQ(grid_st.xy(i + 1, ymid).x() - grid_st.xy(i, ymid).x(), delta_inner, 1.e-10); + } + for (int j = 0; j < grid_st.ny() - 1; ++j) { + EXPECT_APPROX_EQ(grid_st.xy(xmid, j + 1).y() - grid_st.xy(xmid, j).y(), delta_inner, 1.e-10); + } +} + CASE("var_ratio = 1.13") { ///< definition of grid auto grid = RegularGrid{grid::LinearSpacing{xrange_outer[0], xrange_outer[1], nx}, grid::LinearSpacing{yrange_outer[0], yrange_outer[1], ny}, make_var_ratio_projection(1.13)}; + ///< check over regular grid points stretched using new atlas object and check using look-up table for (idx_t j = 0; j < grid.ny(); ++j) { for (idx_t i = 0; i < grid.nx(); ++i) { @@ -330,7 +507,7 @@ CASE("var_ratio_rot = 1.13") { config.set("north_pole", north_pole); ///< definition of grid I have to rotate this - auto proj_st = make_var_ratio_projection_rot(1.13, {-176., 40.}); + auto proj_st = make_var_ratio_projection_rot(1.13, north_pole); auto grid = RegularGrid{grid::LinearSpacing{xrange_outer[0], xrange_outer[1], nx}, grid::LinearSpacing{yrange_outer[0], yrange_outer[1], ny}, proj_st}; @@ -422,16 +599,12 @@ CASE("var_ratio_rot_inv = 1.13") { auto grid_nr = RegularGrid{grid::LinearSpacing{xrange_outer[0], xrange_outer[1], nx}, grid::LinearSpacing{yrange_outer[0], yrange_outer[1], ny}, proj_st_nr}; - std::vector north_pole = {-176., 40.}; - config.set("north_pole", north_pole); ///< definition of grid I have to rotate this auto proj_st = make_var_ratio_projection_rot(1.13, {-176., 40.}); auto grid = RegularGrid{grid::LinearSpacing{xrange_outer[0], xrange_outer[1], nx}, grid::LinearSpacing{yrange_outer[0], yrange_outer[1], ny}, proj_st}; - Rotation rotation(config); - ///< Test the inverse for (idx_t j = 0; j < grid.ny(); ++j) { for (idx_t i = 0; i < grid.nx(); ++i) {