diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 3001357f7..e79acd6fb 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -721,6 +721,7 @@ util/Rotation.h util/Registry.h util/SphericalPolygon.cc util/SphericalPolygon.h +util/SquareMatrix.h util/UnitSphere.h util/vector.h util/VectorOfAbstract.h diff --git a/src/atlas/grid/CubedSphereGrid.h b/src/atlas/grid/CubedSphereGrid.h index 659abf64b..a0fddc603 100644 --- a/src/atlas/grid/CubedSphereGrid.h +++ b/src/atlas/grid/CubedSphereGrid.h @@ -295,9 +295,19 @@ class CubedSphereGrid : public Grid { // Return the size of the cubed sphere grid, where N is the number of grid boxes along the edge of a tile inline int N() const { return grid_->N(); } - // Return the number of tiles + /// @brief return tiles object. inline atlas::grid::CubedSphereTiles tiles() const { return grid_->tiles(); } + /// @brief return cubed sphere projection object. + inline const projection::detail::CubedSphereProjectionBase& cubedSphereProjection() const { + + const auto projPtr = + dynamic_cast + (projection().get()); + + return *projPtr; + }; + temporary::IterateTIJ tij() const { return temporary::IterateTIJ(*grid_); } const std::string& stagger() const { return grid_->stagger(); } diff --git a/src/atlas/grid/Tiles.cc b/src/atlas/grid/Tiles.cc index bc8665a14..6e06a8999 100644 --- a/src/atlas/grid/Tiles.cc +++ b/src/atlas/grid/Tiles.cc @@ -72,5 +72,13 @@ std::ostream& operator<<(std::ostream& os, const CubedSphereTiles& t) { return os; } +const PointXY& CubedSphereTiles::tileCentre(size_t t) const { + return get()->tileCentre(t); +} + +const JacobianXY& CubedSphereTiles::tileJacobian(size_t t) const { + return get()->tileJacobian(t); +} + } // namespace grid } // namespace atlas diff --git a/src/atlas/grid/Tiles.h b/src/atlas/grid/Tiles.h index 723e932d1..98d049218 100644 --- a/src/atlas/grid/Tiles.h +++ b/src/atlas/grid/Tiles.h @@ -28,6 +28,7 @@ class Hash; //--------------------------------------------------------------------------------------------------------------------- namespace atlas { +class JacobianXY; class PointXY; class PointLonLat; @@ -85,6 +86,13 @@ class CubedSphereTiles : DOXYGEN_HIDE(public util::ObjectHandle register_builder(FV3CubedSphereTiles::static_type()); } diff --git a/src/atlas/grid/detail/tiles/FV3Tiles.h b/src/atlas/grid/detail/tiles/FV3Tiles.h index 9a248f54e..ebb51ffe5 100644 --- a/src/atlas/grid/detail/tiles/FV3Tiles.h +++ b/src/atlas/grid/detail/tiles/FV3Tiles.h @@ -46,6 +46,11 @@ class FV3CubedSphereTiles : public CubedSphereTiles { virtual void print(std::ostream&) const override; + virtual const PointXY& tileCentre(size_t t) const override; + + virtual const JacobianXY& tileJacobian(size_t t) const override; + + private: }; diff --git a/src/atlas/grid/detail/tiles/LFRicTiles.cc b/src/atlas/grid/detail/tiles/LFRicTiles.cc index 8705eb244..e9545f48b 100644 --- a/src/atlas/grid/detail/tiles/LFRicTiles.cc +++ b/src/atlas/grid/detail/tiles/LFRicTiles.cc @@ -26,27 +26,24 @@ namespace detail { namespace { -static constexpr bool debug = false; // constexpr so compiler can optimize `if ( debug ) { ... }` out +constexpr bool debug = false; // constexpr so compiler can optimize `if ( debug ) { ... }` out +constexpr double epsilon = 1.e-12; -using atlas::projection::detail::ProjectionUtilities; +using projection::detail::ProjectionUtilities; -static bool is_tiny(const double& x) { - constexpr double epsilon = 1.e-12; +bool is_tiny(const double& x) { return (std::abs(x) < epsilon); } -static bool is_same(const double& x, const double& y, const double& tol = 1.0) { - constexpr double epsilon = 1.e-12; +bool is_same(const double& x, const double& y, const double& tol = 1.0) { return (std::abs(x - y) < epsilon * tol); } -static bool is_less(const double& lhs, const double& rhs) { - constexpr double epsilon = 1.e-12; +bool is_less(const double& lhs, const double& rhs) { return lhs < rhs - epsilon; } -static bool is_geq(const double& lhs, const double& rhs) { - constexpr double epsilon = 1.e-12; +bool is_geq(const double& lhs, const double& rhs) { return lhs >= rhs - epsilon; } @@ -56,38 +53,22 @@ void sphericalToCartesian(const double lonlat[], double xyz[]) { ProjectionUtilities::sphericalToCartesian(lonlat, xyz, crd_sys, radius); } -atlas::PointXY rotatePlus90AboutPt(const atlas::PointXY& xy, const atlas::PointXY& origin) { - return atlas::PointXY{-xy.y() + origin.x() + origin.y(), xy.x() - origin.x() + origin.y()}; +PointXY rotatePlus90AboutPt(const PointXY& xy, const PointXY& origin) { + return PointXY{-xy.y() + origin.x() + origin.y(), xy.x() - origin.x() + origin.y()}; } -atlas::PointXY rotateMinus90AboutPt(const atlas::PointXY& xy, const atlas::PointXY& origin) { - return atlas::PointXY{xy.y() + origin.x() - origin.y(), -xy.x() + origin.x() + origin.y()}; +PointXY rotateMinus90AboutPt(const PointXY& xy, const PointXY& origin) { + return PointXY{xy.y() + origin.x() - origin.y(), -xy.x() + origin.x() + origin.y()}; } -atlas::PointXY rotatePlus180AboutPt(const atlas::PointXY& xy, const atlas::PointXY& origin) { - return atlas::PointXY{2.0 * origin.x() - xy.x(), 2.0 * origin.y() - xy.y()}; +PointXY rotatePlus180AboutPt(const PointXY& xy, const PointXY& origin) { + return PointXY{2.0 * origin.x() - xy.x(), 2.0 * origin.y() - xy.y()}; } } // anonymous namespace // constructor -LFRicCubedSphereTiles::LFRicCubedSphereTiles(const eckit::Parametrisation&) { - botLeftTile_ = - std::array{atlas::PointXY{0., -45.}, atlas::PointXY{90, -45}, atlas::PointXY{180., -45.}, - atlas::PointXY{270, -45}, atlas::PointXY{0., 45.}, atlas::PointXY{0, -135.}}; - - botRightTile_ = - std::array{atlas::PointXY{90., -45.}, atlas::PointXY{180., -45}, atlas::PointXY{270., -45.}, - atlas::PointXY{360., -45}, atlas::PointXY{90., 45.}, atlas::PointXY{90., -135.}}; - - topLeftTile_ = - std::array{atlas::PointXY{0., 45.}, atlas::PointXY{90, 45}, atlas::PointXY{180., 45.}, - atlas::PointXY{270, 45}, atlas::PointXY{0., 135.}, atlas::PointXY{0, -45.}}; - - topRightTile_ = - std::array{atlas::PointXY{90., 45.}, atlas::PointXY{180., 45}, atlas::PointXY{270., 45.}, - atlas::PointXY{360., 45}, atlas::PointXY{90., 135.}, atlas::PointXY{90., -45.}}; -} +LFRicCubedSphereTiles::LFRicCubedSphereTiles(const eckit::Parametrisation&) {} std::array, 2> LFRicCubedSphereTiles::xy2abOffsets() const { return {{{0., 1., 2., 3., 0., 0.}, {1., 1., 1., 1., 2., 0.}}}; @@ -97,78 +78,78 @@ std::array, 2> LFRicCubedSphereTiles::ab2xyOffsets() const return {{{0., 90., 180., 270., 0., 0.}, {-45., -45., -45., -45., 45., -135.}}}; } -static void tile0Rotate(double xyz[]) { +void tile0Rotate(double xyz[]) { // Face 0, no rotation. } -static void tile1Rotate(double xyz[]) { +void tile1Rotate(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[YY]; xyz[YY] = xyz_in[XX]; } -static void tile2Rotate(double xyz[]) { +void tile2Rotate(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[XX]; xyz[YY] = -xyz_in[YY]; } -static void tile3Rotate(double xyz[]) { +void tile3Rotate(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = xyz_in[YY]; xyz[YY] = -xyz_in[XX]; } -static void tile4Rotate(double xyz[]) { +void tile4Rotate(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = xyz_in[ZZ]; xyz[ZZ] = -xyz_in[XX]; } -static void tile5Rotate(double xyz[]) { +void tile5Rotate(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[ZZ]; xyz[ZZ] = xyz_in[XX]; } -static void tile0RotateInverse(double xyz[]) { +void tile0RotateInverse(double xyz[]) { // Face 0, no rotation. } -static void tile1RotateInverse(double xyz[]) { +void tile1RotateInverse(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = xyz_in[YY]; xyz[YY] = -xyz_in[XX]; } -static void tile2RotateInverse(double xyz[]) { +void tile2RotateInverse(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[XX]; xyz[YY] = -xyz_in[YY]; } -static void tile3RotateInverse(double xyz[]) { +void tile3RotateInverse(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[YY]; xyz[YY] = xyz_in[XX]; } -static void tile4RotateInverse(double xyz[]) { +void tile4RotateInverse(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = -xyz_in[ZZ]; xyz[ZZ] = xyz_in[XX]; } -static void tile5RotateInverse(double xyz[]) { +void tile5RotateInverse(double xyz[]) { double xyz_in[3]; std::copy(xyz, xyz + 3, xyz_in); xyz[XX] = xyz_in[ZZ]; @@ -380,8 +361,8 @@ void LFRicCubedSphereTiles::enforceXYdomain(double xy[]) const { // input is the xy value as PointXY that is a continous "cross " extension in terms of xy space from the tile in question // the output is an xy value that lives on the standard "|---" shape -atlas::PointXY LFRicCubedSphereTiles::tileCubePeriodicity(const atlas::PointXY& xyExtended, - const atlas::idx_t t) const { +PointXY LFRicCubedSphereTiles::tileCubePeriodicity(const PointXY& xyExtended, + const idx_t t) const { // xy space is a function of tiles--Tile 0) // xy space for Tile 1 // // // y ^ // y ^ @@ -484,24 +465,24 @@ atlas::PointXY LFRicCubedSphereTiles::tileCubePeriodicity(const atlas::PointXY& << " is not in the cross extension of tile " << t); } - atlas::PointXY withinRange = xyExtended; + PointXY withinRange = xyExtended; enforceWrapAround(t, withinRange); - atlas::PointXY finalXY = withinRange; - atlas::PointXY tempXY = withinRange; + PointXY finalXY = withinRange; + PointXY tempXY = withinRange; switch (t) { case 0: - finalXY = (withinRange.y() > 135.0) ? rotatePlus180AboutPt(withinRange, atlas::PointXY{135.0, 90.0}) + finalXY = (withinRange.y() > 135.0) ? rotatePlus180AboutPt(withinRange, PointXY{135.0, 90.0}) : withinRange; break; case 1: if ((withinRange.x() >= 90.0) && (withinRange.x() <= 180.0)) { if (withinRange.y() >= 45.0) { - tempXY = rotatePlus90AboutPt(withinRange, atlas::PointXY{90.0, 45.0}); + tempXY = rotatePlus90AboutPt(withinRange, PointXY{90.0, 45.0}); if (withinRange.y() > 135.0) { - finalXY = rotatePlus90AboutPt(tempXY, atlas::PointXY{0.0, 45.0}); + finalXY = rotatePlus90AboutPt(tempXY, PointXY{0.0, 45.0}); finalXY.x() += 360.0; } else { @@ -509,22 +490,22 @@ atlas::PointXY LFRicCubedSphereTiles::tileCubePeriodicity(const atlas::PointXY& } } else if (withinRange.y() < -45.0) { - finalXY = rotateMinus90AboutPt(withinRange, atlas::PointXY{90.0, -45.0}); + finalXY = rotateMinus90AboutPt(withinRange, PointXY{90.0, -45.0}); } } break; case 2: if ((withinRange.x() >= 180.0) && (withinRange.x() <= 270.0)) { if (withinRange.y() > 135.0) { - finalXY = rotatePlus180AboutPt(tempXY, atlas::PointXY{135.0, 90.0}); + finalXY = rotatePlus180AboutPt(tempXY, PointXY{135.0, 90.0}); } else if (withinRange.y() >= 45.0) { - tempXY = rotatePlus90AboutPt(withinRange, atlas::PointXY{180.0, 45.0}); - finalXY = rotatePlus90AboutPt(tempXY, atlas::PointXY{90.0, 45.0}); + tempXY = rotatePlus90AboutPt(withinRange, PointXY{180.0, 45.0}); + finalXY = rotatePlus90AboutPt(tempXY, PointXY{90.0, 45.0}); } else if (withinRange.y() < -45.0) { - tempXY = rotateMinus90AboutPt(withinRange, atlas::PointXY{180.0, -45.0}); - finalXY = rotateMinus90AboutPt(tempXY, atlas::PointXY{90.0, -45.0}); + tempXY = rotateMinus90AboutPt(withinRange, PointXY{180.0, -45.0}); + finalXY = rotateMinus90AboutPt(tempXY, PointXY{90.0, -45.0}); } } break; @@ -532,47 +513,47 @@ atlas::PointXY LFRicCubedSphereTiles::tileCubePeriodicity(const atlas::PointXY& case 3: if (((withinRange.x() >= 270.0) && (withinRange.x() <= 360.0)) || withinRange.x() == 0.0) { if (withinRange.y() > 135.0) { - finalXY = rotatePlus180AboutPt(tempXY, atlas::PointXY{225.0, 90.0}); + finalXY = rotatePlus180AboutPt(tempXY, PointXY{225.0, 90.0}); } else if (withinRange.y() >= 45.0) { - finalXY = rotateMinus90AboutPt(tempXY, atlas::PointXY{360.0, 45.0}); + finalXY = rotateMinus90AboutPt(tempXY, PointXY{360.0, 45.0}); finalXY.x() -= 360.0; } else if (withinRange.y() < -45.0) { - finalXY = rotatePlus90AboutPt(tempXY, atlas::PointXY{360.0, -45.0}); + finalXY = rotatePlus90AboutPt(tempXY, PointXY{360.0, -45.0}); finalXY.x() -= 360.0; } } break; case 4: if (withinRange.y() > 135.0) { - finalXY = rotatePlus180AboutPt(tempXY, atlas::PointXY{135.0, 90.0}); + finalXY = rotatePlus180AboutPt(tempXY, PointXY{135.0, 90.0}); } else if ((withinRange.y() >= 45.0) && (withinRange.y() <= 135.0)) { if ((withinRange.x() > 90.0) && (withinRange.x() <= 180.0)) { - finalXY = rotateMinus90AboutPt(withinRange, atlas::PointXY{90.0, 45.0}); + finalXY = rotateMinus90AboutPt(withinRange, PointXY{90.0, 45.0}); } if ((withinRange.x() > 180.0) && (withinRange.x() <= 270.0)) { - finalXY = rotatePlus180AboutPt(withinRange, atlas::PointXY{135.0, 0.0}); + finalXY = rotatePlus180AboutPt(withinRange, PointXY{135.0, 0.0}); } if ((withinRange.x() > 270.0) && (withinRange.x() <= 360.0)) { - finalXY = rotatePlus90AboutPt(withinRange, atlas::PointXY{360.0, 45.0}); + finalXY = rotatePlus90AboutPt(withinRange, PointXY{360.0, 45.0}); } } break; case 5: if (withinRange.y() > 135.0) { - finalXY = rotatePlus180AboutPt(tempXY, atlas::PointXY{135.0, 90.0}); + finalXY = rotatePlus180AboutPt(tempXY, PointXY{135.0, 90.0}); } else if ((withinRange.y() <= -45.0) && (withinRange.y() >= -135.0)) { if ((withinRange.x() > 90.0) && (withinRange.x() <= 180.0)) { - finalXY = rotatePlus90AboutPt(withinRange, atlas::PointXY{90.0, -45.0}); + finalXY = rotatePlus90AboutPt(withinRange, PointXY{90.0, -45.0}); } if ((withinRange.x() > 180.0) && (withinRange.x() <= 270.0)) { - finalXY = rotatePlus180AboutPt(withinRange, atlas::PointXY{135.0, 0.0}); + finalXY = rotatePlus180AboutPt(withinRange, PointXY{135.0, 0.0}); } if ((withinRange.x() > 270.0) && (withinRange.x() <= 360.0)) { - finalXY = rotateMinus90AboutPt(withinRange, atlas::PointXY{360.0, -45.0}); + finalXY = rotateMinus90AboutPt(withinRange, PointXY{360.0, -45.0}); } } break; @@ -628,39 +609,39 @@ void LFRicCubedSphereTiles::print(std::ostream& os) const { } -bool LFRicCubedSphereTiles::withinCross(const atlas::idx_t tiles, const atlas::PointXY& withinRange) const { +bool LFRicCubedSphereTiles::withinCross(const idx_t tiles, const PointXY& withinRange) const { std::size_t t = static_cast(tiles); - return !((withinRange.x() < botLeftTile_[t].x() && withinRange.y() < botLeftTile_[t].y()) || - (withinRange.x() > botRightTile_[t].x() && withinRange.y() < botRightTile_[t].y()) || - (withinRange.x() > topRightTile_[t].x() && withinRange.y() > topRightTile_[t].y()) || - (withinRange.x() < topLeftTile_[t].x() && withinRange.y() > topLeftTile_[t].y())); + return !((withinRange.x() < botLeftTile(t).x() && withinRange.y() < botLeftTile(t).y()) || + (withinRange.x() > botRightTile(t).x() && withinRange.y() < botRightTile(t).y()) || + (withinRange.x() > topRightTile(t).x() && withinRange.y() > topRightTile(t).y()) || + (withinRange.x() < topLeftTile(t).x() && withinRange.y() > topLeftTile(t).y())); } -void LFRicCubedSphereTiles::enforceWrapAround(const atlas::idx_t t, atlas::PointXY& withinRange) const { +void LFRicCubedSphereTiles::enforceWrapAround(const idx_t t, PointXY& withinRange) const { if (withinRange.x() < 0.0) { - atlas::PointXY temp = withinRange; + PointXY temp = withinRange; temp.x() += 360; if (withinCross(t, temp)) { withinRange = temp; } } if (withinRange.x() > 360.0) { - atlas::PointXY temp = withinRange; + PointXY temp = withinRange; temp.x() -= 360; if (withinCross(t, temp)) { withinRange = temp; } } if (withinRange.y() <= -135.0) { - atlas::PointXY temp = withinRange; + PointXY temp = withinRange; temp.y() += 360; if (withinCross(t, temp)) { withinRange = temp; } } if (withinRange.y() > 225.0) { - atlas::PointXY temp = withinRange; + PointXY temp = withinRange; temp.y() -= 360; if (withinCross(t, temp)) { withinRange = temp; @@ -670,6 +651,52 @@ void LFRicCubedSphereTiles::enforceWrapAround(const atlas::idx_t t, atlas::Point return; } +const PointXY& LFRicCubedSphereTiles::tileCentre(size_t t) const { + return tileCentres_[t]; +} + +const JacobianXY& LFRicCubedSphereTiles::tileJacobian(size_t t) const { + return tileJacobians_[t]; +} + +PointXY LFRicCubedSphereTiles::botLeftTile(size_t t) { + return tileCentres_[t] + PointXY{-45., -45.}; +} + +PointXY LFRicCubedSphereTiles::botRightTile(size_t t) { + return tileCentres_[t] + PointXY{45., -45.}; +} + +PointXY LFRicCubedSphereTiles::topLeftTile(size_t t) { + return tileCentres_[t] + PointXY{-45., 45.}; +} + +PointXY LFRicCubedSphereTiles::topRightTile(size_t t) { + return tileCentres_[t] + PointXY{45., 45.}; +} + + +// Centre of each tile in xy space. +const std::array LFRicCubedSphereTiles::tileCentres_ { + PointXY{45., 0.}, + PointXY{135., 0.}, + PointXY{225., 0.}, + PointXY{315., 0.}, + PointXY{45., 90.}, + PointXY{45., -90.} +}; + +// Jacobian of xy space with respect to curvilinear coordinates for each tile. +const std::array LFRicCubedSphereTiles::tileJacobians_{ + JacobianXY{{1., 0.}, {0., 1.}}, + JacobianXY{{1., 0.}, {0., 1.}}, + JacobianXY{{0., -1.}, {1., 0.}}, + JacobianXY{{0., -1.}, {1., 0.}}, + JacobianXY{{1., 0.}, {0., 1.}}, + JacobianXY{{0., 1.}, {-1., 0.}} +}; + + namespace { static CubedSphereTilesBuilder register_builder(LFRicCubedSphereTiles::static_type()); } diff --git a/src/atlas/grid/detail/tiles/LFRicTiles.h b/src/atlas/grid/detail/tiles/LFRicTiles.h index c1c3aef81..1846a856b 100644 --- a/src/atlas/grid/detail/tiles/LFRicTiles.h +++ b/src/atlas/grid/detail/tiles/LFRicTiles.h @@ -43,20 +43,32 @@ class LFRicCubedSphereTiles : public CubedSphereTiles { virtual void enforceXYdomain(double xy[]) const override; - virtual atlas::PointXY tileCubePeriodicity(const atlas::PointXY& xyExtended, - const atlas::idx_t tile) const override; + virtual PointXY tileCubePeriodicity(const PointXY& xyExtended, + const idx_t tile) const override; virtual void print(std::ostream&) const override; + virtual const PointXY& tileCentre(size_t t) const override; + + virtual const JacobianXY& tileJacobian(size_t t) const override; + private: - std::array botLeftTile_; - std::array botRightTile_; - std::array topLeftTile_; - std::array topRightTile_; - bool withinCross(const atlas::idx_t t, const atlas::PointXY& withinRange) const; + static PointXY botLeftTile(size_t t); + static PointXY botRightTile(size_t t); + static PointXY topLeftTile(size_t t); + static PointXY topRightTile(size_t t); + + bool withinCross(const idx_t t, const PointXY& withinRange) const; + + void enforceWrapAround(const idx_t t, PointXY& withinRange) const; + + + // Centre of each tile in xy-space. + static const std::array tileCentres_; + // Jacobian of xy with respect to tile curvilinear coordinates. + static const std::array tileJacobians_; - void enforceWrapAround(const atlas::idx_t t, atlas::PointXY& withinRange) const; }; } // namespace detail diff --git a/src/atlas/grid/detail/tiles/Tiles.cc b/src/atlas/grid/detail/tiles/Tiles.cc index 3d824452a..ea70de579 100644 --- a/src/atlas/grid/detail/tiles/Tiles.cc +++ b/src/atlas/grid/detail/tiles/Tiles.cc @@ -22,7 +22,7 @@ namespace grid { namespace detail { const CubedSphereTiles* CubedSphereTiles::create() { - // default: FV3 version (for now) + // default: LFRic version (for now) util::Config params; params.set("type", "cubedsphere_lfric"); return CubedSphereTiles::create(params); diff --git a/src/atlas/grid/detail/tiles/Tiles.h b/src/atlas/grid/detail/tiles/Tiles.h index 74c9f304f..bcf03af63 100644 --- a/src/atlas/grid/detail/tiles/Tiles.h +++ b/src/atlas/grid/detail/tiles/Tiles.h @@ -16,6 +16,7 @@ #include "atlas/library/config.h" #include "atlas/util/Config.h" #include "atlas/util/Object.h" +#include "atlas/util/SquareMatrix.h" #include "atlas/util/Point.h" namespace eckit { @@ -53,6 +54,10 @@ class CubedSphereTiles : public util::Object { virtual void enforceXYdomain(double xy[]) const = 0; + virtual const PointXY& tileCentre(size_t t) const = 0; + + virtual const JacobianXY& tileJacobian(size_t t) const = 0; + idx_t size() const { return 6; } virtual atlas::PointXY tileCubePeriodicity(const atlas::PointXY& xyExtended, const atlas::idx_t tile) const = 0; diff --git a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc index b91b3466a..7d5782417 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc @@ -139,8 +139,6 @@ class IJ { IJ(idx_t i, idx_t j) : i_(i), j_(j) {} idx_t i() const {return i_;} idx_t j() const {return j_;} - idx_t& i() {return i_;} - idx_t& j() {return j_;} IJ operator+(const IJ& ij) const {return IJ{i() + ij.i(), j() + ij.j()};} IJ operator-(const IJ& ij) const {return IJ{i() - ij.i(), j() - ij.j()};} private: diff --git a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc index 453f01101..2cc7325da 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc @@ -274,7 +274,7 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons const idx_t nCellsTotal = nCellsArray - 6 * 4 * nHalo * nHalo; // Projection and jacobian. - const auto* const csProjection = castProjection(csGrid.projection().get()); + const auto& csProjection = csGrid.cubedSphereProjection(); const auto jacobian = NeighbourJacobian(csGrid); // Get partition information. @@ -487,7 +487,7 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons // This will only determine if tGlobal does not match t. // This is cheaper than determining the correct tGlobal. - idx_t tGlobal = csProjection->getCubedSphereTiles().indexFromXY(xy.data()); + idx_t tGlobal = csProjection.getCubedSphereTiles().indexFromXY(xy.data()); if (tGlobal == t) { // Node is an owner. @@ -773,7 +773,7 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons nodesXy(nodeLocalIdx, YY) = xyLocal.y(); // Set lon-lat. - const PointLonLat lonLat = csProjection->lonlat(xyGlobal); + const PointLonLat lonLat = csProjection.lonlat(xyGlobal); nodesLonLat(nodeLocalIdx, LON) = lonLat.lon(); nodesLonLat(nodeLocalIdx, LAT) = lonLat.lat(); @@ -884,7 +884,7 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons cellsXy(cellLocalIdx, YY) = xyLocal.y(); // Set lon-lat. - const PointLonLat lonLat = csProjection->lonlat(xyGlobal); + const PointLonLat lonLat = csProjection.lonlat(xyGlobal); cellsLonLat(cellLocalIdx, LON) = lonLat.lon(); cellsLonLat(cellLocalIdx, LAT) = lonLat.lat(); diff --git a/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc b/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc index c2bd5db9c..2ed0290d2 100644 --- a/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc @@ -98,9 +98,9 @@ void NodalCubedSphereMeshGenerator::generate(const Grid& grid, const grid::Distr } // Get tiles - auto csprojection = castProjection(csgrid.projection().get()); + const auto& csprojection = csgrid.cubedSphereProjection(); // grid must use FV3Tiles class. - if (csprojection->getCubedSphereTiles().type() != "cubedsphere_fv3") { + if (csprojection.getCubedSphereTiles().type() != "cubedsphere_fv3") { throw_Exception("NodalCubedSphereMeshGenerator only works with FV3 tiles", Here()); } diff --git a/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.cc b/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.cc index 733e37018..921e8600d 100644 --- a/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.cc +++ b/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.cc @@ -10,69 +10,12 @@ #include "atlas/grid/Iterator.h" #include "atlas/projection/detail/CubedSphereProjectionBase.h" + namespace atlas { namespace meshgenerator { namespace detail { namespace cubedsphere { - -// ----------------------------------------------------------------------------- -// Projection cast -// ----------------------------------------------------------------------------- - -const CubedSphereProjectionBase* castProjection(const ProjectionImpl* projectionPtr) { - const auto* const csProjectionPtr = dynamic_cast(projectionPtr); - - if (!csProjectionPtr) { - throw_Exception("Cannot cast ProjectionImpl* to CubedSphereProjectionBase*.", Here()); - } - return csProjectionPtr; -} - -// ----------------------------------------------------------------------------- -// Jacobian2 class -// ----------------------------------------------------------------------------- - -Jacobian2::Jacobian2(double df0_by_dx0, double df0_by_dx1, double df1_by_dx0, double df1_by_dx1): - df0_by_dx0_(df0_by_dx0), df0_by_dx1_(df0_by_dx1), df1_by_dx0_(df1_by_dx0), df1_by_dx1_(df1_by_dx1) {} - -Jacobian2::Jacobian2(const Point2& f00, const Point2& f10, const Point2& f01, double dx0, double dx1): - Jacobian2((f10[0] - f00[0]) / dx0, (f01[0] - f00[0]) / dx1, (f10[1] - f00[1]) / dx0, (f01[1] - f00[1]) / dx1) {} - -Jacobian2::Jacobian2(const Point2& f00, const Point2& f10, const Point2& f01): Jacobian2(f00, f10, f01, 1., 1.) {} - -double Jacobian2::det() const { - return df0_by_dx0_ * df1_by_dx1_ - df0_by_dx1_ * df1_by_dx0_; -} - -Jacobian2 Jacobian2::operator*(double a) const { - return Jacobian2(df0_by_dx0_ * a, df0_by_dx1_ * a, df1_by_dx0_ * a, df1_by_dx1_ * a); -} - -Point2 Jacobian2::operator*(const Point2& dx) const { - return Point2(dx[0] * df0_by_dx0_ + df0_by_dx1_ * dx[1], dx[0] * df1_by_dx0_ + df1_by_dx1_ * dx[1]); -} - -Jacobian2 Jacobian2::operator*(const Jacobian2& Jb) const { - return Jacobian2(df0_by_dx0_ * Jb.df0_by_dx0_ + df0_by_dx1_ * Jb.df1_by_dx0_, - df0_by_dx0_ * Jb.df0_by_dx1_ + df0_by_dx1_ * Jb.df1_by_dx1_, - df1_by_dx0_ * Jb.df0_by_dx0_ + df1_by_dx1_ * Jb.df1_by_dx0_, - df1_by_dx0_ * Jb.df0_by_dx1_ + df1_by_dx1_ * Jb.df1_by_dx1_); -} - -Jacobian2 Jacobian2::inverse() const { - return Jacobian2(df1_by_dx1_, -df0_by_dx1_, -df1_by_dx0_, df0_by_dx0_) * (1. / det()); -} - -Jacobian2 Jacobian2::sign() const { - const double smallNumber = det() * std::numeric_limits::epsilon(); - const auto signValue = [&](double number) -> double { - return std::abs(number) < smallNumber ? 0. : number < 0. ? -1. : 1.; - }; - return Jacobian2(signValue(df0_by_dx0_), signValue(df0_by_dx1_), signValue(df1_by_dx0_), signValue(df1_by_dx1_)); -} - - // ----------------------------------------------------------------------------- // NeighbourJacobian class // ----------------------------------------------------------------------------- @@ -89,7 +32,7 @@ NeighbourJacobian::NeighbourJacobian(const CubedSphereGrid& csGrid) { } // Get projection. - csProjection_ = castProjection(csGrid.projection().get()); + csProjection_ = &csGrid.cubedSphereProjection(); // Get tiles. const auto& csTiles = csProjection_->getCubedSphereTiles(); @@ -101,78 +44,42 @@ NeighbourJacobian::NeighbourJacobian(const CubedSphereGrid& csGrid) { // Get cell width. const double cellWidth = 90. / N_; - // Get xy of points (i = 0, j = 0), (i = 1, j = 0) and (i = 0, j = 1) on tiles. - std::array xy00; - std::array xy10; - std::array xy01; - - // Loop over grid points. - auto tijIt = csGrid.tij().begin(); - for (const PointXY& xy : csGrid.xy()) { - const auto t = static_cast((*tijIt).t()); - const idx_t i = (*tijIt).i(); - const idx_t j = (*tijIt).j(); - - - if (i == 0 && j == 0) { - xy00[t] = xy; - } - else if (i == 1 && j == 0) { - xy10[t] = xy; - } - else if (i == 0 && j == 1) { - xy01[t] = xy; - } - ++tijIt; - } for (size_t t = 0; t < 6; ++t) { // Calculate tile Jacobians. - dxy_by_dij_[t] = Jacobian2(xy00[t], xy10[t], xy01[t]); - - // Rescale by cell width (gains an extra couple of decimal places of precision). - dxy_by_dij_[t] = dxy_by_dij_[t].sign() * cellWidth; + dxy_by_dij_[t] = csTiles.tileJacobian(t) * cellWidth; // Get inverse. dij_by_dxy_[t] = dxy_by_dij_[t].inverse(); // Set xy00. Grid point needs moving to (i = 0, j = 0). - xy00_[t] = xy00[t] + dxy_by_dij_[t] * PointIJ(-0.5, -0.5); - - // Get other three corners so we can work out xy min/max. - const PointXY xyN0 = xy00_[t] + dxy_by_dij_[t] * PointIJ(N_, 0); - const PointXY xyNN = xy00_[t] + dxy_by_dij_[t] * PointIJ(N_, N_); - const PointXY xy0N = xy00_[t] + dxy_by_dij_[t] * PointIJ(0, N_); + xy00_[t] = csTiles.tileCentre(t) + dxy_by_dij_[t] * PointIJ(-0.5 * N_, -0.5 * N_); // Get xy min/max. - std::tie(xyMin_[t].x(), xyMax_[t].x()) = std::minmax({xy00_[t].x(), xyN0.x(), xyNN.x(), xy0N.x()}); - std::tie(xyMin_[t].y(), xyMax_[t].y()) = std::minmax({xy00_[t].y(), xyN0.y(), xyNN.y(), xy0N.y()}); - - // Round to nearest degree. - xyMin_[t].x() = std::round(xyMin_[t].x()); - xyMax_[t].x() = std::round(xyMax_[t].x()); - xyMin_[t].y() = std::round(xyMin_[t].y()); - xyMax_[t].y() = std::round(xyMax_[t].y()); + xyMin_[t].x() = csTiles.tileCentre(t).x() + 45.; + xyMax_[t].x() = csTiles.tileCentre(t).x() - 45.; + xyMin_[t].y() = csTiles.tileCentre(t).y() + 45.; + xyMax_[t].y() = csTiles.tileCentre(t).y() - 45.; // Neighbour assignment lambda. const auto neighbourAssignment = [&](TileEdge::k k) -> void { // Shift points in to neighbouring tiles. - PointIJ ijDisplacement; + PointIJ Dij; switch (k) { case TileEdge::LEFT: { - ijDisplacement = PointIJ(-2, 0); + Dij = PointIJ(-2, 0); break; } case TileEdge::BOTTOM: { - ijDisplacement = PointIJ(0, -2); + Dij = PointIJ(0, -2); break; } case TileEdge::RIGHT: { - ijDisplacement = PointIJ(N_, 0); + Dij = PointIJ(N_, 0); break; } case TileEdge::TOP: { - ijDisplacement = PointIJ(0, N_); + Dij = PointIJ(0, N_); break; } case TileEdge::UNDEFINED: { @@ -180,13 +87,15 @@ NeighbourJacobian::NeighbourJacobian(const CubedSphereGrid& csGrid) { } } - // Convert displacement from ij to xy. - const PointXY xyDisplacement = dxy_by_dij_[t] * ijDisplacement; + // half-index displacements. + const auto dij00 = PointIJ(0.5, 0.5); + const auto dij10 = PointIJ(1.5, 0.5); + const auto dij01 = PointIJ(0.5, 1.5); // Get neighbour xy points in xy space local to tile. - const PointXY xy00Local = xy00[t] + xyDisplacement; - const PointXY xy10Local = xy10[t] + xyDisplacement; - const PointXY xy01Local = xy01[t] + xyDisplacement; + const PointXY xy00Local = xy00_[t] + dxy_by_dij_[t] * (Dij + dij00); + const PointXY xy10Local = xy00_[t] + dxy_by_dij_[t] * (Dij + dij10); + const PointXY xy01Local = xy00_[t] + dxy_by_dij_[t] * (Dij + dij01); // Convert from local xy to global xy. const PointXY xy00Global = csTiles.tileCubePeriodicity(xy00Local, static_cast(t)); @@ -197,7 +106,7 @@ NeighbourJacobian::NeighbourJacobian(const CubedSphereGrid& csGrid) { neighbours_[t].t_[k] = csTiles.indexFromXY(xy00Global.data()); // Set Jacobian of global xy with respect to local ij. - auto dxyGlobal_by_dij = Jacobian2(xy00Global, xy10Global, xy01Global); + auto dxyGlobal_by_dij = JacobianXY(xy00Global, xy10Global, xy01Global); // Rescale by cell width (gains an extra couple of decimal places of precision). dxyGlobal_by_dij = dxyGlobal_by_dij.sign() * cellWidth; @@ -222,7 +131,7 @@ NeighbourJacobian::NeighbourJacobian(const CubedSphereGrid& csGrid) { PointXY NeighbourJacobian::xy(const PointIJ& ij, idx_t t) const { // Get jacobian. - const Jacobian2& jac = dxy_by_dij_[static_cast(t)]; + const JacobianXY& jac = dxy_by_dij_[static_cast(t)]; const PointXY& xy00 = xy00_[static_cast(t)]; // Return ij @@ -231,7 +140,7 @@ PointXY NeighbourJacobian::xy(const PointIJ& ij, idx_t t) const { PointIJ NeighbourJacobian::ij(const PointXY& xy, idx_t t) const { // Get jacobian. - const Jacobian2& jac = dij_by_dxy_[static_cast(t)]; + const JacobianXY& jac = dij_by_dxy_[static_cast(t)]; const PointXY& xy00 = xy00_[static_cast(t)]; // Return ij @@ -290,7 +199,7 @@ PointTXY NeighbourJacobian::xyLocalToGlobal(const PointXY& xyLocal, idx_t tLocal // Get reference points and jacobian. const PointXY& xy00Local_ = neighbours_[static_cast(tLocal)].xy00Local_[k]; const PointXY& xy00Global_ = neighbours_[static_cast(tLocal)].xy00Global_[k]; - const Jacobian2& jac = neighbours_[static_cast(tLocal)].dxyGlobal_by_dxyLocal_[k]; + const JacobianXY& jac = neighbours_[static_cast(tLocal)].dxyGlobal_by_dxyLocal_[k]; // Get t. tGlobal = neighbours_[static_cast(tLocal)].t_[k]; diff --git a/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.h b/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.h index bcaff0e8a..452656845 100644 --- a/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.h +++ b/src/atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.h @@ -9,6 +9,7 @@ #include "atlas/library/config.h" #include "atlas/runtime/Exception.h" +#include "atlas/util/SquareMatrix.h" #include "atlas/util/Point.h" namespace atlas { @@ -55,9 +56,6 @@ struct TileEdge { }; }; -/// Cast Projection to CubedSphereProjectionBase. -const CubedSphereProjectionBase* castProjection(const ProjectionImpl* projectionPtr); - /// Class to store (i, j) indices as a Point2 coordinate. class PointIJ : public Point2 { public: @@ -115,61 +113,6 @@ class PointTIJ : public std::pair { const PointIJ& ij() const { return second; } }; -/// \brief Jacobian class for 2 dimensional vector fields. -/// -/// \details Wrapper class for an Eigen matrix which stores the partial -/// of f(x). Objects can be constructed directly from the four partial -/// derivatives, or by supplying three Point2 objects with the -/// following relative positions: -/// -/// ^ -/// | *f(X0, X1 + dx1) -/// | -/// x1 | -/// | -/// | *f(X0, X1) *f(X0 + dx0, X1) -/// +----------------------------> -/// x0 -class Jacobian2 { -public: - /// Default constructor. - Jacobian2() = default; - - /// Partial derivative constructor. - Jacobian2(double df0_by_dx0, double df0_by_dx1, double df1_by_dx0, double df1_by_dx1); - - /// Discrete point constructor (explicit dx). - Jacobian2(const Point2& f00, const Point2& f10, const Point2& f01, double dx0, double dx1); - - /// Discrete point contructor (implicit dx). - Jacobian2(const Point2& f00, const Point2& f10, const Point2& f01); - - /// Determinant of Jacobian. - double det() const; - - /// Jacobian-scalar multiplication. - Jacobian2 operator*(double a) const; - - /// Jacobian-vector multiplication. - Point2 operator*(const Point2& dx) const; - - /// Jacobian-Jacobian multiplication. - Jacobian2 operator*(const Jacobian2& Jb) const; - - /// Inverse Jacobian (partial derivatives of x(f)). - Jacobian2 inverse() const; - - /// Get signed elements of matrix (i.e, 0, +1, -1). - Jacobian2 sign() const; - -private: - // Data storage. - double df0_by_dx0_; - double df0_by_dx1_; - double df1_by_dx0_; - double df1_by_dx1_; -}; - /// \brief Class to convert between ij and xy on a tile and its four /// surrounding neighbours. /// @@ -217,10 +160,10 @@ class NeighbourJacobian { idx_t N_{}; // Jacobian of xy with respect to ij for each tile. - std::array dxy_by_dij_{}; + std::array dxy_by_dij_{}; // Jacobian of ij with respect to xy for each tile. - std::array dij_by_dxy_{}; + std::array dij_by_dxy_{}; // Lower-left xy position on each tile. std::array xy00_{}; @@ -237,7 +180,7 @@ class NeighbourJacobian { std::array t_{}; // Jacobian of remote xy with respect to local xy. - std::array dxyGlobal_by_dxyLocal_{}; + std::array dxyGlobal_by_dxyLocal_{}; // Lower left most local xy position on neighbour tiles. std::array xy00Local_; diff --git a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc index d8748feb3..49309f47e 100644 --- a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc +++ b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.cc @@ -21,12 +21,42 @@ #include "atlas/util/Config.h" #include "atlas/util/Constants.h" #include "atlas/util/CoordinateEnums.h" +#include "atlas/util/SquareMatrix.h" namespace { static constexpr bool debug = false; // constexpr so compiler can optimize `if ( debug ) { ... }` out static constexpr double deg2rad = atlas::util::Constants::degreesToRadians(); static constexpr double rad2deg = atlas::util::Constants::radiansToDegrees(); + +// Define small number relative to 360. +constexpr double epsilon = std::numeric_limits::epsilon() * 360.; + +// Define some "fuzzy" comparison operators. + +// a is approximately equal to b. +bool equal(double a, double b) { + return std::abs(a - b) <= epsilon; +} +// a is less than b. +bool lessThan(double a, double b) { + return a < b && !equal(a, b); +} +// a is greater than b. +bool greaterThan(double a, double b) { + return a > b && !equal(a, b); +} +// a is less than or approximately equal to b. +bool lessEqual(double a, double b) { + return a < b || equal(a, b); +} +// a is greater than or approximately equal to b. +bool greaterEqual(double a, double b) { + return a > b || equal(a, b); +} + + + } // namespace namespace atlas { @@ -39,6 +69,105 @@ namespace detail { CubedSphereEquiAnglProjection::CubedSphereEquiAnglProjection(const eckit::Parametrisation& params): CubedSphereProjectionBase(params) {} +// ------------------------------------------------------------------------------------------------- + +void CubedSphereEquiAnglProjection::xy2alphabeta(double crd[], idx_t t) const { + + // Get tile centre. + const auto& xyCentre = getCubedSphereTiles().tileCentre(static_cast(t)); + + // Check that xy coordinate is within valid "+" shaped halo region. + const auto inCross = [&](const double crd[]) -> bool { + return (greaterEqual(crd[XX], xyCentre[XX] - 45.) && lessEqual(crd[XX], xyCentre[XX] + 45.)) || + (greaterEqual(crd[YY], xyCentre[YY] - 45.) && lessEqual(crd[YY], xyCentre[YY] + 45.)); + }; + if (!inCross(crd)) { + auto sStream = std::stringstream(); + sStream << "xy coordinate (" << crd[0] << ", " << crd[1] << ") is not in range for tile " << t << "."; + ATLAS_THROW_EXCEPTION(sStream.str()); + } + + // Get alphaBeta Jacobian. + const auto alphabetaJacobian = + getCubedSphereTiles().tileJacobian(static_cast(t)).inverse(); + + // Set (alpha, beta) coord. + const Point2 alphabeta = alphabetaJacobian * (Point2(crd) - xyCentre); + crd[0] = alphabeta[0]; + crd[1] = alphabeta[1]; + + // Define correction. + const auto correction = [](const double crd[]) -> double { + return rad2deg * std::atan(std::tan(crd[0] * deg2rad) * + std::tan(crd[1] * deg2rad)); + }; + + // Correct halo (alpha, beta) coord. + if (lessThan(crd[0], -45.)) { + // Left. + crd[1] = -correction(crd); + } + else if (greaterThan(crd[0], 45.)) { + // Right. + crd[1] = correction(crd); + } + else if (lessThan(crd[1], -45.)) { + // Bottom. + crd[0] = -correction(crd); + } + else if (greaterThan(crd[1], 45.)) { + // Top. + crd[0] = correction(crd); + } +} + +// ------------------------------------------------------------------------------------------------- + +void CubedSphereEquiAnglProjection::alphabeta2xy(double crd[], idx_t t) const { + + + // Define correction. + const auto correction1 = [](const double crd[]) -> double { + return rad2deg * std::atan(std::tan(crd[1] * deg2rad) / + std::tan(crd[0] * deg2rad)); + }; + const auto correction2 = [](const double crd[]) -> double { + return rad2deg * std::atan(std::tan(crd[0] * deg2rad) / + std::tan(crd[1] * deg2rad)); + }; + + // Correct halo (alpha, beta) coord. + if (lessThan(crd[0], -45.) && greaterThan(crd[1], crd[0]) && lessEqual(crd[1], -crd[0])) { + // Left trapezium. + crd[1] = -correction1(crd); + } + else if (greaterThan(crd[0], 45.) && greaterEqual(crd[1], -crd[0]) && lessThan(crd[1], crd[0])) { + // Right trapezium. + crd[1] = correction1(crd); + } + else if (lessThan(crd[1], -45.) && greaterEqual(crd[0], crd[1]) && lessThan(crd[0], -crd[1])) { + // Bottom trapezium. + crd[0] = -correction2(crd); + } + else if (greaterThan(crd[1], 45.) && greaterThan(crd[0], -crd[1]) && lessEqual(crd[0], crd[1])) { + // Top trapezium. + crd[0] = correction2(crd); + } + + // Get tile centre. + const auto& xyCentre = getCubedSphereTiles().tileCentre(static_cast(t)); + + // Get xy Jacobian. + const auto xyJacobian = + getCubedSphereTiles().tileJacobian(static_cast(t)); + + // Set xy coord. + const Point2 xy = xyJacobian * Point2(crd) + xyCentre; + crd[XX] = xy[XX]; + crd[YY] = xy[YY]; +} + + // ------------------------------------------------------------------------------------------------- void CubedSphereEquiAnglProjection::lonlat2xy(double crd[]) const { diff --git a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h index 1d69870b2..1797abea8 100644 --- a/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h +++ b/src/atlas/projection/detail/CubedSphereEquiAnglProjection.h @@ -27,6 +27,12 @@ class CubedSphereEquiAnglProjection final : public CubedSphereProjectionBase { static std::string static_type() { return "cubedsphere_equiangular"; } std::string type() const override { return static_type(); } + /// @brief Convert (x, y) coordinate to (alpha, beta) on tile t. + void xy2alphabeta(double crd[], idx_t t) const override; + + /// @brief Convert (alpha, beta) coordinate to (x, y) on tile t. + void alphabeta2xy(double crd[], idx_t t) const override; + // projection and inverse projection void xy2lonlat(double crd[]) const override; void lonlat2xy(double crd[]) const override; diff --git a/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc b/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc index 925020d5b..0ca6e1e7c 100644 --- a/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc +++ b/src/atlas/projection/detail/CubedSphereEquiDistProjection.cc @@ -35,6 +35,18 @@ CubedSphereEquiDistProjection::CubedSphereEquiDistProjection(const eckit::Parame // ------------------------------------------------------------------------------------------------- +void CubedSphereEquiDistProjection::xy2alphabeta(double crd[], idx_t t) const { + throw_NotImplemented("xy2alphabeta not implemented for CubedSphereEquiDistProjection", Here()); +} + +// ------------------------------------------------------------------------------------------------- + +void CubedSphereEquiDistProjection::alphabeta2xy(double crd[], idx_t t) const { + throw_NotImplemented("alphabeta2xy not implemented for CubedSphereEquiDistProjection", Here()); +} + +// ------------------------------------------------------------------------------------------------- + void CubedSphereEquiDistProjection::lonlat2xy(double crd[]) const { if (debug) { Log::info() << "equidist lonlat2xy start : lonlat = " << crd[LON] << " " << crd[LAT] << std::endl; diff --git a/src/atlas/projection/detail/CubedSphereEquiDistProjection.h b/src/atlas/projection/detail/CubedSphereEquiDistProjection.h index 6fc59f1a7..136a91690 100644 --- a/src/atlas/projection/detail/CubedSphereEquiDistProjection.h +++ b/src/atlas/projection/detail/CubedSphereEquiDistProjection.h @@ -27,6 +27,12 @@ class CubedSphereEquiDistProjection final : public CubedSphereProjectionBase { static std::string static_type() { return "cubedsphere_equidistant"; } std::string type() const override { return static_type(); } + /// @brief Convert (x, y) coordinate to (alpha, beta) on tile t. + void xy2alphabeta(double crd[], idx_t t) const override; + + /// @brief Convert (alpha, beta) coordinate to (x, y) on tile t. + void alphabeta2xy(double crd[], idx_t t) const override; + // projection and inverse projection void xy2lonlat(double crd[]) const override; void lonlat2xy(double crd[]) const override; diff --git a/src/atlas/projection/detail/CubedSphereProjectionBase.cc b/src/atlas/projection/detail/CubedSphereProjectionBase.cc index ba05d72b6..a8295be4b 100644 --- a/src/atlas/projection/detail/CubedSphereProjectionBase.cc +++ b/src/atlas/projection/detail/CubedSphereProjectionBase.cc @@ -232,6 +232,22 @@ void CubedSphereProjectionBase::alphabetat2xy(const idx_t t, const double ab[], // ------------------------------------------------------------------------------------------------- +Point2 CubedSphereProjectionBase::xy2alphabeta(const Point2& xy, idx_t t) const { + auto alphabeta = Point2(xy); + xy2alphabeta(alphabeta.data(), t); + return alphabeta; +} + +// ------------------------------------------------------------------------------------------------- + +Point2 CubedSphereProjectionBase::alphabeta2xy(const Point2& alphabeta, idx_t t) const { + auto xy = Point2(alphabeta); + alphabeta2xy(xy.data(), t); + return xy; +} + +// ------------------------------------------------------------------------------------------------- + } // namespace detail } // namespace projection } // namespace atlas diff --git a/src/atlas/projection/detail/CubedSphereProjectionBase.h b/src/atlas/projection/detail/CubedSphereProjectionBase.h index 17c34a861..29b00cac9 100644 --- a/src/atlas/projection/detail/CubedSphereProjectionBase.h +++ b/src/atlas/projection/detail/CubedSphereProjectionBase.h @@ -30,6 +30,31 @@ class CubedSphereProjectionBase : public ProjectionImpl { atlas::grid::CubedSphereTiles getCubedSphereTiles() const { return tiles_; }; + /// @brief Convert (x, y) coordinate to (alpha, beta) on tile t. + /// + /// @details Converts the Atlas xy coordinates to the angular coordinates + /// described of tile t, described by by Ronchi et al. (1996, + /// Journal of Computational Physics, 124, 93). + /// Note that the xy coordinate must lie within the domain + /// (x <= xmin && x >= xmax) || (y <= ymin && y >= ymax) where + /// xmin, xmax, ymin and ymax are the boundaries of tile t. + ///@{ + Point2 xy2alphabeta(const Point2& xy, idx_t t) const; + virtual void xy2alphabeta(double crd[], idx_t t) const = 0; + ///@} + + /// @brief Convert (alpha, beta) coordinate to (x, y) on tile t. + /// + /// @details Performs the inverse of xy2alpha beta. Note that the result is + /// degenerate when abs(alpha) > 45 && abs(beta) > 45 && + /// abs(alpha) == abs(beta). In these circumstances, the method + /// will return the anticlockwise-most of the two possible + /// values. + ///@{ + Point2 alphabeta2xy(const Point2& alphabeta, idx_t t) const; + virtual void alphabeta2xy(double crd[], idx_t) const = 0; + ///@} + protected: // projection and inverse projection void xy2lonlat_post(double xyz[], const idx_t t, double crd[]) const; diff --git a/src/atlas/redistribution/detail/RedistributeGeneric.cc b/src/atlas/redistribution/detail/RedistributeGeneric.cc index 53254341a..39cf4d2ff 100644 --- a/src/atlas/redistribution/detail/RedistributeGeneric.cc +++ b/src/atlas/redistribution/detail/RedistributeGeneric.cc @@ -215,35 +215,42 @@ std::pair, std::vector> getUidIntersection(const std::ve // Iterate over a field, in the order of an index list, and apply a functor to // each element. -// Rank 1 overload. -template -void iterateField(const std::vector& idxList, array::ArrayView& fieldView, const Functor& f) { - for (const idx_t i : idxList) { - f(fieldView(i)); +// Recursive ForEach to visit all elements of field. +template +struct ForEach { + template + static void apply(const std::vector& idxList, array::ArrayView& fieldView, + const Functor& f, Idxs... idxs) { + // Iterate over dimension Dim of array. + for (idx_t idx = 0; idx < fieldView.shape(Dim); ++idx) { + ForEach::apply(idxList, fieldView, f, idxs..., idx); + } } -} +}; -// Rank 2 overload. -template -void iterateField(const std::vector& idxList, array::ArrayView& fieldView, const Functor f) { - for (const idx_t i : idxList) { - for (idx_t j = 0; j < fieldView.shape(1); ++j) { - f(fieldView(i, j)); +// Beginning of recursion when Dim == 0. +template +struct ForEach { + template + static void apply(const std::vector& idxList, array::ArrayView& fieldView, + const Functor& f, Idxs... idxs) { + // Iterate over dimension 0 of array in order defined by idxList. + for (idx_t idx : idxList) { + ForEach::apply(idxList, fieldView, f, idxs..., idx); } } -} +}; -// Rank 3 overload. -template -void iterateField(const std::vector& idxList, array::ArrayView& fieldView, const Functor f) { - for (const idx_t i : idxList) { - for (idx_t j = 0; j < fieldView.shape(1); ++j) { - for (idx_t k = 0; k < fieldView.shape(2); ++k) { - f(fieldView(i, j, k)); - } - } +// End of recursion when Dim == Rank. +template +struct ForEach { + template + static void apply(const std::vector& idxList, array::ArrayView& fieldView, + const Functor& f, Idxs... idxs) { + // Apply functor. + f(fieldView(idxs...)); } -} +}; } // namespace @@ -297,25 +304,23 @@ void RedistributeGeneric::execute(const FieldSet& sourceFieldSet, FieldSet& targ // Determine datatype. void RedistributeGeneric::do_execute(const Field& sourceField, Field& targetField) const { + + // Available datatypes defined in array/LocalView.cc switch (sourceField.datatype().kind()) { case array::DataType::KIND_REAL64: { - do_execute(sourceField, targetField); - break; + return do_execute(sourceField, targetField); } case array::DataType::KIND_REAL32: { - do_execute(sourceField, targetField); - break; + return do_execute(sourceField, targetField); } case array::DataType::KIND_INT64: { - do_execute(sourceField, targetField); - break; + return do_execute(sourceField, targetField); } case array::DataType::KIND_INT32: { - do_execute(sourceField, targetField); - break; + return do_execute(sourceField, targetField); } default: { - throw_NotImplemented("No implementation for data type " + sourceField.datatype().str(), Here()); + ATLAS_THROW_EXCEPTION("No implementation for data type " + sourceField.datatype().str()); } } } @@ -323,21 +328,20 @@ void RedistributeGeneric::do_execute(const Field& sourceField, Field& targetFiel // Determine rank. template void RedistributeGeneric::do_execute(const Field& sourceField, Field& targetField) const { + + // Available ranks defined in array/LocalView.cc switch (sourceField.rank()) { - case 1: { - do_execute(sourceField, targetField); - break; - } - case 2: { - do_execute(sourceField, targetField); - break; - } - case 3: { - do_execute(sourceField, targetField); - break; - } + case 1: {return do_execute(sourceField, targetField);} + case 2: {return do_execute(sourceField, targetField);} + case 3: {return do_execute(sourceField, targetField);} + case 4: {return do_execute(sourceField, targetField);} + case 5: {return do_execute(sourceField, targetField);} + case 6: {return do_execute(sourceField, targetField);} + case 7: {return do_execute(sourceField, targetField);} + case 8: {return do_execute(sourceField, targetField);} + case 9: {return do_execute(sourceField, targetField);} default: { - throw_NotImplemented("No implementation for rank " + std::to_string(sourceField.rank()), Here()); + ATLAS_THROW_EXCEPTION("No implementation for rank " + std::to_string(sourceField.rank())); } } } @@ -380,14 +384,17 @@ void RedistributeGeneric::do_execute(const Field& sourceField, Field& targetFiel auto recvBufferIt = recvBuffer.cbegin(); // Copy sourceField to sendBuffer. - iterateField(sourceLocalIdx_, sourceView, [&](const Value& elem) { *sendBufferIt++ = elem; }); + ForEach::apply(sourceLocalIdx_, sourceView, + [&](const Value& elem){*sendBufferIt++ = elem;}); // Perform MPI communication. mpi::comm().allToAllv(sendBuffer.data(), sendCounts.data(), sendDisps.data(), recvBuffer.data(), recvCounts.data(), recvDisps.data()); // Copy recvBuffer to targetField. - iterateField(targetLocalIdx_, targetView, [&](Value& elem) { elem = *recvBufferIt++; }); + ForEach::apply(targetLocalIdx_, targetView, + [&](Value& elem){elem = *recvBufferIt++;}); + } namespace { diff --git a/src/atlas/util/SquareMatrix.h b/src/atlas/util/SquareMatrix.h new file mode 100644 index 000000000..4bd91392c --- /dev/null +++ b/src/atlas/util/SquareMatrix.h @@ -0,0 +1,162 @@ +/* + * (C) Crown Copyright 2022 Met Office + * + * 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. + */ + +#pragma once + +/// @file SquareMatrix.h +/// +/// This file contains classes and functions for working with small square matrices. + +#include + +#include "atlas/util/Point.h" + +namespace atlas { + +/// @brief Simple 2x2 matrix with static memory allocation. +/// +/// @details Matrix class which can be used in conjunction with other matrices +/// and Point2 classes. + +template +class SquareMatrix2 { + // Data storage. + std::array, 2> data_{}; + +public: + + /// @brief Default constructor. + SquareMatrix2() = default; + + /// @brief List constructor. + SquareMatrix2(std::initializer_list> list) { + // Get pointer to first element of data_. + Value* elemPtr = data_.front().data(); + for (const std::initializer_list& subList : list) { + for (const Value& elem : subList) { + *elemPtr++ = elem; + } + } + } + + /// @brief std::array constructor. + SquareMatrix2(const std::array, 2>& arr) : data_{arr} {} + + /// @brief return data array. + std::array, 2>& data() {return data_;} + + /// @brief return const data array. + const std::array, 2> data() const {return data_;} + + /// @brief get row. + const std::array& operator[](size_t i) const {return data_[i];} + + /// @brief set row. + std::array& operator[](size_t i) {return data_[i];} + + /// @brief Determinant of matrix. + Value det() const { + return data_[0][0] * data_[1][1] - data_[0][1] * data_[1][0]; + } + + /// @brief matrix-scalar multiplication. + SquareMatrix2 operator*(Value a) const { + return SquareMatrix2{{data_[0][0] * a, data_[0][1] * a}, + {data_[1][0] * a, data_[1][1] * a}}; + } + + /// @brief matrix-vector multiplication. + Point2 operator*(const Point2& x) const { + return Point2{data_[0][0] * x[0] + data_[0][1] * x[1], data_[1][0] * x[0] + data_[1][1] * x[1]}; + } + + /// @brief matrix-matrix multiplication. + SquareMatrix2 operator*(const SquareMatrix2& B) const { + return SquareMatrix2{{data_[0][0] * B[0][0] + data_[0][1] * B[1][0], + data_[0][0] * B[0][1] + data_[0][1] * B[1][1]}, + {data_[1][0] * B[0][0] + data_[1][1] * B[1][0], + data_[1][0] * B[0][1] + data_[1][1] * B[1][1]}}; + } + + /// @brief Inverse matrix. + SquareMatrix2 inverse() const { + return SquareMatrix2{{data_[1][1], -data_[0][1]}, + {-data_[1][0], data_[0][0]}} * (1. / det()); + } + + /// @brief Get signed elements of matrix (i.e., 0, +1 or -1). + SquareMatrix2 sign() const { + const Value smallNumber = det() * std::numeric_limits::epsilon(); + const auto signValue = [&](Value number) -> Value { + return std::abs(number) < smallNumber ? 0. : number < 0. ? -1. : 1.; + }; + return SquareMatrix2{{signValue(data_[0][0]), signValue(data_[0][1])}, + {signValue(data_[1][0]), signValue(data_[1][1])}}; + } + + + +}; + +/// @brief ostream insertion operator. +template +inline std::ostream& operator<<(std::ostream& out, const SquareMatrix2& A) { + + out << "{{" << A[0][0] << "," << A[0][1] << "},{" + << A[1][0] << "," << A[1][1] << "}}"; + + return out; +} + + +/// @brief Jacobian class for a 2D vector field. +/// +/// @details Includes a finite difference constructor which takes three f(x) +/// values arranged as follows: +/// +/// ^ +/// | *f(X0, X1 + dx1) +/// | +/// x1 | +/// | +/// | *f(X0, X1) *f(X0 + dx0, X1) +/// +----------------------------> +/// x0 +class JacobianXY : public SquareMatrix2 { +public: +using SquareMatrix2::SquareMatrix2; + + + /// @brief Converting copy constructor. + JacobianXY(const SquareMatrix2& mat) : SquareMatrix2(mat) {}; + + /// @brief Converting move constructor. + JacobianXY(SquareMatrix2&& mat) : SquareMatrix2(mat) {}; + + /// @brief Converting copy assignment. + JacobianXY& operator=(const SquareMatrix2& mat) { + this->data() = mat.data(); + return *this; + } + + /// @brief Converting move assignment. + JacobianXY& operator=(SquareMatrix2&& mat) { + this->data() = std::move(mat.data()); + return *this; + } + + /// @brief Finite difference constructor. + /// + /// @param f00 = f(X0 , X1 ) + /// @param f10 = f(X0 + dx0, X1 ) + /// @param f01 = f(X0, , X1 + dx1) + JacobianXY(const Point2& f00, const Point2& f10, const Point2& f01, double dx0 = 1., double dx1 = 1.) : + SquareMatrix2{{(f10[0] - f00[0]) / dx0, (f01[0] - f00[0]) / dx1}, + {(f10[1] - f00[1]) / dx0, (f01[1] - f00[1]) / dx1}} {} +}; + +} // namespace atlas diff --git a/src/tests/mesh/test_cubedsphere_meshgen.cc b/src/tests/mesh/test_cubedsphere_meshgen.cc index 117bf4b24..698b39ccd 100644 --- a/src/tests/mesh/test_cubedsphere_meshgen.cc +++ b/src/tests/mesh/test_cubedsphere_meshgen.cc @@ -366,7 +366,7 @@ CASE("cubedsphere_dual_mesh_test") { const auto targetPartitioner = grid::MatchingPartitioner(sourceMesh); // Set target grid, mesh and functionspace. - const auto targetGrid = Grid("CS-LFR-C-48"); + const auto targetGrid = Grid("CS-LFR-C-24"); const auto targetMesh = MeshGenerator("cubedsphere_dual", util::Config("halo", 3)).generate(targetGrid, targetPartitioner); const auto targetFunctionSpace = functionspace::NodeColumns(targetMesh); auto targetField = @@ -417,7 +417,7 @@ CASE("cubedsphere_dual_mesh_test") { const idx_t nLevels = 5; // Set grid, mesh and functionspace. - const auto grid = Grid("CS-LFR-C-48"); + const auto grid = Grid("CS-LFR-C-24"); const auto mesh = MeshGenerator("cubedsphere_dual", util::Config("halo", 3)).generate(grid); const auto functionSpace = functionspace::CellColumns(mesh); diff --git a/src/tests/projection/CMakeLists.txt b/src/tests/projection/CMakeLists.txt index b44a1ef90..c2c744282 100644 --- a/src/tests/projection/CMakeLists.txt +++ b/src/tests/projection/CMakeLists.txt @@ -9,7 +9,6 @@ foreach(test test_bounding_box test_projection_LAEA - test_projection_cubed_sphere test_projection_variable_resolution test_rotation ) @@ -28,3 +27,11 @@ endif() if( HAVE_FCTEST ) add_fctest( TARGET atlas_fctest_projection SOURCES fctest_projection.F90 LINKER_LANGUAGE Fortran LIBS atlas_f ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) endif() + +ecbuild_add_test( TARGET atlas_test_cubedsphere_projection + MPI 6 + CONDITION eckit_HAVE_MPI AND MPI_SLOTS GREATER_EQUAL 6 + SOURCES test_cubedsphere_projection.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) diff --git a/src/tests/projection/test_projection_cubed_sphere.cc b/src/tests/projection/test_cubedsphere_projection.cc similarity index 57% rename from src/tests/projection/test_projection_cubed_sphere.cc rename to src/tests/projection/test_cubedsphere_projection.cc index f19ea633e..0e6fc5af6 100644 --- a/src/tests/projection/test_projection_cubed_sphere.cc +++ b/src/tests/projection/test_cubedsphere_projection.cc @@ -1,33 +1,96 @@ /* - * (C) Crown Copyright 2021 MetOffice + * (C) Crown Copyright 2022 Met Office * * This software is licensed under the terms of the Apache Licence Version 2.0 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. - * In applying this licence, ECMWF does not waive the privileges and immunities - * granted to it by virtue of its status as an intergovernmental organisation - * nor does it submit to any jurisdiction. */ -#include -#include +#include "atlas/array/MakeView.h" #include "atlas/grid.h" -#include "atlas/grid/Grid.h" -#include "atlas/projection/Projection.h" +#include "atlas/grid/CubedSphereGrid.h" +#include "atlas/mesh.h" +#include "atlas/meshgenerator.h" +#include "atlas/projection/detail/CubedSphereProjectionBase.h" +#include "atlas/output/Gmsh.h" #include "atlas/util/Constants.h" -#include "atlas/util/Point.h" - - -#include "atlas/grid/Tiles.h" -#include "atlas/grid/detail/tiles/FV3Tiles.h" -#include "atlas/grid/detail/tiles/LFRicTiles.h" -#include "atlas/grid/detail/tiles/Tiles.h" +#include "atlas/util/SquareMatrix.h" #include "tests/AtlasTestEnvironment.h" - namespace atlas { namespace test { +// Small number relative to 360. +constexpr double epsilon = std::numeric_limits::epsilon() * 360.; + +void testProjection(const std::string& gridType, const std::string& meshType, + const std::string& outputID) { + + // Create grid. + const auto grid = CubedSphereGrid(gridType); + + // Get projection. + const auto& csProjection = grid.cubedSphereProjection(); + + // Output tile centres and Jacobains. + for (size_t i = 0; i < 6; ++i) { + Log::info() << "Tile " << i << std::endl + << "Centre:" << std::endl + << csProjection.getCubedSphereTiles().tileCentre(i) << std::endl + << "Jacobian:" << std::endl + << csProjection.getCubedSphereTiles().tileJacobian(i) << std::endl << std::endl; + } + + // Create mesh. + auto mesh = MeshGenerator(meshType, util::Config("halo", 3)).generate(grid); + + // Output mesh. + const auto gmshConfig = + util::Config("coordinates", "xy") | util::Config("ghost", true) | util::Config("info", true); + auto gmsh = output::Gmsh(outputID + "_before.msh", gmshConfig); + gmsh.write(mesh); + + + + // "Hack" mesh xy coordinates. + auto xyView = array::make_view(mesh.nodes().xy()); + const auto tijView = array::make_view(mesh.nodes().field("tij")); + for (idx_t i = 0; i < mesh.nodes().size(); ++i) { + const auto t = tijView(i, 0); + + const auto xy = PointXY(xyView(i, XX), xyView(i, YY)); + const auto alphabeta = csProjection.xy2alphabeta(xy, t); -//----------------------------------------------------------------------------- + // Inverse function is degenerate when abs(alpha) == abs(beta) and + // abs(alpha) > 45. + const bool degenerate = std::abs(alphabeta[0]) > 45. && + approx_eq(std::abs(alphabeta[0]), std::abs(alphabeta[1])); + + // Check inverse function. + if (!degenerate) { + const auto newXy = csProjection.alphabeta2xy(alphabeta, t); + EXPECT_APPROX_EQ(xy[XX], newXy[XX], epsilon); + EXPECT_APPROX_EQ(xy[YY], newXy[YY], epsilon); + } + + // overwrite mesh xy field. + xyView(i, XX) = alphabeta[0]; + xyView(i, YY) = alphabeta[1]; + + } + + // Output mesh updated mesh. + gmsh = output::Gmsh(outputID + "_after.msh", gmshConfig); + gmsh.write(mesh); + +} + + + +CASE("cubedsphere_xy_to_alphabeta_test") { + + testProjection("CS-LFR-13", "cubedsphere", "cs_primal"); + testProjection("CS-LFR-13", "cubedsphere_dual", "cs_dual"); + +} CASE("test_tiles") { int resolution(2); @@ -131,7 +194,6 @@ CASE("test_projection_cubedsphere_xy_latlon") { } } -//----------------------------------------------------------------------------- } // namespace test } // namespace atlas