From 20094caed6a00b1adbc09692bd1254f1d18ccfe5 Mon Sep 17 00:00:00 2001 From: Slavko Brdar Date: Wed, 8 Mar 2023 11:57:19 +0100 Subject: [PATCH 1/3] New algorithm for ConvexSphericalPolygon intersection; interpolation fails --- src/atlas/util/ConvexSphericalPolygon.cc | 547 +++++++---------- src/atlas/util/ConvexSphericalPolygon.h | 44 +- src/tests/util/test_convexsphericalpolygon.cc | 559 +++++++++++++----- 3 files changed, 643 insertions(+), 507 deletions(-) diff --git a/src/atlas/util/ConvexSphericalPolygon.cc b/src/atlas/util/ConvexSphericalPolygon.cc index fe266f87f..40b3584c7 100644 --- a/src/atlas/util/ConvexSphericalPolygon.cc +++ b/src/atlas/util/ConvexSphericalPolygon.cc @@ -19,9 +19,9 @@ #include "atlas/util/ConvexSphericalPolygon.h" #include "atlas/util/CoordinateEnums.h" #include "atlas/util/NormaliseLongitude.h" +#include "atlas/library/FloatingPointExceptions.h" -#define DEBUG_OUTPUT 0 -#define DEBUG_OUTPUT_DETAIL 0 +// #define DEBUG_OUTPUT 1 namespace atlas { namespace util { @@ -32,14 +32,6 @@ namespace { constexpr double EPS = std::numeric_limits::epsilon(); constexpr double EPS2 = EPS * EPS; -constexpr double TOL = 1.e4 * EPS; // two points considered "same" -constexpr double TOL2 = TOL * TOL; - -enum IntersectionType -{ - NO_INTERSECT = -100, - OVERLAP -}; double distance2(const PointXYZ& p1, const PointXYZ& p2) { double dx = p2[0] - p1[0]; @@ -52,38 +44,35 @@ double norm2(const PointXYZ& p) { return p[0] * p[0] + p[1] * p[1] + p[2] * p[2]; } -bool approx_eq(const double& v1, const double& v2, const double& tol) { - return std::abs(v1 - v2) < tol; +bool approx_eq(const double& v1, const double& v2) { + return std::abs(v1 - v2) <= EPS; } -bool approx_eq(const PointXYZ& v1, const PointXYZ& v2, const double& tol) { +bool approx_eq(const PointXYZ& v1, const PointXYZ& v2) { //return approx_eq( v1[0], v2[0], t ) && approx_eq( v1[1], v2[1], t ) && approx_eq( v1[2], v2[2], t ); - return distance2(v1, v2) < tol * tol; -} - -bool approx_eq_null(const PointXYZ& v1, const double& tol) { - //return approx_eq( v1[0], 0., t ) && approx_eq( v1[1], 0., t ) && approx_eq( v1[2], 0., t ); - double n = v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]; - return n < tol * tol; + return distance2(v1, v2) <= EPS2; } void lonlat2xyz(const PointLonLat& lonlat, PointXYZ& xyz) { eckit::geometry::Sphere::convertSphericalToCartesian(1., lonlat, xyz); } -void xyz2lonlat(const PointXYZ xyz, PointLonLat& lonlat) { +void xyz2lonlat(const PointXYZ& xyz, PointLonLat& lonlat) { eckit::geometry::Sphere::convertCartesianToSpherical(1., xyz, lonlat); } -double norm_max(const PointXYZ& p, const PointXYZ& q) { - double n01 = std::max(std::abs(p[0] - q[0]), std::abs(p[1] - q[1])); - return std::max(n01, std::abs(p[2] - q[2])); +PointLonLat xyz2lonlat(const PointXYZ& xyz) { + PointLonLat lonlat; + eckit::geometry::Sphere::convertCartesianToSpherical(1., xyz, lonlat); + return lonlat; } -template +#if 0 +// NOTE: StackVector is not used +template struct StackVector { private: - using Wrapped = std::array; + using Wrapped = std::array; public: using reference = typename Wrapped::reference; @@ -101,7 +90,7 @@ struct StackVector { #endif void push_back(const T& value) { wrapped_[size_++] = value; - ATLAS_ASSERT(size_ < ConvexSphericalPolygon::MAX_SIZE); + ATLAS_ASSERT(size_ < MAX_SIZE); } size_t size() const { return size_; } @@ -109,38 +98,7 @@ struct StackVector { size_t size_{0}; Wrapped wrapped_; }; - -struct PolygonEdgeIntersection { - static constexpr int BEGIN = 1; - static constexpr int END = 2; - static constexpr int INSIDE = 3; - - PolygonEdgeIntersection(const ConvexSphericalPolygon& polygon, int edge_index, const PointXYZ& point) { - auto matches = [](const PointXYZ& p1, const PointXYZ& p2) { - return (distance2(p1, p2) < 1e-16); - // We would like this to be TOL2 instead, but that gives bad results - }; - - ATLAS_ASSERT(edge_index >= 0); - ATLAS_ASSERT(edge_index < polygon.size()); - - const int node_index = edge_index; - - if (matches(point, polygon[node_index])) { - location = BEGIN; - } - else if (matches(point, polygon[polygon.next(node_index)])) { - location = END; - } - else { - location = INSIDE; - } - } - bool isPointAtBegin() const { return location == BEGIN; } - bool isPointAtEnd() const { return location == END; } - bool isPointInside() const { return location == INSIDE; } - int location; -}; +#endif } // namespace @@ -153,13 +111,14 @@ ConvexSphericalPolygon::ConvexSphericalPolygon(const PointLonLat points[], size_ size_t isp = 1; for (size_t i = 1; i < size_ - 1; ++i) { lonlat2xyz(points[i], sph_coords_[isp]); - if (approx_eq(sph_coords_[isp], sph_coords_[isp - 1], TOL)) { +// Log::info() << " d : " << PointXYZ::distance(sph_coords_[isp], sph_coords_[isp - 1]) << std::endl; + if (approx_eq(sph_coords_[isp], sph_coords_[isp - 1])) { continue; } ++isp; } lonlat2xyz(points[size_ - 1], sph_coords_[isp]); - if (approx_eq(sph_coords_[isp], sph_coords_[0], TOL) or approx_eq(sph_coords_[isp], sph_coords_[isp - 1], TOL)) { + if (approx_eq(sph_coords_[isp], sph_coords_[0]) or approx_eq(sph_coords_[isp], sph_coords_[isp - 1])) { } else { ++isp; @@ -178,13 +137,13 @@ ConvexSphericalPolygon::ConvexSphericalPolygon(const PointXYZ points[], size_t s size_t isp = 1; for (size_t i = 1; i < size_ - 1; ++i) { sph_coords_[isp] = points[i]; - if (approx_eq(sph_coords_[isp], sph_coords_[isp - 1], TOL)) { + if (approx_eq(sph_coords_[isp], sph_coords_[isp - 1])) { continue; } ++isp; } sph_coords_[isp] = points[size_ - 1]; - if (approx_eq(sph_coords_[isp], sph_coords_[0], TOL) or approx_eq(sph_coords_[isp], sph_coords_[isp - 1], TOL)) { + if (approx_eq(sph_coords_[isp], sph_coords_[0]) or approx_eq(sph_coords_[isp], sph_coords_[isp - 1])) { } else { ++isp; @@ -196,74 +155,8 @@ ConvexSphericalPolygon::ConvexSphericalPolygon(const PointXYZ points[], size_t s } } - -void ConvexSphericalPolygon::simplify() { - ATLAS_ASSERT(size_ < MAX_SIZE); - if (size_ < 3) { - size_ = 0; - valid_ = false; - return; - } - idx_t isp = 0; - idx_t i = 0; - idx_t j; - idx_t k; - bool search_3pts = true; - auto& points = sph_coords_; - for (; i < size_ && search_3pts; ++i) { - const PointXYZ& P0 = points[i]; - for (j = i + 1; j < size_ && search_3pts; ++j) { - const PointXYZ& P1 = points[j]; - if (approx_eq(P0, P1, 1.e-10)) { - continue; - } - for (k = j + 1; k < size_ && search_3pts; ++k) { - const PointXYZ& P2 = points[k]; - if (approx_eq(P1, P2, TOL) or approx_eq(P0, P2, TOL)) { - continue; - } - if (GreatCircleSegment{P0, P1}.inLeftHemisphere(P2, -EPS)) { - sph_coords_[isp++] = P0; - sph_coords_[isp++] = P1; - sph_coords_[isp++] = P2; - search_3pts = false; - } - } - } - } - if (search_3pts) { - valid_ = false; - size_ = 0; - return; - } - for (; k < size_ - 1; ++k) { - if (approx_eq(points[k], sph_coords_[isp - 1], TOL) or - (not GreatCircleSegment{sph_coords_[isp - 2], sph_coords_[isp - 1]}.inLeftHemisphere(points[k], -EPS))) { - continue; - } - sph_coords_[isp] = points[k]; - isp++; - } - const PointXYZ& Pl2 = sph_coords_[isp - 2]; - const PointXYZ& Pl1 = sph_coords_[isp - 1]; - const PointXYZ& P0 = sph_coords_[0]; - const PointXYZ& P = points[size_ - 1]; - if ((not approx_eq(P, P0, EPS)) and (not approx_eq(P, Pl1, EPS)) and - GreatCircleSegment{Pl2, Pl1}.inLeftHemisphere(P, -EPS) and - GreatCircleSegment{Pl1, P}.inLeftHemisphere(P0, -EPS)) { - sph_coords_[isp] = P; - ++isp; - } - size_ = isp; - valid_ = size_ > 2; - - computed_area_ = false; - computed_radius_ = false; - computed_centroid_ = false; -} - void ConvexSphericalPolygon::compute_centroid() const { - const auto triangles = triangulate(radius()); + const auto triangles = triangulate(); area_ = triangles.area(); computed_area_ = true; @@ -288,9 +181,9 @@ bool ConvexSphericalPolygon::validate() { int nni = next(ni); const PointXYZ& P = sph_coords_[i]; const PointXYZ& nextP = sph_coords_[ni]; - ATLAS_ASSERT(std::abs(PointXYZ::dot(P, P) - 1.) < 10. * EPS); - ATLAS_ASSERT(not approx_eq(P, PointXYZ::mul(nextP, -1.), TOL)); - valid_ = valid_ && GreatCircleSegment{P, nextP}.inLeftHemisphere(sph_coords_[nni], -EPS); + ATLAS_ASSERT(std::abs(PointXYZ::dot(P, P) - 1.) < 10 * EPS); + ATLAS_ASSERT(not approx_eq(P, PointXYZ::mul(nextP, -1.))); + valid_ = valid_ && GreatCircleSegment{P, nextP}.inLeftHemisphere(sph_coords_[nni], -0.5*EPS); } } return valid_; @@ -300,8 +193,16 @@ bool ConvexSphericalPolygon::equals(const ConvexSphericalPolygon& plg, const dou if (size_ == 0 and plg.size_ == 0) { return true; } - if ((not plg.valid_) || (not valid_) || size_ != plg.size()) { - Log::info() << " ConvexSphericalPolygon::equals == not compatible\n"; + if (not valid_) { + Log::info() << " ConvexSphericalPolygon::equals : this polygon is not valid\n"; + return false; + } + if (not plg.valid_) { + Log::info() << " ConvexSphericalPolygon::equals : other polygon passed as argument is not valid\n"; + return false; + } + if (size_ != plg.size()) { + Log::info() << " ConvexSphericalPolygon::equals : incompatible sizes: " << size_ << " != " << plg.size() << "\n"; return false; } int offset = 0; @@ -313,7 +214,7 @@ bool ConvexSphericalPolygon::equals(const ConvexSphericalPolygon& plg, const dou } } if (offset == size_) { - Log::info() << "ConvexSphericalPolygon::equals == no point equal\n"; + Log::info() << "ConvexSphericalPolygon::equals : no point equal\n"; return false; } @@ -321,82 +222,38 @@ bool ConvexSphericalPolygon::equals(const ConvexSphericalPolygon& plg, const dou int idx = (offset + j) % size_; auto dist2 = distance2(plg.sph_coords_[j], sph_coords_[idx]); if (dist2 > le2) { - Log::info() << " ConvexSphericalPolygon::equals == point distance " << std::sqrt(dist2) << "\n"; + Log::info() << " ConvexSphericalPolygon::equals : point distance " << std::sqrt(dist2) << " < " << le2 << "(= "<< deg_prec << " deg)" << "\n"; return false; } } return true; } -// note: unit sphere! -// I. Todhunter (1886), Paragr. 99 -ConvexSphericalPolygon::SubTriangles ConvexSphericalPolygon::triangulate(const double cell_radius) const { +// cf. Folke Eriksson, "On the Measure of Solid Angles", Mathematics Magazine, Vol. 63, No. 3, pp. 184-187 (1990) +ConvexSphericalPolygon::SubTriangles ConvexSphericalPolygon::triangulate() const { SubTriangles triangles; if (size_ < 3) { return triangles; } - size_t itri{0}; - if (cell_radius < 1.e-6) { // plane area - for (int i = 1; i < size_ - 1; i++) { - const PointXYZ pl = sph_coords_[i] - sph_coords_[0]; - const PointXYZ pr = sph_coords_[i + 1] - sph_coords_[0]; - triangles[itri].centroid = PointXYZ::normalize(sph_coords_[0] + sph_coords_[i] + sph_coords_[i + 1]); - triangles[itri].area = 0.5 * PointXYZ::norm(PointXYZ::cross(pl, pr)); - ++itri; - } - } - else { // spherical area - const PointXYZ& a = sph_coords_[0]; - for (size_t i = 1; i < size_ - 1; i++) { - const PointXYZ& b = sph_coords_[i]; - const PointXYZ& c = sph_coords_[i + 1]; - auto ab = PointXYZ::cross(a, b); - auto bc = PointXYZ::cross(b, c); - auto ca = PointXYZ::cross(c, a); - const double ab_norm = PointXYZ::norm(ab); - const double bc_norm = PointXYZ::norm(bc); - const double ca_norm = PointXYZ::norm(ca); - if (ab_norm < EPS or bc_norm < EPS or ca_norm < EPS) { - continue; - } - double abc = -PointXYZ::dot(ab, bc) / (ab_norm * bc_norm); - double bca = -PointXYZ::dot(bc, ca) / (bc_norm * ca_norm); - double cab = -PointXYZ::dot(ca, ab) / (ca_norm * ab_norm); - if (abc <= -1.) { - abc = M_PI; - } - else if (abc < 1.) { - abc = std::acos(abc); - } - else { - abc = 0.; - } - if (bca <= -1.) { - bca = M_PI; - } - else if (bca < 1.) { - bca = std::acos(bca); - } - else { - bca = 0.; - } - if (cab <= -1.) { - cab = M_PI; - } - else if (cab < 1.) { - cab = std::acos(cab); - } - else { - cab = 0.; - } - triangles[itri].centroid = PointXYZ::normalize(a + b + c); - triangles[itri].area = abc + bca + cab - M_PI; - ++itri; - } + const PointXYZ& a = sph_coords_[0]; + for (size_t i = 1; i < size_ - 1; i++) { + const PointXYZ& b = sph_coords_[i]; + const PointXYZ& c = sph_coords_[i + 1]; + triangles[itri].centroid = PointXYZ::normalize(a + b + c); +// if (PointXYZ::distance(a, b) + PointXYZ::distance(b, c) < 1e-10) { +// triangles[itri].area = 0.5 * PointXYZ::norm(PointXYZ::cross(b - a, c - b)); +// } +// else { + auto abc = PointXYZ::dot(a, b) + PointXYZ::dot(b, c) + PointXYZ::dot(c, a); + auto a_bc = PointXYZ::dot(a, PointXYZ::cross(b, c)); + triangles[itri].area = 2. * std::atan(std::abs(a_bc) / (1. + abc)); +// } + ++itri; } triangles.size() = itri; return triangles; + } @@ -408,95 +265,93 @@ double ConvexSphericalPolygon::SubTriangles::area() const { return area; } -// @return lowest point id of this polygon's segment intersecting [s1,s2)) -int ConvexSphericalPolygon::intersect(const int start, const GreatCircleSegment& s, PointXYZ& I) const { - for (int i = start; i < size_; i++) { - const int id0 = i; - const int id1 = next(i); - const GreatCircleSegment p(sph_coords_[id0], sph_coords_[id1]); - I = s.intersect(p); - if (I[0] == 0 && I[1] == 0 && I[2] == 0) { - // intersection not on [p1,p2) - continue; - } - if (I[0] == 1 && I[1] == 1) { - return OVERLAP; - } - return id0; - } - return NO_INTERSECT; -} - - -void ConvexSphericalPolygon::clip(const GreatCircleSegment& great_circle) { +void ConvexSphericalPolygon::clip(const GreatCircleSegment& great_circle, std::ostream* out, double pointsSameEPS) { ATLAS_ASSERT(valid_); - ATLAS_ASSERT(distance2(great_circle.first(), great_circle.second()) > TOL2); - + ATLAS_ASSERT(not approx_eq(great_circle.first(), great_circle.second())); + std::vector clipped_sph_coords; + clipped_sph_coords.reserve(ConvexSphericalPolygon::MAX_SIZE); auto invalidate_this_polygon = [&]() { size_ = 0; valid_ = false; area_ = 0.; }; - - // Count and mark all vertices to be possibly considered in clipped polygon - StackVector vertex_in(size_); - int num_vertices_in = 0; +#if DEBUG_OUTPUT + int add_point_num = 0; +#endif + bool first_in = great_circle.inLeftHemisphere(sph_coords_[0], -1.5 * EPS, out); for (int i = 0; i < size_; i++) { - vertex_in[i] = great_circle.inLeftHemisphere(sph_coords_[i], -10. * TOL); - num_vertices_in += vertex_in[i] ? 1 : 0; - } - - PointXYZ i1; - const int f1 = intersect(0, great_circle, i1); - const bool segment_only_touches_last_point = (f1 == size_ - 1); - if (f1 == OVERLAP || f1 == NO_INTERSECT || segment_only_touches_last_point) { - if (num_vertices_in < 3) { - invalidate_this_polygon(); - } - return; - } - PolygonEdgeIntersection intersection_1(*this, f1, i1); - - PointXYZ i2; // second intersection point - auto start2 = [&]() { return f1 + 1 + (intersection_1.isPointAtEnd() ? 1 : 0); }; - const int f2 = intersect(start2(), great_circle, i2); - if (f2 == OVERLAP || f2 == NO_INTERSECT) { - if (num_vertices_in < 3) { - invalidate_this_polygon(); + int in = (i+1) % size_; + bool second_in = great_circle.inLeftHemisphere(sph_coords_[in], -1.5 * EPS, out); +#if DEBUG_OUTPUT + if (out) { + out->precision(18); + *out << " ** first: " << xyz2lonlat(sph_coords_[i]) << ", in ? " << first_in << std::endl; + *out << " second: " << xyz2lonlat(sph_coords_[in]) << ", in ? " << second_in << std::endl; } - return; - } - PolygonEdgeIntersection intersection_2(*this, f2, i2); - - // Create new vector of clipped coordinates - StackVector clipped_sph_coords; - { - auto keep_vertex = [&](int index) { clipped_sph_coords.push_back(sph_coords_[index]); }; - auto insert_point = [&](const PointXYZ& p) { clipped_sph_coords.push_back(p); }; - - for (int i = 0; i <= f1; i++) { - if (vertex_in[i]) { - keep_vertex(i); +#endif + if (first_in and second_in) { + clipped_sph_coords.emplace_back(sph_coords_[in]); +#if DEBUG_OUTPUT + if (out) { + *out << " ((" << ++add_point_num << ")) both_in add second: " << xyz2lonlat(sph_coords_[in]) << std::endl; } +#endif } - if ((not vertex_in[f1] and intersection_1.isPointAtBegin()) or - (not vertex_in[next(f1)] and intersection_1.isPointAtEnd()) or intersection_1.isPointInside()) { - insert_point(i1); + else if (not first_in and not second_in) { + // continue to update first_in } - for (int i = f1 + 1; i <= f2; i++) { - if (vertex_in[i]) { - keep_vertex(i); + else { + const GreatCircleSegment segment(sph_coords_[i], sph_coords_[in]); + PointXYZ ip = great_circle.intersect(segment, out, pointsSameEPS); +#if DEBUG_OUTPUT + if (out) { + *out << " ip : " << xyz2lonlat(ip) << std::endl; } - } - if ((not vertex_in[f2] and intersection_2.isPointAtBegin()) or - (not vertex_in[next(f2)] and intersection_2.isPointAtEnd()) or intersection_2.isPointInside()) { - insert_point(i2); - } - for (int i = f2 + 1; i < size_; i++) { - if (vertex_in[i]) { - keep_vertex(i); +#endif + if (ip[0] == 1 and ip[1] == 1 and ip[2] == 1) { + // consider the segments parallel +#if DEBUG_OUTPUT + if (out) { + *out << " ((" << ++add_point_num << ")) ip=(1,1,1) add second: " + << xyz2lonlat(sph_coords_[in]) << std::endl; + } +#endif + clipped_sph_coords.emplace_back(sph_coords_[in]); + first_in = second_in; + continue; + } + if (second_in) { + int inn = (in+1) % size_; + const GreatCircleSegment segment_n(sph_coords_[in], sph_coords_[inn]); + if (segment.inLeftHemisphere(ip, -1.5 * EPS, out) and + segment_n.inLeftHemisphere(ip, -1.5 * EPS, out) and + (PointXYZ::distance(ip, sph_coords_[in]) > pointsSameEPS)) { + clipped_sph_coords.emplace_back(ip); +#if DEBUG_OUTPUT + if (out) { + *out << " ((" << ++add_point_num << ")) second_in add ip: " << xyz2lonlat(ip) << std::endl; + } +#endif + } + clipped_sph_coords.emplace_back(sph_coords_[in]); +#if DEBUG_OUTPUT + if (out) { + *out << " ((" << ++add_point_num << ")) second_in add second: " << xyz2lonlat(sph_coords_[in]) << std::endl; + } +#endif + } + else { + if (PointXYZ::distance(ip, sph_coords_[i]) > pointsSameEPS) { + clipped_sph_coords.emplace_back(ip); +#if DEBUG_OUTPUT + if (out) { + *out << " ((" << ++add_point_num << ")) first_in add ip: " << xyz2lonlat(ip) << std::endl; + } +#endif + } } } + first_in = second_in; } // Update polygon @@ -516,22 +371,68 @@ void ConvexSphericalPolygon::clip(const GreatCircleSegment& great_circle) { // intersect a polygon with this polygon // @param[in] pol clipping polygon // @param[out] intersecting polygon -ConvexSphericalPolygon ConvexSphericalPolygon::intersect(const ConvexSphericalPolygon& plg) const { - ConvexSphericalPolygon intersection = *this; +ConvexSphericalPolygon ConvexSphericalPolygon::intersect(const ConvexSphericalPolygon& plg, std::ostream* out, double pointsSameEPS) const { + + bool fpe_disabled = atlas::library::disable_floating_point_exception(FE_INVALID); + auto restore_fpe = [fpe_disabled] { + if (fpe_disabled) { + atlas::library::enable_floating_point_exception(FE_INVALID); + } + }; + + // the larger area polygon is the intersector + ConvexSphericalPolygon intersection; + ConvexSphericalPolygon intersector; + std::string intor_id = "P1"; + std::string inted_id = "P2"; + + bool this_intersector = (area() > plg.area()); + if (approx_eq(area(), plg.area())) { + PointXYZ dc = centroid() - plg.centroid(); + if (dc[0] > 0. or (dc[0] == 0. and (dc[1] > 0. or (dc[1] == 0. and dc[2] > 0.)))) { + this_intersector = true; + } + } + if (this_intersector) { + intersector = *this; + intersection = plg; + } + else { + intor_id = "P2"; + inted_id = "P1"; + intersector = plg; + intersection = *this; + } +#if DEBUG_OUTPUT + if (out) { + *out << inted_id << " : "; + print(*out); + *out << std::endl << intor_id << " : "; + plg.print(*out); + *out << std::endl; + } +#endif if (intersection.valid_) { - for (size_t i = 0; i < plg.size_; i++) { - const PointXYZ& s1 = plg.sph_coords_[i]; - const PointXYZ& s2 = plg.sph_coords_[(i != plg.size_ - 1) ? i + 1 : 0]; - intersection.clip(GreatCircleSegment(s1, s2)); + for (size_t i = 0; i < intersector.size_; i++) { + const PointXYZ& s1 = intersector.sph_coords_[i]; + const PointXYZ& s2 = intersector.sph_coords_[(i != intersector.size_ - 1) ? i + 1 : 0]; +#if DEBUG_OUTPUT + if (out) { + *out << std::endl << "Clip with [" << intor_id << "_" << i << ", " << intor_id << "_" + << (i+1) % intersector.size_ << "]" << std::endl; + } +#endif + intersection.clip(GreatCircleSegment(s1, s2), out, pointsSameEPS); if (not intersection.valid_) { + restore_fpe(); return intersection; } } } - intersection.simplify(); intersection.computed_area_ = false; intersection.computed_radius_ = false; intersection.computed_centroid_ = false; + restore_fpe(); return intersection; } @@ -548,72 +449,64 @@ void ConvexSphericalPolygon::print(std::ostream& out) const { out << "}"; } +std::string ConvexSphericalPolygon::json(int precision) const { + std::stringstream ss; + if( precision ) { + ss.precision(16); + } + ss << "["; + for (size_t i = 0; i < size(); ++i) { + if (i > 0) { + ss << ","; + } + PointLonLat ip_ll; + xyz2lonlat(sph_coords_[i], ip_ll); + ss << "[" << ip_ll.lon() << "," << ip_ll.lat() << "]"; + } + ss << "]"; + return ss.str(); +} + + double ConvexSphericalPolygon::compute_radius() const { double radius{0.}; if (valid_) { - PointXYZ centroid; - centroid = sph_coords_[0]; - size_t isp = 1; - for (size_t i = 1; i < size_; ++i) { - if (approx_eq(sph_coords_[isp], sph_coords_[isp - 1], TOL)) { - continue; - } - centroid = centroid + sph_coords_[isp]; - ++isp; + if (not computed_centroid_) { + compute_centroid(); } - centroid = PointXYZ::div(centroid, PointXYZ::norm(centroid)); - for (size_t i = 0; i < size_; ++i) { - radius = std::max(radius, PointXYZ::distance(sph_coords_[i], centroid)); + radius = std::max(radius, PointXYZ::distance(sph_coords_[i], centroid_)); } } return radius; } -bool ConvexSphericalPolygon::GreatCircleSegment::contains(const PointXYZ& p) const { - /* - * @brief Point-on-segment test on great circle segments - * @param[in] P given point in (x,y,z) coordinates - * @return - */ - constexpr double eps_large = 1.e3 * EPS; - - // Case where p is one of the endpoints - double pp1n2 = distance2(p, p1_); - double pp2n2 = distance2(p, p2_); - if (pp1n2 < EPS2 or pp2n2 < EPS2) { - return true; - } - - PointXYZ p12 = cross(); - double p12n2 = norm2(p12); - double p12n = std::sqrt(p12n2); - p12 /= p12n; - if (std::abs(PointXYZ::dot(p, p12)) > eps_large) { - return false; - } - double pp = PointXYZ::distance(p1_, p2_); - double pp1 = PointXYZ::distance(p, p1_); - double pp2 = PointXYZ::distance(p, p2_); - return (std::min(pp - pp1, pp - pp2) > -eps_large); -} - -PointXYZ ConvexSphericalPolygon::GreatCircleSegment::intersect(const GreatCircleSegment& p) const { +// 'this' great circle's intersection with the segment 'p': [p.first(), p.second()) +PointXYZ ConvexSphericalPolygon::GreatCircleSegment::intersect(const GreatCircleSegment& p, std::ostream* out, double pointsSameEPS) const { const auto& s = *this; PointXYZ sp = PointXYZ::cross(s.cross(), p.cross()); double sp_norm = PointXYZ::norm(sp); - if (sp_norm > EPS) { + bool gcircles_distinct = (sp_norm > EPS); +#if DEBUG_OUTPUT + if (out) { + *out << " Great circles distinct ? " << sp_norm << " > " << EPS << " ? " + << gcircles_distinct << std::endl; + } +#endif + if (gcircles_distinct) { sp /= sp_norm; - if (p.contains(sp)) { + auto sp_2 = sp * -1.; + double d = distance2(p.first(), p.second()); + double d1 = std::max(distance2(sp, p.first()), distance2(sp, p.second())); + double d2 = std::max(distance2(sp_2, p.first()), distance2(sp_2, p.second())); + if (d1 < d2) { return sp; } - sp *= -1.; - if (p.contains(sp)) { - return sp; + else { + return sp_2; } - return PointXYZ(0, 0, 0); } else { return PointXYZ(1, 1, 1); diff --git a/src/atlas/util/ConvexSphericalPolygon.h b/src/atlas/util/ConvexSphericalPolygon.h index 9a5783fa0..7b3338fd5 100644 --- a/src/atlas/util/ConvexSphericalPolygon.h +++ b/src/atlas/util/ConvexSphericalPolygon.h @@ -16,6 +16,9 @@ #include "atlas/util/Point.h" #include "atlas/util/Polygon.h" +#include "atlas/runtime/Log.h" + +#define DEBUG_OUTPUT 1 namespace atlas { namespace util { @@ -31,17 +34,20 @@ class ConvexSphericalPolygon { class GreatCircleSegment { public: GreatCircleSegment(const PointXYZ& p1, const PointXYZ& p2): p1_(p1), p2_(p2), cross_(PointXYZ::cross(p1, p2)) {} - bool contains(const PointXYZ&) const; - - PointXYZ intersect(const GreatCircleSegment&) const; - // Hemisphere is defined by segment when walking from first() to second() - // Positive offset: distance into left hemisphere, e.g. to exclude segment itself with tolerance - // Negative offset: distance into right hemisphere, e.g. to include segment with tolerance - bool inLeftHemisphere(const PointXYZ& P, const double offset = 0.) const { - return (PointXYZ::dot(cross(), P) > offset); + // For a given segment, the "left" hemisphere is defined on the left of the segment when walking from first() to second() + inline bool inLeftHemisphere(const PointXYZ& P, const double offset = 0., std::ostream* out = NULL) const { +#if DEBUG_OUTPUT + if (out) { + *out << " inLeftHemi: " <= " << offset << " ? " + << (PointXYZ::dot(cross(), P) >= offset) << std::endl; + } +#endif + return (PointXYZ::dot(cross(), P) >= offset); // has to have = included } + PointXYZ intersect(const GreatCircleSegment&, std::ostream* f = NULL, double pointsSameEPS = std::numeric_limits::epsilon()) const; + const PointXYZ& first() const { return p1_; } const PointXYZ& second() const { return p2_; } @@ -96,7 +102,7 @@ class ConvexSphericalPolygon { return radius_; } - ConvexSphericalPolygon intersect(const ConvexSphericalPolygon& pol) const; + ConvexSphericalPolygon intersect(const ConvexSphericalPolygon& pol, std::ostream* f = nullptr, double pointsEqualEPS = std::numeric_limits::epsilon()) const; /* * @brief check if two spherical polygons area equal @@ -107,6 +113,8 @@ class ConvexSphericalPolygon { void print(std::ostream&) const; + std::string json(int precision = 0) const; + friend std::ostream& operator<<(std::ostream& out, const ConvexSphericalPolygon& p) { p.print(out); return out; @@ -136,29 +144,15 @@ class ConvexSphericalPolygon { double compute_radius() const; - SubTriangles triangulate(const double radius) const; + SubTriangles triangulate() const; - void clip(const GreatCircleSegment&); + void clip(const GreatCircleSegment&, std::ostream* f = nullptr, double pointsSameEPS = std::numeric_limits::epsilon()); /* * @return true:polygon is convex */ bool validate(); - /* - * @brief Segment-sph_polygon intersection - * @param[in] s1, s2 segment endpoints in (x,y,z) coordinates - * @param[in] start start with polygon segments [pol[start],pol[start+1]],... - * @param[out] ip intersection point or nullptr - * @return -1: overlap with one of polygon edges, - * 0: no_intersect, - * 1 + (id of this polygon's segment intersecting [s1,s2)): otherwise - */ - int intersect(const int start, const GreatCircleSegment&, PointXYZ& ip) const; - - /// Makes the polygon convex and skips consequtive node that is too close to previous - void simplify(); - private: std::array sph_coords_; mutable PointXYZ centroid_; diff --git a/src/tests/util/test_convexsphericalpolygon.cc b/src/tests/util/test_convexsphericalpolygon.cc index f4e158855..f7897a29e 100644 --- a/src/tests/util/test_convexsphericalpolygon.cc +++ b/src/tests/util/test_convexsphericalpolygon.cc @@ -2,6 +2,8 @@ #include #include #include +#include +#include #include "atlas/util/ConvexSphericalPolygon.h" #include "atlas/util/Geometry.h" @@ -13,6 +15,7 @@ namespace atlas { namespace test { using ConvexSphericalPolygon = util::ConvexSphericalPolygon; +const double EPS = std::numeric_limits::epsilon(); util::ConvexSphericalPolygon make_polygon(const std::initializer_list& list) { return util::ConvexSphericalPolygon{std::vector(list)}; @@ -29,17 +32,112 @@ util::ConvexSphericalPolygon make_polygon() { return util::ConvexSphericalPolygon{}; } + +template +std::string to_json(const It& begin, const It& end, int precision = 0) { + std::stringstream ss; + ss << "[\n"; + size_t size = std::distance(begin,end); + size_t c=0; + for( auto it = begin; it != end; ++it, ++c ) { + ss << " " << it->json(precision); + if( c < size-1 ) { + ss << ",\n"; + } + } + ss << "\n]"; + return ss.str(); +} + +template +std::string to_json(const ConvexSphericalPolygonContainer& polygons, int precision = 0) { + return to_json(polygons.begin(),polygons.end(),precision); +} +std::string to_json(std::initializer_list&& polygons, int precision = 0) { + return to_json(polygons.begin(),polygons.end(),precision); +} + +void check_intersection(const ConvexSphericalPolygon& plg1, const ConvexSphericalPolygon& plg2, const ConvexSphericalPolygon& iplg_sol, double pointsSameEPS = 5.e6 * EPS, std::ostream* out = nullptr) { + auto iplg = plg1.intersect(plg2, out, pointsSameEPS); + Log::info().indent(); + Log::info() << "plg1 area : " << plg1.area() << "\n"; + Log::info() << "plg2 area : " << plg2.area() << "\n"; + Log::info() << "iplg area : " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg1,plg2,iplg},20) << "\n"; + EXPECT(std::min(plg1.area(), plg2.area()) >= iplg.area()); + EXPECT(iplg.equals(iplg_sol, 1.e-8)); + Log::info().unindent(); +} + CASE("test default constructor") { ConvexSphericalPolygon p; EXPECT(bool(p) == false); } +CASE("Size of ConvexSphericalPolygon") { + // This test illustrates that ConvexSphericalPolygon is allocated on the stack completely, + // as sizeof(ConvexSphericalPolygon) includes space for MAX_SIZE coordinates of type PointXYZ + EXPECT(sizeof(PointXYZ) == sizeof(double) * 3); + size_t expected_size = 0; + expected_size += (1 + ConvexSphericalPolygon::MAX_SIZE) * sizeof(PointXYZ); + expected_size += sizeof(size_t); + expected_size += sizeof(bool); + expected_size += 2 * sizeof(double); + EXPECT(sizeof(ConvexSphericalPolygon) >= expected_size); // greater because compiler may add some padding +} + +CASE("analyse intersect") { + const double du = 0.5; + const double dv = 1.1 * EPS; + const double duc = 0.5 * du; + const double sduc = std::sqrt(1. - 0.25 * du * du); + const double dvc = 1. - 0.5 * dv * dv; + const double sdvc = dv * std::sqrt(1. - 0.25 * dv * dv); + PointXYZ s0p0{sduc, -duc, 0.}; + PointXYZ s0p1{sduc, duc, 0.}; + PointXYZ s1p0{dvc * sduc, -dvc * duc, -sdvc}; + PointXYZ s1p1{dvc * sduc, dvc * duc, sdvc}; + + EXPECT_APPROX_EQ(dv, PointXYZ::norm(s0p0 - s1p0), EPS); + EXPECT_APPROX_EQ(du, PointXYZ::norm(s0p0 - s0p1), EPS); + EXPECT_APPROX_EQ(dv, PointXYZ::norm(s0p1 - s1p1), EPS); + + ConvexSphericalPolygon::GreatCircleSegment s1(s0p0, s0p1); + ConvexSphericalPolygon::GreatCircleSegment s2(s1p0, s1p1); + + // analytical solution + PointXYZ Isol{1., 0., 0.}; + + // test "intersection" + PointXYZ I12 = s1.intersect(s2); + PointXYZ I21 = s2.intersect(s1); + EXPECT_APPROX_EQ(std::abs(PointXYZ::norm(I12) - 1.), 0., EPS); + EXPECT_APPROX_EQ(PointXYZ::norm(I12 - Isol), 0., EPS); + EXPECT_APPROX_EQ(std::abs(PointXYZ::norm(I21) - 1.), 0., EPS); + EXPECT_APPROX_EQ(PointXYZ::norm(I21 - Isol), 0., EPS); +} + +CASE("test_json_format") { + auto plg = make_polygon({{0., 45.}, {0., 0.}, {90., 0.}, {90., 45.}}); + EXPECT_EQ(plg.json(), "[[0,45],[0,0],[90,0],[90,45]]"); + + std::vector plg_g = { + make_polygon({{0, 60}, {0, 50}, {40, 60}}), //0 + make_polygon({{0, 60}, {0, 50}, {20, 60}}), + make_polygon({{10, 60}, {10, 50}, {30, 60}}), //3 + }; + Log::info() << to_json(plg_g) << std::endl; + +} + CASE("test_spherical_polygon_area") { auto plg1 = make_polygon({{0., 90.}, {0., 0.}, {90., 0.}}); EXPECT_APPROX_EQ(plg1.area(), M_PI_2); auto plg2 = make_polygon({{0., 45.}, {0., 0.}, {90., 0.}, {90., 45.}}); auto plg3 = make_polygon({{0., 90.}, {0., 45.}, {90., 45.}}); - Log::info() << "area diff: " << plg1.area() - plg2.area() - plg3.area() << std::endl; + Log::info().indent(); + EXPECT(plg1.area() - plg2.area() - plg3.area() <= EPS); + Log::info().unindent(); EXPECT_APPROX_EQ(std::abs(plg1.area() - plg2.area() - plg3.area()), 0, 1e-15); } @@ -50,23 +148,23 @@ CASE("test_spherical_polygon_intersection") { std::array plg_f = {make_polygon({{0, 70}, {0, 60}, {40, 60}, {40, 70}}), make_polygon({{0, 90}, {0, 0}, {40, 0}})}; std::array plg_g = { - make_polygon({{0, 60}, {0, 50}, {40, 60}}), //0 - make_polygon({{0, 60}, {0, 50}, {20, 60}}), - make_polygon({{10, 60}, {10, 50}, {30, 60}}), //3 - make_polygon({{40, 80}, {0, 60}, {40, 60}}), - make_polygon({{0, 80}, {0, 60}, {40, 60}}), //5 - make_polygon({{20, 80}, {0, 60}, {40, 60}}), - make_polygon({{20, 70}, {0, 50}, {40, 50}}), //7 - make_polygon({{0, 90}, {0, 60}, {40, 60}}), - make_polygon({{-10, 80}, {-10, 50}, {50, 80}}), //9 - make_polygon({{0, 80}, {0, 50}, {40, 50}, {40, 80}}), - make_polygon({{0, 65}, {20, 55}, {40, 60}, {20, 65}}), //11 - make_polygon({{20, 65}, {0, 60}, {20, 55}, {40, 60}}), - make_polygon({{10, 63}, {20, 55}, {30, 63}, {20, 65}}), //13 - make_polygon({{20, 75}, {0, 70}, {5, 5}, {10, 0}, {20, 0}, {40, 70}}), - make_polygon({{0, 50}, {0, 40}, {5, 45}}), //15 - make_polygon({{0, 90}, {0, 80}, {20, 0}, {40, 80}}), - make_polygon({{0, 65}, {0, 55}, {40, 65}, {40, 75}}), //17 + make_polygon({{0, 60}, {0, 50}, {40, 60}}), // 0 + make_polygon({{0, 60}, {0, 50}, {20, 60}}), // 1 + make_polygon({{10, 60}, {10, 50}, {30, 60}}), // 2 + make_polygon({{40, 80}, {0, 60}, {40, 60}}), // 3 + make_polygon({{0, 80}, {0, 60}, {40, 60}}), // 4 + make_polygon({{20, 80}, {0, 60}, {40, 60}}), // 5 + make_polygon({{20, 70}, {0, 50}, {40, 50}}), // 6 + make_polygon({{0, 90}, {0, 60}, {40, 60}}), // 7 + make_polygon({{-10, 80}, {-10, 50}, {50, 80}}), // 8 + make_polygon({{0, 80}, {0, 50}, {40, 50}, {40, 80}}), // 9 + make_polygon({{0, 65}, {20, 55}, {40, 60}, {20, 65}}), // 10 + make_polygon({{20, 65}, {0, 60}, {20, 55}, {40, 60}}), // 11 + make_polygon({{10, 63}, {20, 55}, {30, 63}, {20, 65}}), // 12 + make_polygon({{20, 75}, {0, 70}, {5, 5}, {10, 0}, {20, 0}, {40, 70}}), // 13 + make_polygon({{0, 50}, {0, 40}, {5, 45}}), // 14 + make_polygon({{0, 90}, {0, 80}, {20, 0}, {40, 80}}), // 15 + make_polygon({{0, 65}, {0, 55}, {40, 65}, {40, 75}}), // 16 }; std::array plg_i = { make_polygon(), //0 @@ -103,154 +201,208 @@ CASE("test_spherical_polygon_intersection") { make_polygon({{0, 50}, {0, 40}, {5, 45}}), make_polygon({{0, 90}, {0, 80}, {20, 0}, {40, 80}}), //32 make_polygon({{0, 65}, {0, 55}, {40, 65}, {40, 75}})}; + Log::info().indent(); for (int i = 0; i < nplg_f; i++) { for (int j = 0; j < nplg_g; j++) { - Log::info() << "\n(" << i * nplg_g + j << ") Intersecting polygon\n " << plg_f[i] << std::endl; - Log::info() << "with polygon\n " << plg_g[j] << std::endl; - auto plg_fg = plg_f[i].intersect(plg_g[j]); - auto plg_gf = plg_g[j].intersect(plg_f[i]); - Log::info() << "got polygon\n "; - if (plg_fg) { - Log::info() << plg_fg << std::endl; - Log::info() << " " << plg_gf << std::endl; - EXPECT(plg_fg.equals(plg_gf)); - EXPECT(plg_fg.equals(plg_i[i * nplg_g + j], 0.1)); - Log::info() << "instead of polygon\n "; - Log::info() << plg_i[i * nplg_g + j] << std::endl; - } - else { - Log::info() << " empty" << std::endl; - Log::info() << "instead of polygon\n "; - Log::info() << plg_i[i * nplg_g + j] << std::endl; + auto plg_fg = plg_f[i].intersect(plg_g[j], nullptr, 5.e6 * std::numeric_limits::epsilon()); + bool polygons_equal = plg_i[i * nplg_g + j].equals(plg_fg, 0.1); + EXPECT(polygons_equal); + if( not polygons_equal or (i==0 && j==10)) { + Log::info() << "\nIntersected the polygon plg_f["<= expected_size); // greater because compiler may add some padding -} +CASE("test_spherical_polygon_intersection_stretched") { + // plg1 is an octant + const auto plg1 = make_polygon({{0, 90}, {0, 0}, {90, 0}}); + Log::info().precision(20); -CASE("analyse intersect") { - const double EPS = std::numeric_limits::epsilon(); - const double du = 0.5; - const double dv = 1.1 * EPS; - const double duc = 0.5 * du; - const double sduc = std::sqrt(1. - 0.25 * du * du); - const double dvc = 1. - 0.5 * dv * dv; - const double sdvc = dv * std::sqrt(1. - 0.25 * dv * dv); - PointXYZ s0p0{sduc, -duc, 0.}; - PointXYZ s0p1{sduc, duc, 0.}; - PointXYZ s1p0{dvc * sduc, -dvc * duc, -sdvc}; - PointXYZ s1p1{dvc * sduc, dvc * duc, sdvc}; + const double delta = 1.1 * EPS; + const double u = 1. - delta; + const double v = std::sqrt(1 - u * u); + const double w = delta; + const double a = sqrt(1 - w * w); - EXPECT_APPROX_EQ(dv, PointXYZ::norm(s0p0 - s1p0), EPS); - EXPECT_APPROX_EQ(du, PointXYZ::norm(s0p0 - s0p1), EPS); - EXPECT_APPROX_EQ(dv, PointXYZ::norm(s0p1 - s1p1), EPS); + // plg2 is a stretched polygon + std::vector plg2_points = { + PointXYZ{u, v, 0.}, + PointXYZ{a, 0., w }, + PointXYZ{u, -v, 0.}, + PointXYZ{0., 0., -1.}}; + auto plg2 = util::ConvexSphericalPolygon(plg2_points.data(), plg2_points.size()); + EXPECT(plg2.size() == 4); - ConvexSphericalPolygon::GreatCircleSegment s1(s0p0, s0p1); - ConvexSphericalPolygon::GreatCircleSegment s2(s1p0, s1p1); + auto plg11 = plg1.intersect(plg1); + auto plg12 = plg1.intersect(plg2); + auto plg22 = plg2.intersect(plg2); - // analytical solution - PointXYZ Isol{1., 0., 0.}; + EXPECT(plg11.size() == 3); + EXPECT(plg12.size() == 3); + EXPECT(plg22.size() == 4); - // test "intersection" - PointXYZ I = s1.intersect(s2); - EXPECT_APPROX_EQ(std::abs(PointXYZ::norm(I) - 1.), 0., EPS); - EXPECT_APPROX_EQ(PointXYZ::norm(I - Isol), 0., EPS); + // plg3 is the analytical solution + std::vector plg3_points = { + PointXYZ{u, v, 0.}, + PointXYZ{a, 0, w}, + PointXYZ{1, 0, 0.}}; + auto plg3 = util::ConvexSphericalPolygon(plg3_points.data(), plg3_points.size()); + EXPECT(plg3.size() == 3); + EXPECT(plg12.size() == 3); + auto plg_sol = make_polygon({{1.2074e-06, 0.}, {0., 1.3994e-14}, {0., 0.}}); + EXPECT(plg12.equals(plg_sol, 1e-10)); - // test "contains" - EXPECT(s1.contains(Isol) && s2.contains(Isol)); - EXPECT(s1.contains(I) && s2.contains(I)); + Log::info().indent(); + Log::info() << "delta : " << delta << std::endl; + Log::info() << "plg12.area : " << plg12.area() << std::endl; + Log::info() << "exact intersection area : " << plg3.area() << std::endl; + double error_area = plg12.area() - plg3.area(); + EXPECT(error_area < EPS and error_area >= 0.); + Log::info().unindent(); + Log::info().precision(-1); } -CASE("source_covered") { - const double dd = 0.; - const auto csp0 = make_polygon({{0, 90}, {0, 0}, {90, 0}}); - double dcov1 = csp0.area(); // optimal coverage - double dcov2 = dcov1; // intersection-based coverage - double dcov3 = dcov1; // normalized intersection-based coverage - double darea = 0; // commutative area error in intersection: |area(A^B)-area(B^A)| - - double max_tarea = 0.; - double min_tarea = 0.; - double accumulated_tarea = 0.; - - const int n = 900; - const int m = 900; - const double dlat = 90. / n; +CASE("test_lonlat_pole_problem") { + const auto north_octant = make_polygon({{0, 90}, {0, 0}, {90, 0}}); + const double first_lat = 90. - 1.e+12 * EPS; + const int m = 10000; const double dlon = 90. / m; + std::vector csp(m); + Log::info().indent(); + + ATLAS_TRACE_SCOPE("create polygons") + for (int j = 0; j < m; j++) { + csp[j] = + make_polygon(PointLonLat{0, 90}, PointLonLat{dlon * j, first_lat}, PointLonLat{dlon * (j + 1), first_lat}); + } + + double max_area_overshoot = 0.; + int false_zero = 0; + for (int j = 0; j < m; j++) { + auto csp_i = csp[j].intersect(north_octant); + // intersection area should not be larger than its father polygon's + max_area_overshoot = std::max(max_area_overshoot, csp_i.area() - csp[j].area()); + if (csp_i.area() < EPS) { + false_zero++; + } + } + Log::info() << "False zero area intersection: " << false_zero << std::endl; + Log::info() << "Max area overshoot: " << max_area_overshoot << std::endl; + EXPECT(max_area_overshoot <= m * EPS); + EXPECT(false_zero == 0); + Log::info().unindent(); +} + +CASE("test_thin_elements_area") { + const auto north_octant = make_polygon({{0, 90}, {0, 0}, {90, 0}}); + const auto south_octant = make_polygon({{0,0}, {0, -90},{90, 0}}); + const int n = 3; + const int m = 2500; + const double dlat = 180. / n; + const double dlon = 90. / m; + + ATLAS_ASSERT(n > 1); std::vector csp(n * m); - std::vector tgt_area(n * m); - ATLAS_TRACE_SCOPE("1") + + ATLAS_TRACE_SCOPE("create polygons") for (int j = 0; j < m; j++) { csp[j] = make_polygon(PointLonLat{0, 90}, PointLonLat{dlon * j, 90 - dlat}, PointLonLat{dlon * (j + 1), 90 - dlat}); } - ATLAS_TRACE_SCOPE("2") - for (int i = 1; i < n; i++) { + for (int i = 1; i < n - 1; i++) { for (int j = 0; j < m; j++) { csp[i * m + j] = make_polygon( PointLonLat{dlon * j, 90 - dlat * i}, PointLonLat{dlon * j, 90 - dlat * (i + 1)}, PointLonLat{dlon * (j + 1), 90 - dlat * (i + 1)}, PointLonLat{dlon * (j + 1), 90 - dlat * i}); } } - ConvexSphericalPolygon cspi0; - ConvexSphericalPolygon csp0i; + for (int j = 0; j < m; j++) { + csp[(n - 1) * m + j] = + make_polygon(PointLonLat{dlon * j, 90 - dlat * (n-1)}, PointLonLat{dlon * j, -90.}, + PointLonLat{dlon * (j + 1), 90 - dlat * (n-1)}); + } + + double coverage_north = 0.; // north octant coverage by intersections with "csp"s + double coverage_south = 0.; // south octant coverage by intersections with "csp"s + double coverage_norm = 0.; // area sum of intersections in the north octant normalised to sum up to the area of the north octant + double coverage_csp = 0.; // area sum of all "csp"s + double accumulated_tarea = 0.; + std::vector i_north_area(n * m); - ATLAS_TRACE_SCOPE("3") + ATLAS_TRACE_SCOPE("intersect polygons") for (int i = 0; i < n; i++) { for (int j = 0; j < m; j++) { auto ipoly = i * m + j; - cspi0 = csp[ipoly].intersect(csp0); // intersect small csp[...] with large polygon csp0 - // should be approx csp[...] as well - csp0i = csp0.intersect(csp[ipoly]); // opposite: intersect csp0 with small csp[...] - // should be approx csp[...] as well - double a_csp = csp[ipoly].area(); - double a_cspi0 = cspi0.area(); // should be approx csp[...].area() - double a_csp0i = csp0i.area(); // should approx match a_cspi0 - EXPECT_APPROX_EQ(a_cspi0, a_csp, 1.e-10); - EXPECT_APPROX_EQ(a_csp0i, a_csp, 1.e-10); - darea = std::max(darea, a_cspi0 - a_csp0i); // should remain approx zero - dcov1 -= csp[ipoly].area(); - dcov2 -= a_cspi0; - tgt_area[ipoly] = a_cspi0; - accumulated_tarea += tgt_area[ipoly]; + double a_north_csp_i = csp[ipoly].intersect(north_octant).area(); // intersect narrow csp with the north octant + double a_south_csp_i = csp[ipoly].intersect(south_octant).area(); // intersect narrow csp with the south octant + if (i == 0) { + if (n == 2) { + EXPECT_APPROX_EQ(a_north_csp_i, csp[ipoly].area(), EPS); + } + else { + if (a_north_csp_i > csp[ipoly].area()) { + Log::info() << " error: " << a_north_csp_i - csp[ipoly].area() << std::endl; + } + } + } + coverage_north += a_north_csp_i; + coverage_csp += csp[ipoly].area(); + coverage_south += a_south_csp_i; + i_north_area[ipoly] = a_north_csp_i; } } - // normalize weights - double norm_fac = csp0.area() / accumulated_tarea; - ATLAS_TRACE_SCOPE("4") + // normalise weights of the intersection polygons to sum up to the area of the north octant + double norm_fac = north_octant.area() / coverage_north; for (int i = 0; i < n; i++) { for (int j = 0; j < m; j++) { - tgt_area[i * m + j] *= norm_fac; - dcov3 -= tgt_area[i * m + j]; + coverage_norm += i_north_area[i * m + j] * norm_fac; } } + EXPECT(north_octant.area() - M_PI_2 < EPS); // spherical area error for the north octant + EXPECT(south_octant.area() - M_PI_2 < EPS); // spherical area error for the south octant + EXPECT(north_octant.area() - coverage_north < m * EPS); // error polygon intersec. north + EXPECT(south_octant.area() - coverage_south < m * EPS); // error polygon intersec. north + EXPECT(north_octant.area() - coverage_norm < m * EPS); // error polygon intersec. north + auto err = std::max(north_octant.area() - coverage_north, south_octant.area() - coverage_south); + auto err_normed = north_octant.area() - coverage_norm; + Log::info() << "Octant coverting error : " << err << std::endl; + Log::info() << "Octant coverting error normed : " << err_normed << std::endl; +} - Log::info() << " dlat, dlon : " << dlat << ", " << dlon << "\n"; - Log::info() << " max (commutative_area) : " << darea << "\n"; - Log::info() << " dcov1 : " << dcov1 << "\n"; - Log::info() << " dcov2 : " << dcov2 << "\n"; - Log::info() << " dcov3 : " << dcov3 << "\n"; - Log::info() << " accumulated small polygon area : " << accumulated_tarea << "\n"; - Log::info() << " large polygon area : " << csp0.area() << "\n"; - EXPECT_APPROX_EQ(darea, 0., 1.e-10); - EXPECT_APPROX_EQ(accumulated_tarea, csp0.area(), 1.e-8); +CASE("intesection of epsilon-distorted polygons") { + const double eps = 1e-4; // degrees + const auto plg0 = make_polygon({{42.45283019, 50.48004426}, + {43.33333333, 49.70239033}, + {44.1509434, 50.48004426}, + {43.26923077, 51.2558069}}); + const auto plg1 = make_polygon({{42.45283019, 50.48004426 - eps}, + {43.33333333 + eps, 49.70239033}, + {44.1509434, 50.48004426 + eps}, + {43.26923077 - eps, 51.2558069}}); + const auto iplg_sol = make_polygon({{44.15088878324276, 50.48009332686897}, + {43.68455392823953, 50.89443301919586}, + {43.26918271448949, 51.25576215711414}, + {42.86876285000331, 50.87921047438197}, + {42.45288468219661, 50.47999711267543}, + {42.92307395320301, 50.06869923211562}, + {43.33338148668725, 49.70243705225555}, + {43.72937844034824, 50.08295071539503}}); + check_intersection(plg0, plg1, iplg_sol); } -CASE("edge cases") { +CASE("Edge cases in polygon intersections") { Log::info().precision(20); - SECTION("CS-LFR-256 -> H1280 problem polygon intersection") { + + SECTION("CS-LFR-256 -> H1280") { const auto plg0 = make_polygon({{-23.55468749999994, -41.11286269132660}, {-23.20312500000000, -41.18816845938357}, {-23.20312500000000, -40.83947225425061}, @@ -260,14 +412,15 @@ CASE("edge cases") { {-23.23828125000000, -40.81704944888558}, {-23.27343750000000, -40.77762996221442}}); auto iplg = plg0.intersect(plg1); - auto jplg = plg1.intersect(plg0); - EXPECT_APPROX_EQ(iplg.area(), jplg.area(), 2e-12); // can not take 1e-15 - EXPECT_EQ(iplg.size(), 3); - EXPECT_EQ(jplg.size(), 3); - EXPECT(iplg.equals(jplg, 5.e-7)); - Log::info() << "Intersection area difference: " << std::abs(iplg.area() - jplg.area()) << "\n"; + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},16) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); } - SECTION("CS-LFR-16 -> O32 problem polygon intersection") { + SECTION("CS-LFR-16 -> O32") { const auto plg0 = make_polygon({{174.3750000000001, -16.79832945594544}, {174.3750000000001, -11.19720014633353}, {168.7500000000000, -11.03919441545243}, @@ -276,25 +429,29 @@ CASE("edge cases") { {174.1935483870968, -15.34836475949100}, {177.0967741935484, -15.34836475949100}}); auto iplg = plg0.intersect(plg1); - auto jplg = plg1.intersect(plg0); - EXPECT_APPROX_EQ(iplg.area(), jplg.area(), 1e-13); // can not take 1e-15 - EXPECT_EQ(iplg.size(), 3); - EXPECT_EQ(jplg.size(), 3); - EXPECT(iplg.equals(jplg, 1.e-11)); - Log::info() << "Intersection area difference: " << std::abs(iplg.area() - jplg.area()) << "\n"; + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},16) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); } - SECTION("CS-LFR-2 -> O48 problem polygon intersection") { + SECTION("CS-LFR-2 -> O48") { const auto plg0 = make_polygon({{180, -45}, {180, 0}, {135, 0}, {135, -35.26438968}}); const auto plg1 = make_polygon({{180, -23.31573073}, {180, -25.18098558}, {-177.75, -23.31573073}}); auto iplg = plg0.intersect(plg1); - auto jplg = plg1.intersect(plg0); - EXPECT_APPROX_EQ(iplg.area(), jplg.area(), 1e-15); EXPECT_EQ(iplg.size(), 0); - EXPECT_EQ(jplg.size(), 0); - EXPECT(iplg.equals(jplg, 1.e-12)); - Log::info() << "Intersection area difference: " << std::abs(iplg.area() - jplg.area()) << "\n"; + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},10) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + EXPECT_APPROX_EQ(iplg.area(), 0.); + Log::info().unindent(); } - SECTION("H128->H256 problem polygon intersection") { + SECTION("H128->H256") { const auto plg0 = make_polygon({{42.45283019, 50.48004426}, {43.33333333, 49.70239033}, {44.15094340, 50.48004426}, @@ -304,12 +461,13 @@ CASE("edge cases") { {44.15094340, 50.48004426}, {43.71428571, 50.86815893}}); auto iplg = plg0.intersect(plg1); - auto jplg = plg1.intersect(plg0); - EXPECT_APPROX_EQ(iplg.area(), jplg.area(), 1e-15); - EXPECT_EQ(iplg.size(), 4); - EXPECT_EQ(jplg.size(), 4); - EXPECT(iplg.equals(jplg, 1.e-12)); - Log::info() << "Intersection area difference: " << std::abs(iplg.area() - jplg.area()) << "\n"; + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},10) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); } SECTION("intesection of epsilon-distorted polygons") { @@ -323,13 +481,104 @@ CASE("edge cases") { {44.1509434, 50.48004426 + eps}, {43.26923077 - eps, 51.2558069}}); auto iplg = plg0.intersect(plg1); - auto jplg = plg1.intersect(plg0); - EXPECT_APPROX_EQ(iplg.area(), jplg.area(), 1e-10); EXPECT_EQ(iplg.size(), 8); - EXPECT_EQ(jplg.size(), 8); - EXPECT(iplg.equals(jplg, 1.e-8)); - Log::info() << "Intersection area difference: " << std::abs(iplg.area() - jplg.area()) << "\n"; + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},10) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); + } + SECTION("one") { // TODO: remove, do not know what example + auto plg0 = make_polygon({{0,90},{67.463999999999998636,89.999899999085130275},{67.5,89.999899999085130275}}); + auto plg1 = make_polygon({{72,85.760587120443801723},{90,85.760587120443801723},{-90,85.760587120443801723}}); + auto iplg = plg0.intersect(plg1); + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},20) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); } + + SECTION("debug") { // TODO: remove, do not know what example + auto plg0 = make_polygon({{168.599999999999994,-6.88524379355500038},{168.561872909699019,-7.16627415158899961},{168.862876254180634,-7.16627415158899961}}); + auto plg1 = make_polygon({{168.84,-6.84},{168.84,-7.2},{169.2,-7.2},{169.2,-6.84}}); + auto plg2 = make_polygon({{168.48,-6.84},{168.48,-7.2},{168.84,-7.2},{168.84,-6.84}}); + auto iplg = plg0.intersect(plg1); + auto iplg2 = plg0.intersect(plg2); + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg2,iplg,iplg2},20) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); + } + + + SECTION("debug2") { // TODO: remove, do not know what example + auto plg0 = make_polygon({{168.599999999999994, -6.88524379355500038}, + {168.561872909699019, -7.16627415158899961}, + {168.862876254180634, -7.16627415158899961}}); + auto plg1 = make_polygon({{168.48, -6.84}, + {168.48, -7.20}, + {168.84, -7.20}, + {168.84, -6.84}}); + const auto iplg_sol = make_polygon({{168.600000000000, -6.885243793555000}, + {168.561872909699, -7.166274151589000}, + {168.840000000000, -7.166281023991750}, + {168.840000000000, -7.141837487873564}}); + check_intersection(plg0, plg1, iplg_sol); + } + + SECTION("O320c_O320n_tcell-2524029") { + auto plg0 = make_polygon({{ 16.497973615369710, 89.85157892074884}, + { 0.000000000000000, 89.87355342974176}, + {-54.000000000000010, 89.78487690721863}, + { -9.000000000000002, 89.84788464490568}}); + auto plg1 = make_polygon({{ 36.000000000000000, 89.78487690721863}, + {-54.000000000000010, 89.78487690721863}, + {-36.000000000000000, 89.78487690721863}}); + const auto iplg_sol = make_polygon(); + check_intersection(plg0, plg1, iplg_sol); + } + + SECTION("O128c_F128c_tcell-77536") { + auto plg0 = make_polygon({{157.5,-16.49119584364},{157.5,-17.192948774143},{158.203125,-17.192948774143},{158.203125,-16.49119584364}}); + auto plg1 = make_polygon({{157.7064220183486,-16.49119584364},{157.5,-17.192948774143},{158.3333333333333,-17.192948774143}}); + const auto iplg_sol = make_polygon({{157.7066386314724, -16.49143953797217}, + {157.7063506671565, -16.49143933951490}, + {157.5000000000000, -17.19294877414300}, + {158.2031250000000, -17.19294877414292}, + {158.2031250000000, -17.04778221576889}}); + check_intersection(plg0, plg1, iplg_sol); + } + + SECTION("O128c_F128c_tcell") { + auto plg0 = make_polygon({{135,-10.505756145244},{135,-11.906523334954},{136.40625,-11.906523334954},{136.40625,-10.505756145244}}); + auto plg1 = make_polygon({{135.7377049180328,-10.505756145244},{135,-11.906523334954},{136.5,-11.906523334954}}); + const auto iplg_sol = make_polygon({{135.7381225488238, -10.50652773728132}, + {135.7373006953013, -10.50652782622945}, + {135.0000000000000, -11.90652333495400}, + {136.4062500000000, -11.90652333495396}, + {136.4062500000000, -11.73509894855489}}); + check_intersection(plg0, plg1, iplg_sol); + } + + SECTION("O128c_F128c_tcell-2") { + auto plg0 = make_polygon({{134.296875,-10.877172064989},{134.296875,-11.578925065131},{135,-11.578925065131},{135,-10.877172064989}}); + auto plg1 = make_polygon({{134.6153846153846,-10.877172064989},{134.2241379310345,-11.578925065131},{135,-11.578925065131}}); + const auto iplg_sol = make_polygon({{134.6154929040914, -10.87737018871372}, + {134.6152744739456, -10.87737016536156}, + {134.2968750000000, -11.44875997313061}, + {134.2968749999999, -11.57892506513095}, + {135.0000000000000, -11.57892506513100}}); + check_intersection(plg0, plg1, iplg_sol); + } + } //----------------------------------------------------------------------------- From 56454e9e542a656f35acc067e7818d61156b600f Mon Sep 17 00:00:00 2001 From: Slavko Brdar Date: Wed, 8 Mar 2023 14:38:32 +0100 Subject: [PATCH 2/3] Fix and diagnose problems with ConservativeSphericalInterpolation --- ...nservativeSphericalPolygonInterpolation.cc | 1083 +++++++++++------ ...onservativeSphericalPolygonInterpolation.h | 42 +- .../atlas-conservative-interpolation.cc | 53 +- .../test_interpolation_conservative.cc | 86 +- 4 files changed, 827 insertions(+), 437 deletions(-) diff --git a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc index b87c17987..863591212 100644 --- a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc +++ b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc @@ -8,8 +8,10 @@ * nor does it submit to any jurisdiction. */ +#include #include #include +#include // for mkdir #include "ConservativeSphericalPolygonInterpolation.h" @@ -22,6 +24,7 @@ #include "atlas/mesh/actions/BuildNode2CellConnectivity.h" #include "atlas/meshgenerator.h" #include "atlas/parallel/mpi/mpi.h" +#include "atlas/parallel/omp/omp.h" #include "atlas/runtime/Exception.h" #include "atlas/runtime/Log.h" #include "atlas/runtime/Trace.h" @@ -31,6 +34,8 @@ #include "eckit/log/Bytes.h" +#define PRINT_BAD_POLYGONS 0 + namespace atlas { namespace interpolation { namespace method { @@ -40,8 +45,99 @@ using util::ConvexSphericalPolygon; namespace { +template +std::string to_json(const It& begin, const It& end, int precision = 0) { + std::stringstream ss; + ss << "[\n"; + size_t size = std::distance(begin,end); + size_t c=0; + for( auto it = begin; it != end; ++it, ++c ) { + ss << " " << it->json(precision); + if( c < size-1 ) { + ss << ",\n"; + } + } + ss << "\n]"; + return ss.str(); +} + +template +std::string to_json(const ConvexSphericalPolygonContainer& polygons, int precision = 0) { + return to_json(polygons.begin(),polygons.end(),precision); +} + + MethodBuilder __builder("conservative-spherical-polygon"); +template +std::string polygons_to_json(const It& begin, const It& end, int precision = 0) { + std::stringstream ss; + ss << "[\n"; + size_t size = std::distance(begin,end); + size_t c=0; + for( auto it = begin; it != end; ++it, ++c ) { + ss << " " << it->json(precision); + if( c < size-1 ) { + ss << ",\n"; + } + } + ss << "\n]"; + return ss.str(); +} + +template +std::string polygons_to_json(const ConvexSphericalPolygonContainer& polygons, int precision = 0) { + return polygons_to_json(polygons.begin(),polygons.end(),precision); +} + +void dump_polygons_to_json( const ConvexSphericalPolygon& t_csp, + double pointsSameEPS, + const std::vector>& source_polygons, + const std::vector& source_polygons_considered_indices, + const std::string folder, + const std::string name) { + std::vector csp_arr{ t_csp }; + std::vector csp_arr_intersecting {t_csp}; + std::vector intersections; + int count = 1; + for( auto& s_idx : source_polygons_considered_indices ) { + auto s_csp = std::get<0>(source_polygons[s_idx]); + csp_arr.emplace_back( s_csp ); + std::fstream file_plg_debug(folder + name + "_" + std::to_string(count++) + ".debug", std::ios::out); + ConvexSphericalPolygon iplg = t_csp.intersect(s_csp, &file_plg_debug, pointsSameEPS); + file_plg_debug.close(); + if (iplg) { + if( iplg.area() > 0. ) { + csp_arr_intersecting.emplace_back( s_csp ); + intersections.emplace_back( iplg ); + } + } + } + double tot = 0.; + Log::info().indent(); + std::fstream file_info(folder + name + ".info", std::ios::out); + file_info << "List of intersection weights: " << std::endl; + count = 1; + for( auto& iplg : intersections ) { + csp_arr.emplace_back( iplg ); + csp_arr_intersecting.emplace_back( iplg ); + tot += iplg.area() / t_csp.area(); + file_info << "\t" << count++ << ": " << iplg.area() / t_csp.area() << std::endl; + } + file_info << std::endl << name + ": " << 100. * tot << " % covered."<< std::endl << std::endl; + file_info << "Target polygon + candidate source polygons + " << intersections.size() << " intersections in file:" << std::endl; + file_info << "\t" << folder + name + ".candidates\n" << std::endl; + std::fstream file_plg(folder + name + ".candidates", std::ios::out); + file_plg << polygons_to_json(csp_arr, 16); + file_plg.close(); + file_info << "Target polygon + intersecting source polygon + " << intersections.size() << " intersections in file:" << std::endl; + file_info << "\t" << folder + name + ".intersections\n" << std::endl; + file_plg.open(folder + name + ".intersections", std::ios::out); + file_plg << polygons_to_json(csp_arr_intersecting, 16); + file_plg.close(); + Log::info().unindent(); +} + constexpr double unit_sphere_area() { // 4*pi*r^2 with r=1 return 4. * M_PI; @@ -49,7 +145,7 @@ constexpr double unit_sphere_area() { template size_t memory_of(const std::vector& vector) { - return sizeof(T) * vector.size(); + return sizeof(T) * vector.capacity(); } template size_t memory_of(const std::vector>& vector_of_vector) { @@ -66,7 +162,7 @@ size_t memory_of( for (const auto& params : vector_of_params) { mem += memory_of(params.cell_idx); mem += memory_of(params.centroids); - mem += memory_of(params.src_weights); + mem += memory_of(params.weights); mem += memory_of(params.tgt_weights); } return mem; @@ -108,9 +204,6 @@ void sort_and_accumulate_triplets(std::vector& triplets) } } - -} // namespace - int inside_vertices(const ConvexSphericalPolygon& plg1, const ConvexSphericalPolygon& plg2, int& pout) { int points_in = 0; pout = 0; @@ -132,6 +225,8 @@ int inside_vertices(const ConvexSphericalPolygon& plg1, const ConvexSphericalPol return points_in; } +} // namespace + ConservativeSphericalPolygonInterpolation::ConservativeSphericalPolygonInterpolation(const Config& config): Method(config) { config.get("validate", validate_ = false); @@ -321,8 +416,9 @@ std::vector ConservativeSphericalPolygonInterpolation::get_node_neighbour // Create polygons for cell-centred data. Here, the polygons are mesh cells ConservativeSphericalPolygonInterpolation::CSPolygonArray -ConservativeSphericalPolygonInterpolation::get_polygons_celldata(Mesh& mesh) const { +ConservativeSphericalPolygonInterpolation::get_polygons_celldata(FunctionSpace fs) const { CSPolygonArray cspolygons; + auto mesh = extract_mesh(fs); const idx_t n_cells = mesh.cells().size(); cspolygons.resize(n_cells); const auto& cell2node = mesh.cells().node_connectivity(); @@ -331,7 +427,11 @@ ConservativeSphericalPolygonInterpolation::get_polygons_celldata(Mesh& mesh) con const auto& cell_flags = array::make_view(mesh.cells().flags()); const auto& cell_part = array::make_view(mesh.cells().partition()); std::vector pts_ll; + const int fs_halo = functionspace::CellColumns(fs).halo().size(); for (idx_t cell = 0; cell < n_cells; ++cell) { + if( cell_halo(cell) > fs_halo ) { + continue; + } int halo_type = cell_halo(cell); const idx_t n_nodes = cell2node.cols(cell); pts_ll.clear(); @@ -354,12 +454,13 @@ ConservativeSphericalPolygonInterpolation::get_polygons_celldata(Mesh& mesh) con // (cell_centre, edge_centre, cell_vertex, edge_centre) // additionally, subcell-to-node and node-to-subcells mapping are computed ConservativeSphericalPolygonInterpolation::CSPolygonArray -ConservativeSphericalPolygonInterpolation::get_polygons_nodedata(Mesh& mesh, std::vector& csp2node, +ConservativeSphericalPolygonInterpolation::get_polygons_nodedata(FunctionSpace fs, std::vector& csp2node, std::vector>& node2csp, std::array& errors) const { CSPolygonArray cspolygons; csp2node.clear(); node2csp.clear(); + auto mesh = extract_mesh(fs); node2csp.resize(mesh.nodes().size()); const auto nodes_ll = array::make_view(mesh.nodes().lonlat()); const auto& cell2node = mesh.cells().node_connectivity(); @@ -376,9 +477,19 @@ ConservativeSphericalPolygonInterpolation::get_polygons_nodedata(Mesh& mesh, std eckit::geometry::Sphere::convertSphericalToCartesian(1., p_ll, p_xyz); return p_xyz; }; + const auto node_halo = array::make_view(mesh.nodes().halo()); + const auto node_ghost = array::make_view(mesh.nodes().ghost()); + const auto node_flags = array::make_view(mesh.nodes().flags()); + const auto node_part = array::make_view(mesh.nodes().partition()); + idx_t cspol_id = 0; // subpolygon enumeration errors = {0., 0.}; // over/undershoots in creation of subpolygons + const int fs_halo = functionspace::NodeColumns(fs).halo().size(); + for (idx_t cell = 0; cell < mesh.cells().size(); ++cell) { + if( cell_halo(cell) > fs_halo ) { + continue; + } ATLAS_ASSERT(cell < cell2node.rows()); const idx_t n_nodes = cell2node.cols(cell); ATLAS_ASSERT(n_nodes > 2); @@ -409,6 +520,7 @@ ConservativeSphericalPolygonInterpolation::get_polygons_nodedata(Mesh& mesh, std PointLonLat cell_ll = xyz2ll(cell_mid); double loc_csp_area_shoot = ConvexSphericalPolygon(pts_ll).area(); // get ConvexSphericalPolygon for each valid edge + int halo_type; for (int inode = 0; inode < pts_idx.size(); inode++) { int inode_n = next_index(inode, pts_idx.size()); idx_t node = cell2node(cell, inode); @@ -429,18 +541,22 @@ ConservativeSphericalPolygonInterpolation::get_polygons_nodedata(Mesh& mesh, std subpol_pts_ll[1] = xyz2ll(iedge_mid); subpol_pts_ll[2] = pts_ll[inode_n]; subpol_pts_ll[3] = xyz2ll(jedge_mid); - int halo_type = cell_halo(cell); - if (util::Bitflags::view(cell_flags(cell)).check(util::Topology::PERIODIC) and - cell_part(cell) == mpi::rank()) { + halo_type = node_halo(node_n); + + if (util::Bitflags::view(node_flags(node_n)).check(util::Topology::PERIODIC) and + node_part(node_n) == mpi::rank()) { halo_type = -1; } + ConvexSphericalPolygon cspi(subpol_pts_ll); loc_csp_area_shoot -= cspi.area(); cspolygons.emplace_back(cspi, halo_type); cspol_id++; } - errors[0] += std::abs(loc_csp_area_shoot); - errors[1] = std::max(std::abs(loc_csp_area_shoot), errors[1]); + if (halo_type == 0) { + errors[0] += std::abs(loc_csp_area_shoot); + errors[1] = std::max(std::abs(loc_csp_area_shoot), errors[1]); + } } ATLAS_TRACE_MPI(ALLREDUCE) { mpi::comm().allReduceInPlace(&errors[0], 1, eckit::mpi::sum()); @@ -465,7 +581,7 @@ void ConservativeSphericalPolygonInterpolation::do_setup_impl(const Grid& src_gr tgt_fs_ = functionspace::CellColumns(tgt_mesh_, option::halo(0)); } else { - tgt_fs_ = functionspace::NodeColumns(tgt_mesh_, option::halo(0)); + tgt_fs_ = functionspace::NodeColumns(tgt_mesh_, option::halo(1)); } } sharable_data_->tgt_fs_ = tgt_fs_; @@ -561,48 +677,68 @@ void ConservativeSphericalPolygonInterpolation::do_setup(const FunctionSpace& sr tgt_mesh_ = extract_mesh(tgt_fs_); { - // we need src_halo_size >= 2, whereas tgt_halo_size >= 0 is enough - int src_halo_size = 0; - src_mesh_.metadata().get("halo", src_halo_size); - ATLAS_ASSERT(src_halo_size > 1); + // we need src_halo_size >= 2, tgt_halo_size >= 0 for CellColumns + // if target is NodeColumns, we need: + // tgt_halo_size >= 1 and + // src_halo_size large enough to cover the the target halo cells in the first row + int halo_size = 0; + src_mesh_.metadata().get("halo", halo_size); + if (halo_size < 2) { + Log::info() << "WARNING The halo size on source mesh should be at least 2.\n"; + } + if (not tgt_cell_data_) { + Log::info() << "WARNING The source cells should cover the first row of the target halos.\n"; + } + halo_size = 0; + tgt_mesh_.metadata().get("halo", halo_size); + if (not tgt_cell_data_ and halo_size == 0) { + Log::info() << "WARNING The halo size on target mesh should be at least 1 for the target NodeColumns.\n"; + } } CSPolygonArray src_csp; CSPolygonArray tgt_csp; - std::array errors = {0., 0.}; if (compute_cache) { + std::array errors = {0., 0.}; ATLAS_TRACE("Get source polygons"); StopWatch stopwatch; stopwatch.start(); if (src_cell_data_) { - src_csp = get_polygons_celldata(src_mesh_); + src_csp = get_polygons_celldata(src_fs_); } else { src_csp = - get_polygons_nodedata(src_mesh_, sharable_data_->src_csp2node_, sharable_data_->src_node2csp_, errors); + get_polygons_nodedata(src_fs_, sharable_data_->src_csp2node_, sharable_data_->src_node2csp_, errors); } stopwatch.stop(); sharable_data_->timings.source_polygons_assembly = stopwatch.elapsed(); - } - remap_stat_.errors[Statistics::Errors::SRC_PLG_L1] = errors[0]; - remap_stat_.errors[Statistics::Errors::SRC_PLG_LINF] = errors[1]; - if (compute_cache) { + remap_stat_.counts[Statistics::Counts::SRC_PLG] = src_csp.size(); + remap_stat_.errors[Statistics::Errors::SRC_SUBPLG_L1] = errors[0]; + remap_stat_.errors[Statistics::Errors::SRC_SUBPLG_LINF] = errors[1]; + + errors = {0., 0.}; ATLAS_TRACE("Get target polygons"); - StopWatch stopwatch; stopwatch.start(); if (tgt_cell_data_) { - tgt_csp = get_polygons_celldata(tgt_mesh_); + tgt_csp = get_polygons_celldata(tgt_fs_); } else { tgt_csp = - get_polygons_nodedata(tgt_mesh_, sharable_data_->tgt_csp2node_, sharable_data_->tgt_node2csp_, errors); + get_polygons_nodedata(tgt_fs_, sharable_data_->tgt_csp2node_, sharable_data_->tgt_node2csp_, errors); } stopwatch.stop(); sharable_data_->timings.target_polygons_assembly = stopwatch.elapsed(); + remap_stat_.counts[Statistics::Counts::TGT_PLG] = tgt_csp.size(); + remap_stat_.errors[Statistics::Errors::TGT_SUBPLG_L1] = errors[0]; + remap_stat_.errors[Statistics::Errors::TGT_SUBPLG_LINF] = errors[1]; + } + else { + remap_stat_.counts[Statistics::Counts::SRC_PLG] = -1111; + remap_stat_.errors[Statistics::Errors::SRC_SUBPLG_L1] = -1111.; + remap_stat_.errors[Statistics::Errors::SRC_SUBPLG_LINF] = -1111.; + remap_stat_.counts[Statistics::Counts::TGT_PLG] = -1111; + remap_stat_.errors[Statistics::Errors::TGT_SUBPLG_L1] = -1111.; + remap_stat_.errors[Statistics::Errors::TGT_SUBPLG_LINF] = -1111.; } - remap_stat_.counts[Statistics::Counts::SRC_PLG] = src_csp.size(); - remap_stat_.counts[Statistics::Counts::TGT_PLG] = tgt_csp.size(); - remap_stat_.errors[Statistics::Errors::TGT_PLG_L1] = errors[0]; - remap_stat_.errors[Statistics::Errors::TGT_PLG_LINF] = errors[1]; n_spoints_ = src_fs_.size(); n_tpoints_ = tgt_fs_.size(); @@ -611,82 +747,124 @@ void ConservativeSphericalPolygonInterpolation::do_setup(const FunctionSpace& sr intersect_polygons(src_csp, tgt_csp); auto& src_points_ = sharable_data_->src_points_; - auto& tgt_points_ = sharable_data_->tgt_points_; src_points_.resize(n_spoints_); - tgt_points_.resize(n_tpoints_); sharable_data_->src_areas_.resize(n_spoints_); auto& src_areas_v = sharable_data_->src_areas_; - if (src_cell_data_) { - for (idx_t spt = 0; spt < n_spoints_; ++spt) { - const auto& s_csp = std::get<0>(src_csp[spt]); - src_points_[spt] = s_csp.centroid(); - src_areas_v[spt] = s_csp.area(); - } - } - else { - auto& src_node2csp_ = sharable_data_->src_node2csp_; - const auto lonlat = array::make_view(src_mesh_.nodes().lonlat()); - for (idx_t spt = 0; spt < n_spoints_; ++spt) { - if (src_node2csp_[spt].size() == 0) { - // this is a node to which no subpolygon is associated - // maximal twice per mesh we end here, and that is only when mesh has nodes on poles - auto p = PointLonLat{lonlat(spt, 0), lonlat(spt, 1)}; - eckit::geometry::Sphere::convertSphericalToCartesian(1., p, src_points_[spt]); - } - else { - // .. in the other case, start computing the barycentre - src_points_[spt] = PointXYZ{0., 0., 0.}; + { + ATLAS_TRACE("Store src_areas and src_point"); + if (src_cell_data_) { + for (idx_t spt = 0; spt < n_spoints_; ++spt) { + const auto& s_csp = std::get<0>(src_csp[spt]); + src_points_[spt] = s_csp.centroid(); + src_areas_v[spt] = s_csp.area(); } - src_areas_v[spt] = 0.; - for (idx_t isubcell = 0; isubcell < src_node2csp_[spt].size(); ++isubcell) { - idx_t subcell = src_node2csp_[spt][isubcell]; - const auto& s_csp = std::get<0>(src_csp[subcell]); - src_areas_v[spt] += s_csp.area(); - src_points_[spt] = src_points_[spt] + PointXYZ::mul(s_csp.centroid(), s_csp.area()); + } + else { + auto& src_node2csp_ = sharable_data_->src_node2csp_; + const auto lonlat = array::make_view(src_mesh_.nodes().lonlat()); + for (idx_t spt = 0; spt < n_spoints_; ++spt) { + if (src_node2csp_[spt].size() == 0) { + // this is a node to which no subpolygon is associated + // maximal twice per mesh we end here, and that is only when mesh has nodes on poles + auto p = PointLonLat{lonlat(spt, 0), lonlat(spt, 1)}; + eckit::geometry::Sphere::convertSphericalToCartesian(1., p, src_points_[spt]); + } + else { + // .. in the other case, start computing the barycentre + src_points_[spt] = PointXYZ{0., 0., 0.}; + } + src_areas_v[spt] = 0.; + for (idx_t isubcell = 0; isubcell < src_node2csp_[spt].size(); ++isubcell) { + idx_t subcell = src_node2csp_[spt][isubcell]; + const auto& s_csp = std::get<0>(src_csp[subcell]); + src_areas_v[spt] += s_csp.area(); + src_points_[spt] = src_points_[spt] + PointXYZ::mul(s_csp.centroid(), s_csp.area()); + } + double src_point_norm = PointXYZ::norm(src_points_[spt]); + if (src_point_norm == 0.) { + ATLAS_DEBUG_VAR(src_point_norm); + ATLAS_DEBUG_VAR(src_points_[spt]); + ATLAS_DEBUG_VAR(src_node2csp_[spt].size()); + for (idx_t isubcell = 0; isubcell < src_node2csp_[spt].size(); ++isubcell) { + idx_t subcell = src_node2csp_[spt][isubcell]; + ATLAS_DEBUG_VAR(subcell); + const auto& s_csp = std::get<0>(src_csp[subcell]); + s_csp.print(Log::info()); + Log::info() << std::endl; + src_areas_v[spt] += s_csp.area(); + ATLAS_DEBUG_VAR(s_csp.area()); + ATLAS_DEBUG_VAR(s_csp.centroid()); + src_points_[spt] = src_points_[spt] + PointXYZ::mul(s_csp.centroid(), s_csp.area()); + } + Log::info().flush(); + // something went wrong, improvise + src_point_norm = 1.; + } + src_points_[spt] = PointXYZ::div(src_points_[spt], src_point_norm); } - double src_point_norm = PointXYZ::norm(src_points_[spt]); - ATLAS_ASSERT(src_point_norm > 0.); - src_points_[spt] = PointXYZ::div(src_points_[spt], src_point_norm); } } + auto& tgt_points_ = sharable_data_->tgt_points_; + tgt_points_.resize(n_tpoints_); sharable_data_->tgt_areas_.resize(n_tpoints_); auto& tgt_areas_v = sharable_data_->tgt_areas_; - if (tgt_cell_data_) { - for (idx_t tpt = 0; tpt < n_tpoints_; ++tpt) { - const auto& t_csp = std::get<0>(tgt_csp[tpt]); - tgt_points_[tpt] = t_csp.centroid(); - tgt_areas_v[tpt] = t_csp.area(); - } - } - else { - auto& tgt_node2csp_ = sharable_data_->tgt_node2csp_; - const auto lonlat = array::make_view(tgt_mesh_.nodes().lonlat()); - for (idx_t tpt = 0; tpt < n_tpoints_; ++tpt) { - if (tgt_node2csp_[tpt].size() == 0) { - // this is a node to which no subpolygon is associated - // maximal twice per mesh we end here, and that is only when mesh has nodes on poles - auto p = PointLonLat{lonlat(tpt, 0), lonlat(tpt, 1)}; - eckit::geometry::Sphere::convertSphericalToCartesian(1., p, tgt_points_[tpt]); - } - else { - // .. in the other case, start computing the barycentre - tgt_points_[tpt] = PointXYZ{0., 0., 0.}; + { + ATLAS_TRACE("Store src_areas and src_point"); + if (tgt_cell_data_) { + for (idx_t tpt = 0; tpt < n_tpoints_; ++tpt) { + const auto& t_csp = std::get<0>(tgt_csp[tpt]); + tgt_points_[tpt] = t_csp.centroid(); + tgt_areas_v[tpt] = t_csp.area(); } - tgt_areas_v[tpt] = 0.; - for (idx_t isubcell = 0; isubcell < tgt_node2csp_[tpt].size(); ++isubcell) { - idx_t subcell = tgt_node2csp_[tpt][isubcell]; - const auto& t_csp = std::get<0>(tgt_csp[subcell]); - tgt_areas_v[tpt] += t_csp.area(); - tgt_points_[tpt] = tgt_points_[tpt] + PointXYZ::mul(t_csp.centroid(), t_csp.area()); + } + else { + auto& tgt_node2csp_ = sharable_data_->tgt_node2csp_; + const auto lonlat = array::make_view(tgt_mesh_.nodes().lonlat()); + for (idx_t tpt = 0; tpt < n_tpoints_; ++tpt) { + if (tgt_node2csp_[tpt].size() == 0) { + // this is a node to which no subpolygon is associated + // maximal twice per mesh we end here, and that is only when mesh has nodes on poles + auto p = PointLonLat{lonlat(tpt, 0), lonlat(tpt, 1)}; + eckit::geometry::Sphere::convertSphericalToCartesian(1., p, tgt_points_[tpt]); + } + else { + // .. in the other case, start computing the barycentre + tgt_points_[tpt] = PointXYZ{0., 0., 0.}; + } + tgt_areas_v[tpt] = 0.; + for (idx_t isubcell = 0; isubcell < tgt_node2csp_[tpt].size(); ++isubcell) { + idx_t subcell = tgt_node2csp_[tpt][isubcell]; + const auto& t_csp = std::get<0>(tgt_csp[subcell]); + tgt_areas_v[tpt] += t_csp.area(); + tgt_points_[tpt] = tgt_points_[tpt] + PointXYZ::mul(t_csp.centroid(), t_csp.area()); + } + double tgt_point_norm = PointXYZ::norm(tgt_points_[tpt]); + ATLAS_ASSERT(tgt_point_norm > 0.); + if (tgt_point_norm == 0.) { + ATLAS_DEBUG_VAR(tgt_point_norm); + ATLAS_DEBUG_VAR(tgt_points_[tpt]); + ATLAS_DEBUG_VAR(tgt_node2csp_[tpt].size()); + for (idx_t isubcell = 0; isubcell < tgt_node2csp_[tpt].size(); ++isubcell) { + idx_t subcell = tgt_node2csp_[tpt][isubcell]; + ATLAS_DEBUG_VAR(subcell); + const auto& t_csp = std::get<0>(tgt_csp[subcell]); + t_csp.print(Log::info()); + Log::info() << std::endl; + tgt_areas_v[tpt] += t_csp.area(); + ATLAS_DEBUG_VAR(t_csp.area()); + ATLAS_DEBUG_VAR(t_csp.centroid()); + tgt_points_[tpt] = tgt_points_[tpt] + PointXYZ::mul(t_csp.centroid(), t_csp.area()); + } + Log::info().flush(); + // something went wrong, improvise + tgt_point_norm = 1.; + } + tgt_points_[tpt] = PointXYZ::div(tgt_points_[tpt], tgt_point_norm); } - double tgt_point_norm = PointXYZ::norm(tgt_points_[tpt]); - ATLAS_ASSERT(tgt_point_norm > 0.); - tgt_points_[tpt] = PointXYZ::div(tgt_points_[tpt], tgt_point_norm); } } } - if (not matrix_free_) { StopWatch stopwatch; stopwatch.start(); @@ -751,8 +929,10 @@ void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolyg util::KDTree kdt_search; kdt_search.reserve(tgt_csp.size()); double max_tgtcell_rad = 0.; + const int tgt_halo_intersection_depth = (tgt_cell_data_ ? 0 : 1); // if target NodeColumns, one target halo required for subcells around target nodes for (idx_t jcell = 0; jcell < tgt_csp.size(); ++jcell) { - if (std::get<1>(tgt_csp[jcell]) == 0) { + auto tgt_halo_type = std::get<1>(tgt_csp[jcell]); + if (tgt_halo_type <= tgt_halo_intersection_depth) { // and tgt_halo_type != -1) { const auto& t_csp = std::get<0>(tgt_csp[jcell]); kdt_search.insert(t_csp.centroid(), jcell); max_tgtcell_rad = std::max(max_tgtcell_rad, t_csp.radius()); @@ -768,17 +948,11 @@ void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolyg stopwatch_src_already_in.start(); std::set src_cent; - auto polygon_point = [](const ConvexSphericalPolygon& pol) { - PointXYZ p{0., 0., 0.}; - for (int i = 0; i < pol.size(); i++) { - p = p + pol[i]; - } - p /= pol.size(); - return p; - }; auto src_already_in = [&](const PointXYZ& halo_cent) { if (src_cent.find(halo_cent) == src_cent.end()) { - src_cent.insert(halo_cent); + atlas_omp_critical{ + src_cent.insert(halo_cent); + } return false; } return true; @@ -807,79 +981,114 @@ void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolyg tgt_iparam.resize(tgt_csp.size()); } + // the worst target polygon coverage for analysis of intersection + std::pair worst_tgt_overcover; + std::pair worst_tgt_undercover; + worst_tgt_overcover.first = -1; + worst_tgt_overcover.second = -1.; + worst_tgt_undercover.first = -1; + worst_tgt_undercover.second = M_PI; + + // NOTE: polygon vertex points at distance < pointsSameEPS will be replaced with one point + constexpr double pointsSameEPS = 5.e6 * std::numeric_limits::epsilon(); + eckit::Channel blackhole; - eckit::ProgressTimer progress("Intersecting polygons ", src_csp.size(), " cell", double(10), - src_csp.size() > 50 ? Log::info() : blackhole); - for (idx_t scell = 0; scell < src_csp.size(); ++scell, ++progress) { - stopwatch_src_already_in.start(); - if (src_already_in(polygon_point(std::get<0>(src_csp[scell])))) { + eckit::ProgressTimer progress("Intersecting polygons ", src_csp.size() / atlas_omp_get_max_threads(), " (cell/thread)", double(10), + src_csp.size() / atlas_omp_get_max_threads() > 50 ? Log::info() : blackhole); + atlas_omp_parallel_for (idx_t scell = 0; scell < src_csp.size(); ++scell) { + if ( atlas_omp_get_thread_num() == 0 ) { + ++progress; + stopwatch_src_already_in.start(); + } + bool already_in = src_already_in((std::get<0>(src_csp[scell]).centroid())); + if ( atlas_omp_get_thread_num() == 0 ) { stopwatch_src_already_in.stop(); - continue; } - stopwatch_src_already_in.stop(); - - const auto& s_csp = std::get<0>(src_csp[scell]); - const double s_csp_area = s_csp.area(); - double src_cover_area = 0.; - - stopwatch_kdtree_search.start(); - auto tgt_cells = kdt_search.closestPointsWithinRadius(s_csp.centroid(), s_csp.radius() + max_tgtcell_rad); - stopwatch_kdtree_search.stop(); - for (idx_t ttcell = 0; ttcell < tgt_cells.size(); ++ttcell) { - auto tcell = tgt_cells[ttcell].payload(); - const auto& t_csp = std::get<0>(tgt_csp[tcell]); - stopwatch_polygon_intersections.start(); - ConvexSphericalPolygon csp_i = s_csp.intersect(t_csp); - double csp_i_area = csp_i.area(); - stopwatch_polygon_intersections.stop(); - if (validate_) { - // check zero area intersections with inside_vertices - int pout; - if (inside_vertices(s_csp, t_csp, pout) > 2 && csp_i.area() < 3e-16) { - dump_intersection(s_csp, tgt_csp, tgt_cells); + if (not already_in) { + const auto& s_csp = std::get<0>(src_csp[scell]); + const double s_csp_area = s_csp.area(); + double src_cover_area = 0.; + + stopwatch_kdtree_search.start(); + auto tgt_cells = kdt_search.closestPointsWithinRadius(s_csp.centroid(), s_csp.radius() + max_tgtcell_rad); + stopwatch_kdtree_search.stop(); + for (idx_t ttcell = 0; ttcell < tgt_cells.size(); ++ttcell) { + auto tcell = tgt_cells[ttcell].payload(); + const auto& t_csp = std::get<0>(tgt_csp[tcell]); + if( atlas_omp_get_thread_num() == 0 ) { + stopwatch_polygon_intersections.start(); + } + ConvexSphericalPolygon csp_i = s_csp.intersect(t_csp, nullptr, pointsSameEPS); + double csp_i_area = csp_i.area(); + if( atlas_omp_get_thread_num() == 0 ) { + stopwatch_polygon_intersections.stop(); } - } - if (csp_i_area > 0.) { if (validate_) { - tgt_iparam[tcell].cell_idx.emplace_back(scell); - tgt_iparam[tcell].tgt_weights.emplace_back(csp_i_area); + int pout; + // TODO: this can be removed soon + if (inside_vertices(s_csp, t_csp, pout) > 2 && csp_i.area() < 3e-16) { + dump_intersection("Zero area intersections with inside_vertices", s_csp, tgt_csp, tgt_cells); + } + // TODO: assuming intersector search works fine, this should be move under "if (csp_i_area > 0)" + atlas_omp_critical { + tgt_iparam[tcell].cell_idx.emplace_back(scell); + tgt_iparam[tcell].tgt_weights.emplace_back(csp_i_area); + } } - src_iparam_[scell].cell_idx.emplace_back(tcell); - src_iparam_[scell].src_weights.emplace_back(csp_i_area); - double target_weight = csp_i_area / t_csp.area(); - src_iparam_[scell].tgt_weights.emplace_back(target_weight); - src_iparam_[scell].centroids.emplace_back(csp_i.centroid()); - src_cover_area += csp_i_area; - ATLAS_ASSERT(target_weight < 1.1); - ATLAS_ASSERT(csp_i_area / s_csp_area < 1.1); - } - } - const double src_cover_err = std::abs(s_csp_area - src_cover_area); - const double src_cover_err_percent = 100. * src_cover_err / s_csp_area; - if (src_cover_err_percent > 0.1 and std::get<1>(src_csp[scell]) == 0) { - // HACK: source cell at process boundary will not be covered by target cells, skip them - // TODO: mark these source cells beforehand and compute error in them among the processes - - if (validate_) { - if (mpi::size() == 1) { - Log::info() << "WARNING src cell covering error : " << src_cover_err_percent << "%\n"; - dump_intersection(s_csp, tgt_csp, tgt_cells); + if (csp_i_area > 0) { + src_iparam_[scell].cell_idx.emplace_back(tcell); + src_iparam_[scell].weights.emplace_back(csp_i_area); + double target_weight = csp_i_area / t_csp.area(); + src_iparam_[scell].tgt_weights.emplace_back(target_weight); // TODO: tgt_weights vector should be removed for the sake of highres + if (order_ == 2 or not matrix_free_ or not matrixAllocated()) { + src_iparam_[scell].centroids.emplace_back(csp_i.centroid()); + } + src_cover_area += csp_i_area; + + if (validate_) { + // this check is concerned with accuracy of polygon intersections + if (target_weight > 1.1) { + dump_intersection("Intersection larger than target", s_csp, tgt_csp, tgt_cells); + } + if (csp_i_area / s_csp_area > 1.1) { + dump_intersection("Intersection larger than source", s_csp, tgt_csp, tgt_cells); + } + } } } - area_coverage[TOTAL_SRC] += src_cover_err; - area_coverage[MAX_SRC] = std::max(area_coverage[MAX_SRC], src_cover_err); - } - if (src_iparam_[scell].cell_idx.size() == 0) { - num_pol[SRC_NONINTERSECT]++; - } - if (normalise_intersections_ && src_cover_err_percent < 1.) { - double wfactor = s_csp.area() / (src_cover_area > 0. ? src_cover_area : 1.); - for (idx_t i = 0; i < src_iparam_[scell].src_weights.size(); i++) { - src_iparam_[scell].src_weights[i] *= wfactor; - src_iparam_[scell].tgt_weights[i] *= wfactor; + const double src_cover_err = std::abs(s_csp_area - src_cover_area); + const double src_cover_err_percent = 100. * src_cover_err / s_csp_area; + if (src_cover_err_percent > 0.1 and std::get<1>(src_csp[scell]) == 0) { + // NOTE: source cell at process boundary will not be covered by target cells, skip them + // TODO: mark these source cells beforehand and compute error in them among the processes + if (validate_ and mpi::size() == 1) { + dump_intersection("Source cell not exactly covered", s_csp, tgt_csp, tgt_cells); + if (statistics_intersection_) { + atlas_omp_critical{ + area_coverage[TOTAL_SRC] += src_cover_err; + area_coverage[MAX_SRC] = std::max(area_coverage[MAX_SRC], src_cover_err); + } + } + } } - } - num_pol[SRC_TGT_INTERSECT] += src_iparam_[scell].src_weights.size(); + if (src_iparam_[scell].cell_idx.size() == 0 and statistics_intersection_) { + atlas_omp_critical{ + num_pol[SRC_NONINTERSECT]++; + } + } + if (normalise_intersections_ && src_cover_err_percent < 1.) { + double wfactor = s_csp.area() / (src_cover_area > 0. ? src_cover_area : 1.); + for (idx_t i = 0; i < src_iparam_[scell].weights.size(); i++) { + src_iparam_[scell].weights[i] *= wfactor; + src_iparam_[scell].tgt_weights[i] *= wfactor; + } + } + if (statistics_intersection_) { + atlas_omp_critical{ + num_pol[SRC_TGT_INTERSECT] += src_iparam_[scell].weights.size(); + } + } + } // already in } timings.polygon_intersections = stopwatch_polygon_intersections.elapsed(); timings.target_kdtree_search = stopwatch_kdtree_search.elapsed(); @@ -888,48 +1097,151 @@ void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolyg num_pol[TGT] = tgt_csp.size(); ATLAS_TRACE_MPI(ALLREDUCE) { mpi::comm().allReduceInPlace(num_pol.data(), num_pol.size(), eckit::mpi::sum()); - mpi::comm().allReduceInPlace(area_coverage.data(), area_coverage.size(), eckit::mpi::max()); } remap_stat_.counts[Statistics::Counts::INT_PLG] = num_pol[SRC_TGT_INTERSECT]; remap_stat_.counts[Statistics::Counts::UNCVR_SRC] = num_pol[SRC_NONINTERSECT]; - remap_stat_.errors[Statistics::Errors::GEO_L1] = area_coverage[TOTAL_SRC]; - remap_stat_.errors[Statistics::Errors::GEO_LINF] = area_coverage[MAX_SRC]; - - double geo_err_l1 = 0.; - double geo_err_linf = 0.; - for (idx_t scell = 0; scell < src_csp.size(); ++scell) { - const int cell_flag = std::get<1>(src_csp[scell]); - if (cell_flag == -1 or cell_flag > 0) { - // skip periodic & halo cells - continue; + + const std::string polygon_intersection_folder = "polygon_intersection/"; + if (validate_ && mpi::rank() == 0) { + if (mkdir(polygon_intersection_folder.c_str(), 0777) != 0) { + Log::info() << "WARNING Polygon intersection relevant information in is the folder \e[1mpolygon_intersection\e[0m." << std::endl; } - double diff_cell = std::get<0>(src_csp[scell]).area(); - for (idx_t icell = 0; icell < src_iparam_[scell].src_weights.size(); ++icell) { - diff_cell -= src_iparam_[scell].src_weights[icell]; + else { + Log::info() << "WARNING Could not create the folder \e[1mpolygon_intersection\e[0m." << std::endl; } - geo_err_l1 += std::abs(diff_cell); - geo_err_linf = std::max(geo_err_linf, std::abs(diff_cell)); - } - ATLAS_TRACE_MPI(ALLREDUCE) { - mpi::comm().allReduceInPlace(geo_err_l1, eckit::mpi::sum()); - mpi::comm().allReduceInPlace(geo_err_linf, eckit::mpi::max()); } - remap_stat_.errors[Statistics::Errors::GEO_L1] = geo_err_l1 / unit_sphere_area(); - remap_stat_.errors[Statistics::Errors::GEO_LINF] = geo_err_linf; if (validate_) { + double geo_err_l1 = 0.; + double geo_err_linf = -1.; + for (idx_t scell = 0; scell < src_csp.size(); ++scell) { + if (std::get<1>(src_csp[scell]) != 0 ) { + // skip periodic & halo cells + continue; + } + double diff_cell = std::get<0>(src_csp[scell]).area(); + for (idx_t icell = 0; icell < src_iparam_[scell].weights.size(); ++icell) { + diff_cell -= src_iparam_[scell].weights[icell]; + } + geo_err_l1 += std::abs(diff_cell); + geo_err_linf = std::max(geo_err_linf, std::abs(diff_cell)); + } + ATLAS_TRACE_MPI(ALLREDUCE) { + mpi::comm().allReduceInPlace(geo_err_l1, eckit::mpi::sum()); + mpi::comm().allReduceInPlace(geo_err_linf, eckit::mpi::max()); + } + if (mpi::size() == 1) { + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_L1] = geo_err_l1 / unit_sphere_area(); + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_LINF] = geo_err_linf; + } + else { + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_L1] = -1111.; + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_LINF] = -1111.; + } + + std::fstream polygon_intersection_info("polygon_intersection", std::ios::out); + + geo_err_l1 = 0.; + geo_err_linf = -1.; for (idx_t tcell = 0; tcell < tgt_csp.size(); ++tcell) { + if (std::get<1>(tgt_csp[tcell]) != 0) { + // skip periodic & halo cells + continue; + } const auto& t_csp = std::get<0>(tgt_csp[tcell]); + + double tgt_cover_area = 0.; const auto& tiparam = tgt_iparam[tcell]; +#if 0 + // dump polygons in json format + if( tcell == 27609 ) { + dump_polygons_to_json(t_csp, src_csp, tiparam.cell_idx, 1.e-16); + } +#endif for (idx_t icell = 0; icell < tiparam.cell_idx.size(); ++icell) { tgt_cover_area += tiparam.tgt_weights[icell]; } - const double tgt_cover_err_percent = 100. * std::abs(t_csp.area() - tgt_cover_area) / t_csp.area(); - if (tgt_cover_err_percent > 0.1 and std::get<1>(tgt_csp[tcell]) == 0) { - Log::info() << "WARNING tgt cell covering error : " << tgt_cover_err_percent << " %\n"; - dump_intersection(t_csp, src_csp, tiparam.cell_idx); + /* + // TODO: normalise to target cell + double normm = t_csp.area() / (tgt_cover_area > 0. ? tgt_cover_area : t_csp.area()); + for (idx_t icell = 0; icell < tiparam.cell_idx.size(); ++icell) { + idx_t scell = tiparam.cell_idx[icell]; + auto siparam = src_iparam_[scell]; + size_t tgt_intersectors = siparam.cell_idx.size(); + for (idx_t sicell = 0; sicell < tgt_intersectors; sicell++ ) { + if (siparam.cell_idx[icell] == tcell) {; + siparam.weights[icell] *= normm; + siparam.tgt_weights[icell] *= normm; + } + } } + */ + double diff_cell = tgt_cover_area - t_csp.area(); + geo_err_l1 += std::abs(diff_cell); + geo_err_linf = std::max(geo_err_linf, std::abs(diff_cell)); + const double tgt_cover_err = 100. * diff_cell / t_csp.area(); + if (worst_tgt_overcover.second < tgt_cover_err) { + worst_tgt_overcover.second = tgt_cover_err;; + worst_tgt_overcover.first = tcell; + } + if (worst_tgt_undercover.second > tgt_cover_err) { + worst_tgt_undercover.second = tgt_cover_err;; + worst_tgt_undercover.first = tcell; + } + if (std::abs(tgt_cover_err) > 0.5) { + PointLonLat centre_ll; + eckit::geometry::Sphere::convertCartesianToSpherical(1., t_csp.centroid(), centre_ll); + polygon_intersection_info << "WARNING tgt cell " << tcell << " over-covering: \e[1m" << tgt_cover_err << "\e[0m %, cell-centre: " + << centre_ll <<"\n"; + polygon_intersection_info << "source indices: " << tiparam.cell_idx << std::endl; + dump_intersection("Target cell not exaclty covered", t_csp, src_csp, tiparam.cell_idx); + //dump_polygons_to_json(t_csp, src_csp, tiparam.cell_idx, "bad_polygon", 1.e-16); + } + } + ATLAS_TRACE_MPI(ALLREDUCE) { + mpi::comm().allReduceInPlace(geo_err_l1, eckit::mpi::sum()); + mpi::comm().allReduceInPlace(geo_err_linf, eckit::mpi::max()); + } + remap_stat_.errors[Statistics::Errors::TGT_INTERSECTPLG_L1] = geo_err_l1 / unit_sphere_area(); + remap_stat_.errors[Statistics::Errors::TGT_INTERSECTPLG_LINF] = geo_err_linf; + } + else { + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_L1] = -1111.; + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_LINF] = -1111.; + remap_stat_.errors[Statistics::Errors::TGT_INTERSECTPLG_L1] = -1111.; + remap_stat_.errors[Statistics::Errors::TGT_INTERSECTPLG_LINF] = -1111.; + } + + if (validate_) { + std::vector first(mpi::comm().size()); + std::vector second(mpi::comm().size()); + ATLAS_TRACE_MPI(ALLGATHER) { + mpi::comm().allGather(worst_tgt_overcover.first, first.begin(), first.end()); + mpi::comm().allGather(worst_tgt_overcover.second, second.begin(), second.end()); + } + auto max_over = std::max_element(second.begin(), second.end()); + auto rank_over = std::distance(second.begin(), max_over); + Log::info() << "WARNING The worst target polygon over-coveraging: \e[1m" + << *max_over + << "\e[0m %. For details, check the file: worst_target_cell_overcover.info " << std::endl; + if (rank_over == mpi::rank()) { + auto tcell = worst_tgt_overcover.first; + dump_polygons_to_json(std::get<0>(tgt_csp[tcell]), pointsSameEPS, src_csp, tgt_iparam[tcell].cell_idx, polygon_intersection_folder, "worst_target_cell_overcover"); + } + ATLAS_TRACE_MPI(ALLGATHER) { + mpi::comm().allGather(worst_tgt_undercover.first, first.begin(), first.end()); + mpi::comm().allGather(worst_tgt_undercover.second, second.begin(), second.end()); + } + auto min_under = std::min_element(second.begin(), second.end()); + auto rank_under = std::distance(second.begin(), min_under); + Log::info() << "WARNING The worst target polygon under-coveraging: \e[1m" + << *min_under + << "\e[0m %. For details, check the file: worst_target_cell_undercover.info " << std::endl; + + if (rank_under == mpi::rank()) { + auto tcell = worst_tgt_undercover.first; + dump_polygons_to_json(std::get<0>(tgt_csp[tcell]), pointsSameEPS, src_csp, tgt_iparam[tcell].cell_idx, polygon_intersection_folder, "worst_target_cell_undercover"); } } } @@ -943,7 +1255,7 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1 // determine the size of array of triplets used to define the sparse matrix if (src_cell_data_) { for (idx_t scell = 0; scell < n_spoints_; ++scell) { - triplets_size += src_iparam_[scell].centroids.size(); + triplets_size += src_iparam_[scell].cell_idx.size(); } } else { @@ -962,7 +1274,7 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1 if (src_cell_data_ && tgt_cell_data_) { for (idx_t scell = 0; scell < n_spoints_; ++scell) { const auto& iparam = src_iparam_[scell]; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; triplets.emplace_back(tcell, scell, iparam.tgt_weights[icell]); } @@ -974,8 +1286,9 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1 for (idx_t isubcell = 0; isubcell < src_node2csp_[snode].size(); ++isubcell) { const idx_t subcell = src_node2csp_[snode][isubcell]; const auto& iparam = src_iparam_[subcell]; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; + ATLAS_ASSERT(tcell < n_tpoints_); triplets.emplace_back(tcell, snode, iparam.tgt_weights[icell]); } } @@ -985,11 +1298,15 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1 auto& tgt_csp2node_ = data_->tgt_csp2node_; for (idx_t scell = 0; scell < n_spoints_; ++scell) { const auto& iparam = src_iparam_[scell]; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; + if (tnode >= n_tpoints_) { + Log::info() << "tnode, n_tpoints = " << tnode << ", " << n_tpoints_ << std::endl; + ATLAS_ASSERT(false); + } double inv_node_weight = (tgt_areas_v[tnode] > 0. ? 1. / tgt_areas_v[tnode] : 0.); - triplets.emplace_back(tnode, scell, iparam.src_weights[icell] * inv_node_weight); + triplets.emplace_back(tnode, scell, iparam.weights[icell] * inv_node_weight); } } } @@ -1000,16 +1317,42 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1 for (idx_t isubcell = 0; isubcell < src_node2csp_[snode].size(); ++isubcell) { const idx_t subcell = src_node2csp_[snode][isubcell]; const auto& iparam = src_iparam_[subcell]; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; + ATLAS_ASSERT(tnode < n_tpoints_); double inv_node_weight = (tgt_areas_v[tnode] > 0. ? 1. / tgt_areas_v[tnode] : 0.); - triplets.emplace_back(tnode, snode, iparam.src_weights[icell] * inv_node_weight); + triplets.emplace_back(tnode, snode, iparam.weights[icell] * inv_node_weight); + if( tnode >= n_tpoints_) { + Log::info() << tnode << " = tnode, " << n_tpoints_ << " = n_tpoints\n";; + ATLAS_ASSERT(false); + } } } } } sort_and_accumulate_triplets(triplets); // Very expensive!!! (90% of this routine). We need to avoid it + +// if (validate_) { +// std::vector weight_sum(n_tpoints_); +// for( auto& triplet : triplets ) { +// weight_sum[triplet.row()] += triplet.value(); +// } +// if (order_ == 1) { +// // first order should not give overshoots +// double eps = 1e4 * std::numeric_limits::epsilon(); +// for( auto& triplet : triplets ) { +// if (triplet.value() > 1. + eps or triplet.value() < -eps) { +// Log::info() << "target point " << triplet.row() << " weight: " << triplet.value() << std::endl; +// } +// } +// } +// for( size_t row=0; row < n_tpoints_; ++row ) { +// if (std::abs(weight_sum[row] - 1.) > 1e-11) { +// Log::info() << "target weight in row " << row << " differs from 1 by " << std::abs(weight_sum[row] - 1.) << std::endl; +// } +// } +// } return Matrix(n_tpoints_, n_spoints_, triplets); } @@ -1026,28 +1369,30 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 const auto src_halo = array::make_view(src_mesh_.cells().halo()); for (idx_t scell = 0; scell < n_spoints_; ++scell) { const auto nb_cells = get_cell_neighbours(src_mesh_, scell); - triplets_size += (2 * nb_cells.size() + 1) * src_iparam_[scell].centroids.size(); + triplets_size += (2 * nb_cells.size() + 1) * src_iparam_[scell].cell_idx.size(); } triplets.reserve(triplets_size); for (idx_t scell = 0; scell < n_spoints_; ++scell) { - const auto nb_cells = get_cell_neighbours(src_mesh_, scell); const auto& iparam = src_iparam_[scell]; - if (iparam.centroids.size() == 0 && not src_halo(scell)) { + if (iparam.cell_idx.size() == 0 && not src_halo(scell)) { continue; } /* // better conservation after Kritsikis et al. (2017) + // NOTE: ommited here at cost of conservation due to more involved implementation in parallel + // TEMP: use barycentre computed based on the source polygon's vertices, rather then + // the numerical barycentre based on barycentres of the intersections with target cells. + // TODO: for a given source cell, collect the centroids of all its intersections with target cells + // to compute the numerical barycentre of the cell bases on intersection. PointXYZ Cs = {0., 0., 0.}; - for ( idx_t icell = 0; icell < iparam.centroids.size(); ++icell ) { - Cs = Cs + PointXYZ::mul( iparam.centroids[icell], iparam.src_weights[icell] ); + for ( idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell ) { + Cs = Cs + PointXYZ::mul( iparam.centroids[icell], iparam.weights[icell] ); } - const double Cs_norm = PointXYZ::norm( Cs ); - ATLAS_ASSERT( Cs_norm > 0. ); - Cs = PointXYZ::div( Cs, Cs_norm ); */ const PointXYZ& Cs = src_points_[scell]; // compute gradient from cells double dual_area_inv = 0.; std::vector Rsj; + const auto nb_cells = get_cell_neighbours(src_mesh_, scell); Rsj.resize(nb_cells.size()); for (idx_t j = 0; j < nb_cells.size(); ++j) { idx_t nj = next_index(j, nb_cells.size()); @@ -1071,15 +1416,15 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 } // assemble the matrix std::vector Aik; - Aik.resize(iparam.centroids.size()); - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + Aik.resize(iparam.cell_idx.size()); + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { const PointXYZ& Csk = iparam.centroids[icell]; const PointXYZ Csk_Cs = Csk - Cs; Aik[icell] = Csk_Cs - PointXYZ::mul(Cs, PointXYZ::dot(Cs, Csk_Cs)); Aik[icell] = PointXYZ::mul(Aik[icell], iparam.tgt_weights[icell] * dual_area_inv); } if (tgt_cell_data_) { - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { const idx_t tcell = iparam.cell_idx[icell]; for (idx_t j = 0; j < nb_cells.size(); ++j) { idx_t nj = next_index(j, nb_cells.size()); @@ -1093,11 +1438,11 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 } else { auto& tgt_csp2node_ = data_->tgt_csp2node_; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; double inv_node_weight = (tgt_areas_v[tnode] > 0.) ? 1. / tgt_areas_v[tnode] : 0.; - double csp2node_coef = iparam.src_weights[icell] / iparam.tgt_weights[icell] * inv_node_weight; + double csp2node_coef = iparam.weights[icell] / iparam.tgt_weights[icell] * inv_node_weight; for (idx_t j = 0; j < nb_cells.size(); ++j) { idx_t nj = next_index(j, nb_cells.size()); idx_t sj = nb_cells[j]; @@ -1118,7 +1463,7 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 const auto nb_nodes = get_node_neighbours(src_mesh_, snode); for (idx_t isubcell = 0; isubcell < src_node2csp_[snode].size(); ++isubcell) { idx_t subcell = src_node2csp_[snode][isubcell]; - triplets_size += (2 * nb_nodes.size() + 1) * src_iparam_[subcell].centroids.size(); + triplets_size += (2 * nb_nodes.size() + 1) * src_iparam_[subcell].cell_idx.size(); } } triplets.reserve(triplets_size); @@ -1130,14 +1475,14 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 for ( idx_t isubcell = 0; isubcell < src_node2csp_[snode].size(); ++isubcell ) { idx_t subcell = src_node2csp_[snode][isubcell]; const auto& iparam = src_iparam_[subcell]; - for ( idx_t icell = 0; icell < iparam.centroids.size(); ++icell ) { - Cs = Cs + PointXYZ::mul( iparam.centroids[icell], iparam.src_weights[icell] ); + for ( idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell ) { + Cs = Cs + PointXYZ::mul( iparam.centroids[icell], iparam.weights[icell] ); } + const double Cs_norm = PointXYZ::norm( Cs ); + ATLAS_ASSERT( Cs_norm > 0. ); + Cs = PointXYZ::div( Cs, Cs_norm ); } - const double Cs_norm = PointXYZ::norm( Cs ); - ATLAS_ASSERT( Cs_norm > 0. ); - Cs = PointXYZ::div( Cs, Cs_norm ); -*/ + */ const PointXYZ& Cs = src_points_[snode]; // compute gradient from nodes double dual_area_inv = 0.; @@ -1168,19 +1513,19 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 for (idx_t isubcell = 0; isubcell < src_node2csp_[snode].size(); ++isubcell) { idx_t subcell = src_node2csp_[snode][isubcell]; const auto& iparam = src_iparam_[subcell]; - if (iparam.centroids.size() == 0) { + if (iparam.cell_idx.size() == 0) { continue; } std::vector Aik; - Aik.resize(iparam.centroids.size()); - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + Aik.resize(iparam.cell_idx.size()); + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { const PointXYZ& Csk = iparam.centroids[icell]; const PointXYZ Csk_Cs = Csk - Cs; Aik[icell] = Csk_Cs - PointXYZ::mul(Cs, PointXYZ::dot(Cs, Csk_Cs)); Aik[icell] = PointXYZ::mul(Aik[icell], iparam.tgt_weights[icell] * dual_area_inv); } if (tgt_cell_data_) { - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { const idx_t tcell = iparam.cell_idx[icell]; for (idx_t j = 0; j < nb_nodes.size(); ++j) { idx_t nj = next_index(j, nb_nodes.size()); @@ -1194,11 +1539,11 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 } else { auto& tgt_csp2node_ = data_->tgt_csp2node_; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; double inv_node_weight = (tgt_areas_v[tnode] > 1e-15) ? 1. / tgt_areas_v[tnode] : 0.; - double csp2node_coef = iparam.src_weights[icell] / iparam.tgt_weights[icell] * inv_node_weight; + double csp2node_coef = iparam.weights[icell] / iparam.tgt_weights[icell] * inv_node_weight; for (idx_t j = 0; j < nb_nodes.size(); ++j) { idx_t nj = next_index(j, nb_nodes.size()); idx_t sj = nb_nodes[j]; @@ -1218,6 +1563,16 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 return Matrix(n_tpoints_, n_spoints_, triplets); } +void ConservativeSphericalPolygonInterpolation::do_execute(const FieldSet& src_fields, FieldSet& tgt_fields, + Metadata& metadata) const { + std::vector md_array; + md_array.resize( src_fields.size() ); + for (int i = 0; i < src_fields.size(); i++) { // TODO: memory-wise should we openmp this? + do_execute(src_fields[i], tgt_fields[i], md_array[i]); + } + metadata = md_array[0]; // TODO: reduce metadata of a set of variables to a single metadata of a variable? +} + void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_field, Field& tgt_field, Metadata& metadata) const { ATLAS_TRACE("ConservativeMethod::do_execute()"); @@ -1245,8 +1600,8 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel } for (idx_t scell = 0; scell < src_vals.size(); ++scell) { const auto& iparam = src_iparam_[scell]; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { - tgt_vals(iparam.cell_idx[icell]) += iparam.src_weights[icell] * src_vals(scell); + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { + tgt_vals(iparam.cell_idx[icell]) += iparam.weights[icell] * src_vals(scell); } } for (idx_t tcell = 0; tcell < tgt_vals.size(); ++tcell) { @@ -1260,16 +1615,13 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel } else if (order_ == 2) { if (matrix_free_) { - ATLAS_TRACE("matrix_free_order_2"); - const auto& src_iparam_ = data_->src_iparam_; - const auto& tgt_areas_v = data_->tgt_areas_; - if (not src_cell_data_ or not tgt_cell_data_) { ATLAS_NOTIMPLEMENTED; } - + ATLAS_TRACE("matrix_free_order_2"); + const auto& src_iparam_ = data_->src_iparam_; + const auto& tgt_areas_v = data_->tgt_areas_; auto& src_points_ = data_->src_points_; - const auto src_vals = array::make_view(src_field); auto tgt_vals = array::make_view(tgt_field); const auto halo = array::make_view(src_mesh_.cells().halo()); @@ -1277,55 +1629,42 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel tgt_vals(tcell) = 0.; } for (idx_t scell = 0; scell < src_vals.size(); ++scell) { - if (halo(scell)) { + const auto& iparam = src_iparam_[scell]; + if (iparam.cell_idx.size() == 0 or halo(scell)) { continue; } - const auto& iparam = src_iparam_[scell]; - const PointXYZ& P = src_points_[scell]; + const PointXYZ& Cs = src_points_[scell]; PointXYZ grad = {0., 0., 0.}; PointXYZ src_barycenter = {0., 0., 0.}; - auto src_neighbour_cells = get_cell_neighbours(src_mesh_, scell); - double dual_area = 0.; - for (idx_t nb_id = 0; nb_id < src_neighbour_cells.size(); ++nb_id) { - idx_t nnb_id = next_index(nb_id, src_neighbour_cells.size()); - idx_t ncell = src_neighbour_cells[nb_id]; - idx_t nncell = src_neighbour_cells[nnb_id]; - const auto& Pn = src_points_[ncell]; - const auto& Pnn = src_points_[nncell]; - if (ncell != scell && nncell != scell) { - double val = 0.5 * (src_vals(ncell) + src_vals(nncell)) - src_vals(scell); - auto csp = ConvexSphericalPolygon({Pn, Pnn, P}); - if (csp.area() < std::numeric_limits::epsilon()) { - csp = ConvexSphericalPolygon({Pn, P, Pnn}); - } - auto NsNsj = ConvexSphericalPolygon::GreatCircleSegment(P, Pn); - val *= (NsNsj.inLeftHemisphere(Pnn, -1e-16) ? -1 : 1); - dual_area += std::abs(csp.area()); - grad = grad + PointXYZ::mul(PointXYZ::cross(Pn, Pnn), val); - } - else if (ncell != scell) { - ATLAS_NOTIMPLEMENTED; - //double val = 0.5 * ( src_vals( ncell ) - src_vals( scell ) ); - //grad = grad + PointXYZ::mul( PointXYZ::cross( Pn, P ), val ); + auto nb_cells = get_cell_neighbours(src_mesh_, scell); + double dual_area_inv = 0.; + for (idx_t j = 0; j < nb_cells.size(); ++j) { + idx_t nj = next_index(j, nb_cells.size()); + idx_t sj = nb_cells[j]; + idx_t nsj = nb_cells[nj]; + const auto& Csj = src_points_[sj]; + const auto& Cnsj = src_points_[nsj]; + double val = 0.5 * (src_vals(nj) + src_vals(nsj)) - src_vals(j); + if (ConvexSphericalPolygon::GreatCircleSegment(Cs, Csj).inLeftHemisphere(Cnsj, -1e-16)) { + dual_area_inv += ConvexSphericalPolygon({Cs, Csj, Cnsj}).area(); } - else if (nncell != scell) { - ATLAS_NOTIMPLEMENTED; - //double val = 0.5 * ( src_vals( nncell ) - src_vals( scell ) ); - //grad = grad + PointXYZ::mul( PointXYZ::cross( P, Pnn ), val ); + else { + //val *= -1; + dual_area_inv += ConvexSphericalPolygon({Cs, Cnsj, Csj}).area(); } + grad = grad + PointXYZ::mul(PointXYZ::cross(Csj, Cnsj), val); } - if (dual_area > std::numeric_limits::epsilon()) { - grad = PointXYZ::div(grad, dual_area); - } - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { - src_barycenter = src_barycenter + PointXYZ::mul(iparam.centroids[icell], iparam.src_weights[icell]); + dual_area_inv = ((dual_area_inv > 0.) ? 1. / dual_area_inv : 0.); + grad = PointXYZ::mul(grad, dual_area_inv); + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { + src_barycenter = src_barycenter + PointXYZ::mul(iparam.centroids[icell], iparam.weights[icell]); } src_barycenter = PointXYZ::div(src_barycenter, PointXYZ::norm(src_barycenter)); grad = grad - PointXYZ::mul(src_barycenter, PointXYZ::dot(grad, src_barycenter)); ATLAS_ASSERT(std::abs(PointXYZ::dot(grad, src_barycenter)) < 1e-14); - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { tgt_vals(iparam.cell_idx[icell]) += - iparam.src_weights[icell] * + iparam.weights[icell] * (src_vals(scell) + PointXYZ::dot(grad, iparam.centroids[icell] - src_barycenter)); } } @@ -1340,6 +1679,10 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel } stopwatch.stop(); + { + ATLAS_TRACE("halo exchange target"); + tgt_field.haloExchange(); + } auto remap_stat = remap_stat_; if (statistics_conservation_) { @@ -1369,7 +1712,7 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel else { auto& src_node2csp_ = data_->src_node2csp_; for (idx_t spt = 0; spt < src_vals.size(); ++spt) { - if (src_node_ghost(spt) or src_areas_v[spt] < 1e-14) { + if (src_node_halo(spt) or src_node_ghost(spt)) { continue; } err_remap_cons += src_vals(spt) * src_areas_v[spt]; @@ -1386,24 +1729,28 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel } else { for (idx_t tpt = 0; tpt < tgt_vals.size(); ++tpt) { - if (tgt_node_ghost(tpt)) { + if (tgt_node_halo(tpt) or tgt_node_ghost(tpt)) { continue; } err_remap_cons -= tgt_vals(tpt) * tgt_areas_v[tpt]; } } ATLAS_TRACE_MPI(ALLREDUCE) { mpi::comm().allReduceInPlace(&err_remap_cons, 1, eckit::mpi::sum()); } - remap_stat.errors[Statistics::Errors::REMAP_CONS] = std::sqrt(std::abs(err_remap_cons) / unit_sphere_area()); + remap_stat.errors[Statistics::Errors::REMAP_CONS] = err_remap_cons / unit_sphere_area(); metadata.set("conservation_error", remap_stat.errors[Statistics::Errors::REMAP_CONS]); } if (statistics_intersection_) { - metadata.set("polygons.source", remap_stat.counts[Statistics::SRC_PLG]); - metadata.set("polygons.target", remap_stat.counts[Statistics::TGT_PLG]); - metadata.set("polygons.intersections", remap_stat.counts[Statistics::INT_PLG]); + metadata.set("polygons.source", remap_stat.counts[Statistics::Counts::SRC_PLG]); + metadata.set("polygons.target", remap_stat.counts[Statistics::Counts::TGT_PLG]); + metadata.set("polygons.intersections", remap_stat.counts[Statistics::Counts::INT_PLG]); metadata.set("polygons.uncovered_source", remap_stat.counts[Statistics::UNCVR_SRC]); - metadata.set("source_area_error.L1", remap_stat.errors[Statistics::Errors::GEO_L1]); - metadata.set("source_area_error.Linf", remap_stat.errors[Statistics::Errors::GEO_LINF]); + if (validate_) { + metadata.set("source_area_error.L1", remap_stat.errors[Statistics::Errors::SRC_INTERSECTPLG_L1]); + metadata.set("source_area_error.Linf", remap_stat.errors[Statistics::Errors::SRC_INTERSECTPLG_LINF]); + metadata.set("target_area_error.L1", remap_stat.errors[Statistics::Errors::TGT_INTERSECTPLG_L1]); + metadata.set("target_area_error.Linf", remap_stat.errors[Statistics::Errors::TGT_INTERSECTPLG_LINF]); + } } if (statistics_intersection_ || statistics_conservation_) { @@ -1430,15 +1777,16 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel metadata.set("memory.src_node2csp", memory_of(data_->src_node2csp_)); metadata.set("memory.tgt_node2csp", memory_of(data_->tgt_node2csp_)); metadata.set("memory.src_iparam", memory_of(data_->src_iparam_)); - - tgt_field.set_dirty(); } void ConservativeSphericalPolygonInterpolation::print(std::ostream& out) const { out << "ConservativeMethod{"; out << "order:" << order_; - out << ", source:" << (src_cell_data_ ? "cells" : "nodes"); - out << ", target:" << (tgt_cell_data_ ? "cells" : "nodes"); + int halo = 0; + src_mesh_.metadata().get("halo", halo); + out << ", source:" << (src_cell_data_ ? "cells(" : "nodes(") << src_mesh_.grid().name() << ",halo=" << halo << ")"; + tgt_mesh_.metadata().get("halo", halo); + out << ", target:" << (tgt_cell_data_ ? "cells(" : "nodes(") << tgt_mesh_.grid().name() << ",halo=" << halo << ")"; out << ", normalise_intersections:" << normalise_intersections_; out << ", matrix_free:" << matrix_free_; out << ", statistics.intersection:" << statistics_intersection_; @@ -1506,7 +1854,7 @@ void ConservativeSphericalPolygonInterpolation::setup_stat() const { remap_stat_.tgt_area_sum = src_tgt_sums[1]; geo_create_err = std::abs(src_tgt_sums[0] - src_tgt_sums[1]) / unit_sphere_area(); - remap_stat_.errors[Statistics::Errors::GEO_DIFF] = geo_create_err; + remap_stat_.errors[Statistics::Errors::SRCTGT_INTERSECTPLG_DIFF] = geo_create_err; } Field ConservativeSphericalPolygonInterpolation::Statistics::diff(const Interpolation& interpolation, @@ -1543,29 +1891,34 @@ Field ConservativeSphericalPolygonInterpolation::Statistics::diff(const Interpol double diff = src_vals(spt) * src_areas_v[spt]; const auto& iparam = src_iparam_[spt]; if (tgt_cell_data_) { - for (idx_t icell = 0; icell < iparam.src_weights.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.weights.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; if (tgt_cell_halo(tcell) < 1) { - diff -= tgt_vals(iparam.cell_idx[icell]) * iparam.src_weights[icell]; + diff -= tgt_vals(tcell) * iparam.weights[icell]; } } } else { - for (idx_t icell = 0; icell < iparam.src_weights.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.weights.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; - if (tgt_node_halo(tnode) < 1) { - diff -= tgt_vals(tnode) * iparam.src_weights[icell]; + if (not tgt_node_ghost(tnode)) { + diff -= tgt_vals(tnode) * iparam.weights[icell]; } } } - diff_vals(spt) = std::abs(diff) / src_areas_v[spt]; + if (src_areas_v[spt] > 0.) { + diff_vals(spt) = diff / src_areas_v[spt]; + } + else { + Log::info() << " at cell " << spt << " cell-area: " << src_areas_v[spt] << std::endl; + diff_vals(spt) = std::numeric_limits::max(); + } } } else { for (idx_t spt = 0; spt < src_vals.size(); ++spt) { - if (src_node_ghost(spt) or src_areas_v[spt] < 1e-14) { - diff_vals(spt) = 0.; + if (src_node_ghost(spt) or src_node_halo(spt) != 0) { continue; } double diff = src_vals(spt) * src_areas_v[spt]; @@ -1573,38 +1926,50 @@ Field ConservativeSphericalPolygonInterpolation::Statistics::diff(const Interpol for (idx_t subcell = 0; subcell < node2csp.size(); ++subcell) { const auto& iparam = src_iparam_[node2csp[subcell]]; if (tgt_cell_data_) { - for (idx_t icell = 0; icell < iparam.src_weights.size(); ++icell) { - diff -= tgt_vals(iparam.cell_idx[icell]) * iparam.src_weights[icell]; + for (idx_t icell = 0; icell < iparam.weights.size(); ++icell) { + idx_t tcell = iparam.cell_idx[icell]; + if (tgt_cell_halo(tcell) < 1) { + diff -= tgt_vals(tcell) * iparam.weights[icell]; + } } } else { - for (idx_t icell = 0; icell < iparam.src_weights.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.weights.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; - diff -= tgt_vals(tnode) * iparam.src_weights[icell]; + if (not tgt_node_ghost(tnode)) { + diff -= tgt_vals(tnode) * iparam.weights[icell]; + } } } } - diff_vals(spt) = std::abs(diff) / src_areas_v[spt]; + if (src_areas_v[spt] > 0.) { + diff_vals(spt) = diff / src_areas_v[spt]; + } + else { + diff_vals(spt) = std::numeric_limits::max(); + Log::info() << " at cell " << spt << " cell-area: " << src_areas_v[spt] << std::endl; + } } } return diff; } -void ConservativeSphericalPolygonInterpolation::Statistics::accuracy(const Interpolation& interpolation, - const Field target, - std::function func) { +ConservativeSphericalPolygonInterpolation::Metadata ConservativeSphericalPolygonInterpolation::Statistics::accuracy(const Interpolation& interpolation, + const Field target, + std::function func) { auto tgt_vals = array::make_view(target); auto cachable_data_ = ConservativeSphericalPolygonInterpolation::Cache(interpolation).get(); auto tgt_mesh_ = extract_mesh(cachable_data_->src_fs_); auto tgt_cell_data_ = extract_mesh(cachable_data_->tgt_fs_); const auto tgt_cell_halo = array::make_view(tgt_mesh_.cells().halo()); const auto tgt_node_ghost = array::make_view(tgt_mesh_.nodes().ghost()); + const auto tgt_node_halo = array::make_view(tgt_mesh_.nodes().halo()); const auto& tgt_areas_v = cachable_data_->tgt_areas_; + auto& tgt_points_ = cachable_data_->tgt_points_; double err_remap_l2 = 0.; double err_remap_linf = 0.; - auto& tgt_points_ = cachable_data_->tgt_points_; if (tgt_cell_data_) { size_t ncells = std::min(tgt_vals.size(), tgt_mesh_.cells().size()); for (idx_t tpt = 0; tpt < ncells; ++tpt) { @@ -1623,7 +1988,7 @@ void ConservativeSphericalPolygonInterpolation::Statistics::accuracy(const Inter else { size_t nnodes = std::min(tgt_vals.size(), tgt_mesh_.nodes().size()); for (idx_t tpt = 0; tpt < nnodes; ++tpt) { - if (tgt_node_ghost(tpt)) { + if (tgt_node_ghost(tpt) or tgt_node_halo(tpt)) { continue; } auto p = tgt_points_[tpt]; @@ -1638,77 +2003,71 @@ void ConservativeSphericalPolygonInterpolation::Statistics::accuracy(const Inter mpi::comm().allReduceInPlace(&err_remap_l2, 1, eckit::mpi::sum()); mpi::comm().allReduceInPlace(&err_remap_linf, 1, eckit::mpi::max()); } - this->errors[Statistics::Errors::REMAP_L2] = std::sqrt(err_remap_l2 / unit_sphere_area()); - this->errors[Statistics::Errors::REMAP_LINF] = err_remap_linf; + errors[Statistics::Errors::REMAP_L2] = std::sqrt(err_remap_l2 / unit_sphere_area()); + errors[Statistics::Errors::REMAP_LINF] = err_remap_linf; + Metadata metadata; + metadata.set("errors.REMAP_L2", errors[Statistics::Errors::REMAP_L2]); + metadata.set("errors.REMAP_LINF", errors[Statistics::Errors::REMAP_LINF]); + return metadata; } -auto debug_intersection = [](const ConvexSphericalPolygon& plg_1, const ConvexSphericalPolygon& plg_2, - const ConvexSphericalPolygon& iplg, const ConvexSphericalPolygon& jplg) { - const double intersection_comm_err = std::abs(iplg.area() - jplg.area()) / (plg_1.area() > 0 ? plg_1.area() : 1.); - Log::info().indent(); - if (intersection_comm_err > 1e-6) { - Log::info() << "PLG_1 : " << std::setprecision(10) << plg_1 << "\n"; - Log::info() << "area(PLG_1) : " << plg_1.area() << "\n"; - Log::info() << "PLG_2 :" << plg_2 << "\n"; - Log::info() << "PLG_12 : " << iplg << "\n"; - Log::info() << "PLG_21 : " << jplg << "\n"; - Log::info() << "area(PLG_12 - PLG_21) : " << intersection_comm_err << "\n"; - Log::info() << "area(PLG_21) : " << jplg.area() << "\n"; - //ATLAS_ASSERT( false, "SRC.intersect.TGT =/= TGT.intersect.SRC."); - } - int pout; - int pin = inside_vertices(plg_1, plg_2, pout); - if (pin > 2 && iplg.area() < 3e-16) { - Log::info() << " pin : " << pin << ", pout :" << pout << ", total vertices : " << plg_2.size() << "\n"; - Log::info() << "PLG_2 :" << plg_2 << "\n"; - Log::info() << "PLG_12 : " << iplg << "\n"; - Log::info() << "area(PLG_12) : " << iplg.area() << "\n\n"; - //ATLAS_ASSERT( false, "SRC must intersect TGT." ); - } - Log::info().unindent(); -}; - -void ConservativeSphericalPolygonInterpolation::dump_intersection(const ConvexSphericalPolygon& plg_1, +void ConservativeSphericalPolygonInterpolation::dump_intersection(const std::string msg, + const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, const std::vector& plg_2_idx_array) const { +#if PRINT_BAD_POLYGONS double plg_1_coverage = 0.; + std::vector int_areas; + int_areas.resize(plg_2_idx_array.size()); for (int i = 0; i < plg_2_idx_array.size(); ++i) { const auto plg_2_idx = plg_2_idx_array[i]; const auto& plg_2 = std::get<0>(plg_2_array[plg_2_idx]); - auto iplg = plg_1.intersect(plg_2); - auto jplg = plg_2.intersect(plg_1); - debug_intersection(plg_1, plg_2, iplg, jplg); + auto iplg = plg_1.intersect(plg_2, false, 1.e5 * std::numeric_limits::epsilon()); plg_1_coverage += iplg.area(); + int_areas[i] = iplg.area(); } - Log::info().indent(); if (std::abs(plg_1.area() - plg_1_coverage) > 0.01 * plg_1.area()) { - Log::info() << "Polygon coverage incomplete. Printing polygons." << std::endl; - Log::info() << "Polygon 1 : "; + Log::info().indent(); + Log::info() << msg << ", Polygon_1.area: " << plg_1.area() << ", covered: " << plg_1_coverage << std::endl; + Log::info() << "Polygon_1 : "; + Log::info().precision(18); plg_1.print(Log::info()); Log::info() << std::endl << "Printing " << plg_2_idx_array.size() << " covering polygons -->" << std::endl; Log::info().indent(); + plg_1_coverage = plg_1.area(); for (int i = 0; i < plg_2_idx_array.size(); ++i) { const auto plg_2_idx = plg_2_idx_array[i]; const auto& plg_2 = std::get<0>(plg_2_array[plg_2_idx]); - Log::info() << "Polygon " << i + 1 << " : "; plg_2.print(Log::info()); + Log::info() << "\ninters:"; + plg_1_coverage -= int_areas[i]; + auto iplg = plg_1.intersect(plg_2, false, 1.e5 * std::numeric_limits::epsilon()); + Log::info().indent(); + iplg.print(Log::info()); + Log::info() << ", inters.-area: " << iplg.area() << ", plg1-area left: " << plg_1_coverage << std::endl; + Log::info().unindent(); Log::info() << std::endl; } + Log::info() << std::endl; + Log::info().unindent(); Log::info().unindent(); } - Log::info().unindent(); +#endif } template -void ConservativeSphericalPolygonInterpolation::dump_intersection(const ConvexSphericalPolygon& plg_1, +void ConservativeSphericalPolygonInterpolation::dump_intersection(const std::string msg, + const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, const TargetCellsIDs& plg_2_idx_array) const { +#if PRINT_BAD_POLYGONS std::vector idx_array; idx_array.resize(plg_2_idx_array.size()); for (int i = 0; i < plg_2_idx_array.size(); ++i) { idx_array[i] = plg_2_idx_array[i].payload(); } - dump_intersection(plg_1, plg_2_array, idx_array); + dump_intersection(msg, plg_1, plg_2_array, idx_array); +#endif } ConservativeSphericalPolygonInterpolation::Cache::Cache(std::shared_ptr entry): @@ -1749,20 +2108,16 @@ void ConservativeSphericalPolygonInterpolation::Data::print(std::ostream& out) c } void ConservativeSphericalPolygonInterpolation::Statistics::fillMetadata(Metadata& metadata) { - // counts - metadata.set("counts.SRC_PLG", counts[SRC_PLG]); - metadata.set("counts.TGT_PLG", counts[TGT_PLG]); - metadata.set("counts.INT_PLG", counts[INT_PLG]); - metadata.set("counts.UNCVR_SRC", counts[UNCVR_SRC]); - // errors - metadata.set("errors.SRC_PLG_L1", errors[SRC_PLG_L1]); - metadata.set("errors.SRC_PLG_LINF", errors[SRC_PLG_LINF]); - metadata.set("errors.TGT_PLG_L1", errors[TGT_PLG_L1]); - metadata.set("errors.TGT_PLG_LINF", errors[TGT_PLG_LINF]); - metadata.set("errors.GEO_L1", errors[GEO_L1]); - metadata.set("errors.GEO_LINF", errors[GEO_LINF]); - metadata.set("errors.GEO_DIFF", errors[GEO_DIFF]); + metadata.set("errors.SRC_SUBPLG_L1", errors[SRC_SUBPLG_L1]); + metadata.set("errors.SRC_SUBPLG_LINF", errors[SRC_SUBPLG_LINF]); + metadata.set("errors.TGT_SUBPLG_L1", errors[TGT_SUBPLG_L1]); + metadata.set("errors.TGT_SUBPLG_LINF", errors[TGT_SUBPLG_LINF]); + metadata.set("errors.SRC_INTERSECTPLG_L1", errors[SRC_INTERSECTPLG_L1]); + metadata.set("errors.SRC_INTERSECTPLG_LINF", errors[SRC_INTERSECTPLG_LINF]); + metadata.set("errors.TGT_INTERSECTPLG_L1", errors[TGT_INTERSECTPLG_L1]); + metadata.set("errors.TGT_INTERSECTPLG_LINF", errors[TGT_INTERSECTPLG_LINF]); + metadata.set("errors.SRCTGT_INTERSECTPLG_DIFF", errors[SRCTGT_INTERSECTPLG_DIFF]); metadata.set("errors.REMAP_CONS", errors[REMAP_CONS]); metadata.set("errors.REMAP_L2", errors[REMAP_L2]); metadata.set("errors.REMAP_LINF", errors[REMAP_LINF]); @@ -1774,20 +2129,16 @@ ConservativeSphericalPolygonInterpolation::Statistics::Statistics() { } ConservativeSphericalPolygonInterpolation::Statistics::Statistics(const Metadata& metadata): Statistics() { - // counts - metadata.get("counts.SRC_PLG", counts[SRC_PLG]); - metadata.get("counts.TGT_PLG", counts[TGT_PLG]); - metadata.get("counts.INT_PLG", counts[INT_PLG]); - metadata.get("counts.UNCVR_SRC", counts[UNCVR_SRC]); - // errors - metadata.get("errors.SRC_PLG_L1", errors[SRC_PLG_L1]); - metadata.get("errors.SRC_PLG_LINF", errors[SRC_PLG_LINF]); - metadata.get("errors.TGT_PLG_L1", errors[TGT_PLG_L1]); - metadata.get("errors.TGT_PLG_LINF", errors[TGT_PLG_LINF]); - metadata.get("errors.GEO_L1", errors[GEO_L1]); - metadata.get("errors.GEO_LINF", errors[GEO_LINF]); - metadata.get("errors.GEO_DIFF", errors[GEO_DIFF]); + metadata.get("errors.SRC_SUBPLG_L1", errors[SRC_SUBPLG_L1]); + metadata.get("errors.SRC_SUBPLG_LINF", errors[SRC_SUBPLG_LINF]); + metadata.get("errors.TGT_SUBPLG_L1", errors[TGT_SUBPLG_L1]); + metadata.get("errors.TGT_SUBPLG_LINF", errors[TGT_SUBPLG_LINF]); + metadata.get("errors.SRC_INTERSECTPLG_L1", errors[SRC_INTERSECTPLG_L1]); + metadata.get("errors.SRC_INTERSECTPLG_LINF", errors[SRC_INTERSECTPLG_LINF]); + metadata.get("errors.TGT_INTERSECTPLG_L1", errors[TGT_INTERSECTPLG_L1]); + metadata.get("errors.TGT_INTERSECTPLG_LINF", errors[TGT_INTERSECTPLG_LINF]); + metadata.get("errors.SRCTGT_INTERSECTPLG_DIFF", errors[SRCTGT_INTERSECTPLG_DIFF]); metadata.get("errors.REMAP_CONS", errors[REMAP_CONS]); metadata.get("errors.REMAP_L2", errors[REMAP_L2]); metadata.get("errors.REMAP_LINF", errors[REMAP_LINF]); diff --git a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h index bcab9fdba..3dd79099b 100644 --- a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h +++ b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h @@ -25,10 +25,8 @@ class ConservativeSphericalPolygonInterpolation : public Method { struct InterpolationParameters { // one polygon intersection std::vector cell_idx; // target cells used for intersection std::vector centroids; // intersection cell centroids - std::vector src_weights; // intersection cell areas - - // TODO: tgt_weights can be computed on the fly - std::vector tgt_weights; // (intersection cell areas) / (target cell area) + std::vector weights; // intersection cell areas + std::vector tgt_weights; // (intersection cell areas) / (target cell area) // TODO: tgt_weights vector should be removed for the sake of highres }; private: @@ -55,7 +53,6 @@ class ConservativeSphericalPolygonInterpolation : public Method { std::vector> src_node2csp_; std::vector> tgt_node2csp_; - // Timings struct Timings { double source_polygons_assembly{0}; @@ -102,18 +99,20 @@ class ConservativeSphericalPolygonInterpolation : public Method { std::array counts; enum Errors { - SRC_PLG_L1 = 0, // index, over/undershoot in source subpolygon creation - SRC_PLG_LINF, - TGT_PLG_L1, // index, over/untershoot in target subpolygon creation - TGT_PLG_LINF, - GEO_L1, // index, cumulative area mismatch in polygon intersections - GEO_LINF, // index, like GEO_L1 but in L_infinity norm - GEO_DIFF, // index, difference in earth area coverages + SRC_SUBPLG_L1 = 0, // index, \sum_{cell of mesh} {cell.area - \sum_{subpol of cell} subpol.area} + SRC_SUBPLG_LINF, // index, \max_-||- + TGT_SUBPLG_L1, // see above + TGT_SUBPLG_LINF, // see above + SRC_INTERSECTPLG_L1, // index, \sum_{cell of src-mesh} {cell.area - \sum_{intersect of cell} intersect.area} + SRC_INTERSECTPLG_LINF, // index, \max_-||- + TGT_INTERSECTPLG_L1, // see above + TGT_INTERSECTPLG_LINF, // see above + SRCTGT_INTERSECTPLG_DIFF, // index, 1/(unit_sphere.area) ( \sum_{scell} scell.area - \sum{tcell} tcell.area ) REMAP_CONS, // index, error in mass conservation REMAP_L2, // index, error accuracy for given analytical function REMAP_LINF // index, like REMAP_L2 but in L_infinity norm }; - std::array errors; + std::array errors; double tgt_area_sum; double src_area_sum; @@ -123,8 +122,8 @@ class ConservativeSphericalPolygonInterpolation : public Method { Statistics(); Statistics(const Metadata&); - void accuracy(const Interpolation& interpolation, const Field target, - std::function func); + Metadata accuracy(const Interpolation& interpolation, const Field target, + std::function func); // compute difference field of source and target mass @@ -139,6 +138,7 @@ class ConservativeSphericalPolygonInterpolation : public Method { void do_setup(const FunctionSpace& src_fs, const FunctionSpace& tgt_fs) override; void do_setup(const Grid& src_grid, const Grid& tgt_grid, const interpolation::Cache&) override; void do_execute(const Field& src_field, Field& tgt_field, Metadata&) const override; + void do_execute(const FieldSet& src_fields, FieldSet& tgt_fields, Metadata&) const override; void print(std::ostream& out) const override; @@ -160,17 +160,17 @@ class ConservativeSphericalPolygonInterpolation : public Method { void intersect_polygons(const CSPolygonArray& src_csp, const CSPolygonArray& tgt_scp); Matrix compute_1st_order_matrix(); Matrix compute_2nd_order_matrix(); - void dump_intersection(const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, + void dump_intersection(const std::string, const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, const std::vector& plg_2_idx_array) const; template - void dump_intersection(const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, + void dump_intersection(const std::string, const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, const TargetCellsIDs& plg_2_idx_array) const; std::vector sort_cell_edges(Mesh& mesh, idx_t cell_id) const; std::vector sort_node_edges(Mesh& mesh, idx_t cell_id) const; - std::vector get_cell_neighbours(Mesh& mesh, idx_t jcell) const; - std::vector get_node_neighbours(Mesh& mesh, idx_t jcell) const; - CSPolygonArray get_polygons_celldata(Mesh& mesh) const; - CSPolygonArray get_polygons_nodedata(Mesh& mesh, std::vector& csp2node, + std::vector get_cell_neighbours(Mesh&, idx_t jcell) const; + std::vector get_node_neighbours(Mesh&, idx_t jcell) const; + CSPolygonArray get_polygons_celldata(FunctionSpace) const; + CSPolygonArray get_polygons_nodedata(FunctionSpace, std::vector& csp2node, std::vector>& node2csp, std::array& errors) const; diff --git a/src/sandbox/interpolation/atlas-conservative-interpolation.cc b/src/sandbox/interpolation/atlas-conservative-interpolation.cc index 48b6e71eb..6e5199674 100644 --- a/src/sandbox/interpolation/atlas-conservative-interpolation.cc +++ b/src/sandbox/interpolation/atlas-conservative-interpolation.cc @@ -8,19 +8,6 @@ * nor does it submit to any jurisdiction. */ - -// TODO: -// ----- -// * Fix abort encountered with -// mpirun -np 4 atlas-conservative-interpolation --source.grid=O20 --target.grid=H8 --order=2 -// -// QUESTIONS: -// ---------- -// * Why sqrt in ConservativeSphericalPolygon in line -// remap_stat.errors[Statistics::Errors::REMAP_CONS] = std::sqrt(std::abs(err_remap_cons) / unit_sphere_area()); -// used to compute conservation_error - - #include #include #include @@ -35,6 +22,8 @@ #include "atlas/array/MakeView.h" #include "atlas/field.h" #include "atlas/grid.h" +#include "atlas/grid/Spacing.h" +#include "atlas/grid/detail/spacing/CustomSpacing.h" #include "atlas/interpolation/Interpolation.h" #include "atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h" #include "atlas/mesh.h" @@ -73,13 +62,15 @@ class AtlasParallelInterpolation : public AtlasTool { add_option(new eckit::option::Separator("Source/Target options")); add_option(new SimpleOption("source.grid", "source gridname")); + add_option(new SimpleOption("source.partitioner", "source partitioner name (spherical-polygon, lonlat-polygon, brute-force)")); add_option(new SimpleOption("target.grid", "target gridname")); + add_option(new SimpleOption("target.partitioner", "target partitioner name (equal_regions, regular_bands, equal_bands)")); add_option(new SimpleOption("source.functionspace", "source functionspace, to override source grid default")); add_option(new SimpleOption("target.functionspace", "target functionspace, to override target grid default")); add_option(new SimpleOption("source.halo", "default=2")); - add_option(new SimpleOption("target.halo", "default=0")); + add_option(new SimpleOption("target.halo", "default=0 for CellColumns and 1 for NodeColumns")); add_option(new eckit::option::Separator("Interpolation options")); add_option(new SimpleOption("order", "Interpolation order. Supported: 1, 2 (default=1)")); @@ -182,8 +173,22 @@ std::function get_init(const AtlasTool::Args& args) } int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { - auto src_grid = Grid{args.getString("source.grid", "H16")}; - auto tgt_grid = Grid{args.getString("target.grid", "H32")}; + auto get_grid = [](std::string grid_name) { + int grid_number = std::atoi( grid_name.substr(1, grid_name.size()).c_str() ); + if (grid_name.at(0) == 'P') { + Log::info() << "P-grid number: " << grid_number << std::endl; + ATLAS_ASSERT(grid_number > 3); + std::vector y = {90, 89.9999, 0, -90}; + auto xspace = StructuredGrid::XSpace( grid::LinearSpacing(0, 360, grid_number, false) ); + auto yspace = StructuredGrid::YSpace( new grid::spacing::CustomSpacing( y.size(), y.data() ) ); + return StructuredGrid(xspace, yspace); + } + else { + return StructuredGrid{grid_name}; + } + }; + auto src_grid = get_grid(args.getString("source.grid", "H32")); + auto tgt_grid = get_grid(args.getString("target.grid", "H32")); auto create_functionspace = [&](Mesh& mesh, int halo, std::string type) -> FunctionSpace { if (type.empty()) { @@ -199,13 +204,13 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { return functionspace::CellColumns(mesh, option::halo(halo)); } else if (type == "NodeColumns") { - return functionspace::NodeColumns(mesh, option::halo(halo)); + return functionspace::NodeColumns(mesh, option::halo(std::max(1,halo))); } ATLAS_THROW_EXCEPTION("FunctionSpace " << type << " is not recognized."); }; timers.target_setup.start(); - auto tgt_mesh = Mesh{tgt_grid}; + auto tgt_mesh = Mesh{tgt_grid, grid::Partitioner(args.getString("target.partitioner", "regular_bands"))}; auto tgt_functionspace = create_functionspace(tgt_mesh, args.getLong("target.halo", 0), args.getString("target.functionspace", "")); auto tgt_field = tgt_functionspace.createField(); @@ -213,8 +218,8 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { timers.source_setup.start(); auto src_meshgenerator = - MeshGenerator{src_grid.meshgenerator() | option::halo(2) | util::Config("pole_elements", "pentagons")}; - auto src_partitioner = grid::MatchingPartitioner{tgt_mesh}; + MeshGenerator{src_grid.meshgenerator() | option::halo(2) | util::Config("pole_elements", "")}; + auto src_partitioner = grid::MatchingPartitioner{tgt_mesh, util::Config("partitioner",args.getString("source.partitioner", "spherical-polygon"))}; auto src_mesh = src_meshgenerator.generate(src_grid, src_partitioner); auto src_functionspace = create_functionspace(src_mesh, args.getLong("source.halo", 2), args.getString("source.functionspace", "")); @@ -230,14 +235,14 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { for (idx_t n = 0; n < lonlat.shape(0); ++n) { src_view(n) = f(PointLonLat{lonlat(n, LON), lonlat(n, LAT)}); } - src_field.set_dirty(false); + src_field.set_dirty(true); timers.initial_condition.start(); } - timers.interpolation_setup.start(); auto interpolation = Interpolation(option::type("conservative-spherical-polygon") | args, src_functionspace, tgt_functionspace); + Log::info() << interpolation << std::endl; timers.interpolation_setup.stop(); @@ -251,11 +256,13 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { using Statistics = interpolation::method::ConservativeSphericalPolygonInterpolation::Statistics; Statistics stats(metadata); if (args.getBool("statistics.accuracy", false)) { - stats.accuracy(interpolation, tgt_field, get_init(args)); + metadata.set( stats.accuracy(interpolation, tgt_field, get_init(args) ) ); } if (args.getBool("statistics.conservation", false)) { // compute difference field src_conservation_field = stats.diff(interpolation, src_field, tgt_field); + src_conservation_field.set_dirty(true); + src_conservation_field.haloExchange(); } } diff --git a/src/tests/interpolation/test_interpolation_conservative.cc b/src/tests/interpolation/test_interpolation_conservative.cc index 92b426ff0..5cfde0448 100644 --- a/src/tests/interpolation/test_interpolation_conservative.cc +++ b/src/tests/interpolation/test_interpolation_conservative.cc @@ -25,6 +25,7 @@ #include "atlas/meshgenerator.h" #include "atlas/option.h" #include "atlas/util/Config.h" +#include "atlas/util/function/VortexRollup.h" #include "tests/AtlasTestEnvironment.h" @@ -36,17 +37,26 @@ using ConservativeMethod = interpolation::method::ConservativeSphericalPolygonIn using Statistics = ConservativeMethod::Statistics; void do_remapping_test(Grid src_grid, Grid tgt_grid, std::function func, - Statistics& remap_stat_1, Statistics& remap_stat_2) { + Statistics& remap_stat_1, Statistics& remap_stat_2, bool src_cell_data, bool tgt_cell_data) { + std::string src_data_type = (src_cell_data ? "CellColumns(" : "NodeColumns("); + std::string tgt_data_type = (tgt_cell_data ? "CellColumns(" : "NodeColumns("); + Log::info() << "+-----------------------\n"; + Log::info() << src_data_type << src_grid.name() << ") --> " << tgt_data_type << tgt_grid.name() <<")\n"; + Log::info() << "+-----------------------\n"; Log::info().indent(); + // setup conservative remap: compute weights, polygon intersection, etc util::Config config("type", "conservative-spherical-polygon"); config.set("order", 1); config.set("validate", true); config.set("statistics.intersection", true); config.set("statistics.conservation", true); + config.set("src_cell_data", src_cell_data); + config.set("tgt_cell_data", tgt_cell_data); auto conservative_interpolation = Interpolation(config, src_grid, tgt_grid); Log::info() << conservative_interpolation << std::endl; + Log::info() << std::endl; // create source field from analytic function "func" const auto& src_fs = conservative_interpolation.source(); @@ -85,6 +95,7 @@ void do_remapping_test(Grid src_grid, Grid tgt_grid, std::function 1st order constructing new matrix"); @@ -95,14 +106,16 @@ void do_remapping_test(Grid src_grid, Grid tgt_grid, std::function 1st order matrix-free"); cfg.set("matrix_free", true); cfg.set("order", 1); auto interpolation = Interpolation(cfg, src_grid, tgt_grid, cache); Log::info() << interpolation << std::endl; interpolation.execute(src_field, tgt_field); + Log::info() << std::endl; } auto cache_2 = interpolation::Cache{}; { @@ -113,14 +126,16 @@ void do_remapping_test(Grid src_grid, Grid tgt_grid, std::function 2nd order matrix-free"); cfg.set("matrix_free", true); cfg.set("order", 2); auto interpolation = Interpolation(cfg, src_grid, tgt_grid, cache); Log::info() << interpolation << std::endl; interpolation.execute(src_field, tgt_field); + Log::info() << std::endl; } { ATLAS_TRACE("cached -> 2nd order using cached matrix"); @@ -129,10 +144,10 @@ void do_remapping_test(Grid src_grid, Grid tgt_grid, std::function tol) { - auto improvement = [](double& e, double& r) { return 100. * (r - e) / r; }; + //auto improvement = [](double& e, double& r) { return 100. * (r - e) / r; }; + auto improvement = [](double& e, double& r) { return r - e; }; double err; // check polygon intersections - err = remap_stat_1.errors[Statistics::Errors::GEO_DIFF]; - Log::info() << "Polygon area computation improvement: " << improvement(err, tol[0]) << " %" << std::endl; + err = remap_stat_1.errors[Statistics::Errors::SRCTGT_INTERSECTPLG_DIFF]; + Log::info() << "Polygon area computation (new < ref) = (" << err << " < " << tol[0] << ")" << std::endl; EXPECT(err < tol[0]); - err = remap_stat_1.errors[Statistics::Errors::GEO_L1]; - Log::info() << "Polygon intersection improvement : " << improvement(err, tol[1]) << " %" << std::endl; + err = remap_stat_1.errors[Statistics::Errors::TGT_INTERSECTPLG_L1]; + Log::info() << "Polygon intersection (new < ref) = (" << err << " < " << tol[1] << ")" << std::endl; EXPECT(err < tol[1]); // check remap accuracy - err = remap_stat_1.errors[Statistics::Errors::REMAP_L2]; - Log::info() << "1st order accuracy improvement : " << improvement(err, tol[2]) << " %" << std::endl; + err = std::abs(remap_stat_1.errors[Statistics::Errors::REMAP_L2]); + Log::info() << "1st order accuracy (new < ref) = (" << err << " < " << tol[2] << ")" << std::endl; EXPECT(err < tol[2]); - err = remap_stat_2.errors[Statistics::Errors::REMAP_L2]; - Log::info() << "2nd order accuracy improvement : " << improvement(err, tol[3]) << " %" << std::endl; + err = std::abs(remap_stat_2.errors[Statistics::Errors::REMAP_L2]); + Log::info() << "2nd order accuracy (new < ref) = (" << err << " < " << tol[3] << ")" << std::endl; EXPECT(err < tol[3]); // check mass conservation - err = remap_stat_1.errors[Statistics::Errors::REMAP_CONS]; - Log::info() << "1st order conservation improvement : " << improvement(err, tol[4]) << " %" << std::endl; + err = std::abs(remap_stat_1.errors[Statistics::Errors::REMAP_CONS]); + Log::info() << "1st order conservation (new < ref) = (" << err << " < " << tol[4] << ")" << std::endl; EXPECT(err < tol[4]); - err = remap_stat_2.errors[Statistics::Errors::REMAP_CONS]; - Log::info() << "2nd order conservation improvement : " << improvement(err, tol[5]) << " %" << std::endl - << std::endl; + err = std::abs(remap_stat_2.errors[Statistics::Errors::REMAP_CONS]); + Log::info() << "2nd order conservation (new < ref) = (" << err << " < " << tol[5] << ")" << std::endl; EXPECT(err < tol[5]); Log::info().unindent(); } CASE("test_interpolation_conservative") { -#if 1 SECTION("analytic constfunc") { auto func = [](const PointLonLat& p) { return 1.; }; Statistics remap_stat_1; Statistics remap_stat_2; - do_remapping_test(Grid("H47"), Grid("H48"), func, remap_stat_1, remap_stat_2); - check(remap_stat_1, remap_stat_2, {1.e-13, 5.e-8, 2.9e-6, 2.9e-6, 5.5e-5, 5.5e-5}); + bool src_cell_data = true; + bool tgt_cell_data = true; + do_remapping_test(Grid("O128"), Grid("H48"), func, remap_stat_1, remap_stat_2, src_cell_data, tgt_cell_data); + check(remap_stat_1, remap_stat_2, {1.0e-13, 1.0e-13, 1.0e-13, 1.0e-13, 1.0e-13, 1.0e-13}); } - SECTION("analytic Y_2^2 as in Jones(1998)") { + SECTION("vortex_rollup") { auto func = [](const PointLonLat& p) { - double cos = std::cos(0.025 * p[0]); - return 2. + cos * cos * std::cos(2 * 0.025 * p[1]); + return util::function::vortex_rollup(p[0], p[1], 0.5); }; Statistics remap_stat_1; Statistics remap_stat_2; - do_remapping_test(Grid("H47"), Grid("H48"), func, remap_stat_1, remap_stat_2); - check(remap_stat_1, remap_stat_2, {1.e-13, 5.e-8, 4.8e-4, 1.1e-4, 8.9e-5, 1.1e-4}); + + bool src_cell_data = true; + bool tgt_cell_data = true; + do_remapping_test(Grid("O32"), Grid("H24"), func, remap_stat_1, remap_stat_2, src_cell_data, tgt_cell_data); + check(remap_stat_1, remap_stat_2, {1.0e-13, 1.0e-12, 3.0e-3, 6.0e-4, 1.0e-15, 1.0e-8}); + + src_cell_data = true; + tgt_cell_data = false; + do_remapping_test(Grid("O32"), Grid("H24"), func, remap_stat_1, remap_stat_2, src_cell_data, tgt_cell_data); + check(remap_stat_1, remap_stat_2, {1.0e-13, 1.0e-12, 3.0e-3, 6.0e-4, 1.0e-15, 1.0e-8}); + + src_cell_data = false; + tgt_cell_data = true; + do_remapping_test(Grid("O32"), Grid("H24"), func, remap_stat_1, remap_stat_2, src_cell_data, tgt_cell_data); + check(remap_stat_1, remap_stat_2, {1.0e-13, 1.0e-12, 3.0e-3, 6.0e-4, 1.0e-15, 1.0e-8}); + + src_cell_data = false; + tgt_cell_data = false; + do_remapping_test(Grid("O32"), Grid("H24"), func, remap_stat_1, remap_stat_2, src_cell_data, tgt_cell_data); + check(remap_stat_1, remap_stat_2, {1.0e-12, 1.0e-12, 3.0e-3, 6.0e-4, 1.0e-15, 1.0e-8}); } -#endif } } // namespace test From 5233ae580553ed12b29d8680861dd164998baf2e Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Wed, 26 Apr 2023 13:51:27 +0200 Subject: [PATCH 3/3] Configurable eps in compare_pointxyz via ATLAS_COMPAREPOINTXYZ_EPS_FACTOR env var New default 1e8*eps instead of 1e4*eps. For previous behaviour: ATLAS_COMPAREPOINTXYZ_EPS_FACTOR=1e4 For high resolution we want to have it to 1e4, for low resolution 1e8 The detection mechanism of source points needs to be replaced with something deterministic that does not rely on comparing points. --- ...nservativeSphericalPolygonInterpolation.cc | 52 ++++++++++--------- 1 file changed, 27 insertions(+), 25 deletions(-) diff --git a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc index 863591212..d1ea5addf 100644 --- a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc +++ b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc @@ -896,30 +896,6 @@ void ConservativeSphericalPolygonInterpolation::do_setup(const FunctionSpace& sr } } -namespace { -// needed for intersect_polygons only, merely for detecting duplicate points -struct ComparePointXYZ { - bool operator()(const PointXYZ& f, const PointXYZ& s) const { - // eps = ConvexSphericalPolygon::EPS which is the threshold when two points are "same" - double eps = 1e4 * std::numeric_limits::epsilon(); - if (f[0] < s[0] - eps) { - return true; - } - else if (std::abs(f[0] - s[0]) < eps) { - if (f[1] < s[1] - eps) { - return true; - } - else if (std::abs(f[1] - s[1]) < eps) { - if (f[2] < s[2] - eps) { - return true; - } - } - } - return false; - } -}; -} // namespace - void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolygonArray& src_csp, const CSPolygonArray& tgt_csp) { ATLAS_TRACE(); @@ -947,7 +923,33 @@ void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolyg StopWatch stopwatch_polygon_intersections; stopwatch_src_already_in.start(); - std::set src_cent; + + + // needed for intersect_polygons only, merely for detecting duplicate points + // Treshold at which points are considered same + double compare_pointxyz_eps = 1.e8 * std::numeric_limits::epsilon(); + if (::getenv("ATLAS_COMPAREPOINTXYZ_EPS_FACTOR")) { + compare_pointxyz_eps = std::atof(::getenv("ATLAS_COMPAREPOINTXYZ_EPS_FACTOR")) * std::numeric_limits::epsilon(); + } + + auto compare_pointxyz = [eps=compare_pointxyz_eps] (const PointXYZ& f, const PointXYZ& s) -> bool { + if (f[0] < s[0] - eps) { + return true; + } + else if (std::abs(f[0] - s[0]) < eps) { + if (f[1] < s[1] - eps) { + return true; + } + else if (std::abs(f[1] - s[1]) < eps) { + if (f[2] < s[2] - eps) { + return true; + } + } + } + return false; + }; + + std::set src_cent(compare_pointxyz); auto src_already_in = [&](const PointXYZ& halo_cent) { if (src_cent.find(halo_cent) == src_cent.end()) { atlas_omp_critical{