Skip to content

Commit 67ea621

Browse files
wdeconinckwdeconinck
authored
Merge pull request ecmwf#331 from ecmwf/feature/interpolation_from_lam
Increase support for parallel interpolation from LAM grids
2 parents 4b42f2a + 53c0c40 commit 67ea621

10 files changed

Lines changed: 496 additions & 76 deletions

File tree

‎src/atlas/grid/StructuredPartitionPolygon.cc‎

Lines changed: 55 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -39,15 +39,16 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect
3939
throw_Exception("halo must zero or matching the one of the StructuredColumns", Here());
4040
}
4141

42-
bool north_pole_included = 90. - grid.y(0) == 0.;
43-
bool south_pole_included = 90. + grid.y(grid.ny() - 1) == 0;
42+
const int y_numbering = (grid.y().front() < grid.y().back()) ? +1 : -1;
43+
bool jfirst_at_pole = y_numbering * 90. + grid.y(0) == 0.;
44+
bool jlast_at_pole = y_numbering * 90. - grid.y(grid.ny() - 1) == 0;
4445

4546
auto compute_j = [&](const idx_t j) {
4647
if (j < 0) {
47-
return -j - 1 + north_pole_included;
48+
return -j - 1 + jfirst_at_pole;
4849
}
4950
else if (j >= grid.ny()) {
50-
return 2 * grid.ny() - j - 1 - south_pole_included;
51+
return 2 * grid.ny() - j - 1 - jlast_at_pole;
5152
}
5253
return j;
5354
};
@@ -107,14 +108,17 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect
107108
};
108109

109110
double ymax, ymin, xmax, xmin;
111+
if (fs.j_begin() >= fs.j_end()) {
112+
return;
113+
}
110114

111115
// Top
112116
// Top left point
113117
{
114118
idx_t j = _halo ? fs.j_begin_halo() : fs.j_begin();
115119
idx_t i = _halo ? fs.i_begin_halo(j) : fs.i_begin(j);
116120
if (j == 0 && _halo == 0) {
117-
p[YY] = dom.ymax();
121+
p[YY] = y_numbering < 0 ? dom.ymax() : dom.ymin();
118122
}
119123
else {
120124
p[YY] = 0.5 * (compute_y(j - 1) + compute_y(j));
@@ -133,7 +137,7 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect
133137
idx_t j = _halo ? fs.j_begin_halo() : fs.j_begin();
134138
idx_t i = _halo ? fs.i_end_halo(j) - 1 : fs.i_end(j) - 1;
135139
if (j == 0 && _halo == 0) {
136-
p[YY] = dom.ymax();
140+
p[YY] = y_numbering < 0 ? dom.ymax() : dom.ymin();
137141
}
138142
else {
139143
p[YY] = 0.5 * (compute_y(j - 1) + compute_y(j));
@@ -194,7 +198,7 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect
194198
idx_t i = _halo ? fs.i_end_halo(j) - 1 : fs.i_end(j) - 1;
195199

196200
if (j == grid.ny() - 1 && _halo == 0) {
197-
p[YY] = dom.ymin();
201+
p[YY] = y_numbering < 0 ? dom.ymin() : dom.ymax();
198202
}
199203
else {
200204
p[YY] = 0.5 * (compute_y(j) + compute_y(j + 1));
@@ -243,7 +247,7 @@ void compute(const functionspace::FunctionSpaceImpl& _fs, idx_t _halo, std::vect
243247
idx_t j = _halo ? fs.j_end_halo() - 1 : fs.j_end() - 1;
244248
idx_t i = _halo ? fs.i_begin_halo(j) : fs.i_begin(j);
245249
if (j == grid.ny() - 1 && _halo == 0) {
246-
p[YY] = dom.ymin();
250+
p[YY] = y_numbering < 0 ? dom.ymin() : dom.ymax();
247251
}
248252
else {
249253
p[YY] = 0.5 * (compute_y(j) + compute_y(j + 1));
@@ -325,7 +329,7 @@ StructuredPartitionPolygon::StructuredPartitionPolygon(const functionspace::Func
325329
compute(fs, halo, points_, inner_bounding_box_);
326330
auto min = Point2(std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
327331
auto max = Point2(std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest());
328-
for (size_t i = 0; i < inner_bounding_box_.size() - 1; ++i) {
332+
for (int i = 0; i < static_cast<int>(inner_bounding_box_.size()) - 1; ++i) {
329333
min = Point2::componentsMin(min, inner_bounding_box_[i]);
330334
max = Point2::componentsMax(max, inner_bounding_box_[i]);
331335
}
@@ -390,45 +394,52 @@ void StructuredPartitionPolygon::outputPythonScript(const eckit::PathName& filep
390394
"\n" "ax = fig.add_subplot(111,aspect='equal')"
391395
"\n";
392396
}
393-
f << "\n" "verts_" << r << " = [";
394-
for ( size_t i=0; i<points.size(); ++i ) {
395-
f << "\n (" << points[i][XX] << ", " << points[i][YY] << "), ";
396-
}
397-
f << "\n]"
398-
"\n"
399-
"\n" "codes_" << r << " = [Path.MOVETO]"
400-
"\n" "codes_" << r << ".extend([Path.LINETO] * " << ( points.size() - 2 ) << ")"
401-
"\n" "codes_" << r << ".extend([Path.CLOSEPOLY])"
402-
"\n"
403-
"\n" "count_" << r << " = " << count <<
404-
"\n" "count_all_" << r << " = " << count_all <<
405-
"\n";
406-
if ( plot_nodes ) {
407-
f << "\n" "x_" << r << " = [";
408-
for ( idx_t i = 0; i < count; ++i ) {
409-
f << xy( i, XX ) << ", ";
397+
if (points.size() > 0) {
398+
f << "\n" "verts_" << r << " = [";
399+
for ( size_t i=0; i<points.size(); ++i ) {
400+
f << "\n (" << points[i][XX] << ", " << points[i][YY] << "), ";
410401
}
411-
f << "]"
412-
"\n" "y_" << r << " = [";
413-
for ( idx_t i = 0; i < count; ++i ) {
414-
f << xy( i, YY ) << ", ";
402+
f << "\n]"
403+
"\n"
404+
"\n" "codes_" << r << " = [Path.MOVETO]"
405+
"\n" "codes_" << r << ".extend([Path.LINETO] * " << ( points.size() - 2 ) << ")"
406+
"\n" "codes_" << r << ".extend([Path.CLOSEPOLY])"
407+
"\n"
408+
"\n" "count_" << r << " = " << count <<
409+
"\n" "count_all_" << r << " = " << count_all <<
410+
"\n";
411+
if ( plot_nodes ) {
412+
f << "\n" "x_" << r << " = [";
413+
for ( idx_t i = 0; i < count; ++i ) {
414+
f << xy( i, XX ) << ", ";
415+
}
416+
f << "]"
417+
"\n" "y_" << r << " = [";
418+
for ( idx_t i = 0; i < count; ++i ) {
419+
f << xy( i, YY ) << ", ";
420+
}
421+
f << "]";
415422
}
416-
f << "]";
417-
}
418-
f << "\n"
419-
"\n" "c = next(colours)"
420-
"\n" "ax.add_patch(patches.PathPatch(Path(verts_" << r << ", codes_" << r << "), edgecolor=c, facecolor=c, alpha=0.3, lw=1))";
421-
if ( plot_nodes ) {
422-
f << "\n" "if plot_nodes:"
423-
"\n" " ax.scatter(x_" << r << ", y_" << r << ", color=c, marker='o')";
423+
f << "\n"
424+
"\n" "c = next(colours)"
425+
"\n" "ax.add_patch(patches.PathPatch(Path(verts_" << r << ", codes_" << r << "), edgecolor=c, facecolor=c, alpha=0.3, lw=1))";
426+
if ( plot_nodes ) {
427+
f << "\n" "if plot_nodes:"
428+
"\n" " ax.scatter(x_" << r << ", y_" << r << ", color=c, marker='o')";
429+
}
430+
f << "\n";
424431
}
425-
f << "\n";
426432
if ( mpi_rank == mpi_size - 1 ) {
427-
f << "\n" "ax.set_xlim( " << xmin << "-5, " << xmax << "+5)"
428-
"\n" "ax.set_ylim(-90-5, 90+5)"
429-
"\n" "ax.set_xticks([0,45,90,135,180,225,270,315,360])"
430-
"\n" "ax.set_yticks([-90,-45,0,45,90])"
431-
"\n" "plt.grid()"
433+
if (fs_.projection().units() == "degrees") {
434+
f << "\n" "ax.set_xlim( " << xmin << "-5, " << xmax << "+5)"
435+
"\n" "ax.set_ylim(-90-5, 90+5)"
436+
"\n" "ax.set_xticks([0,45,90,135,180,225,270,315,360])"
437+
"\n" "ax.set_yticks([-90,-45,0,45,90])";
438+
}
439+
else {
440+
f << "\n" "ax.autoscale()";
441+
}
442+
f << "\n" "plt.grid()"
432443
"\n" "plt.show()";
433444
}
434445
}
@@ -463,7 +474,6 @@ void StructuredPartitionPolygon::allGather(util::PartitionPolygons& polygons_) c
463474
mypolygon.push_back(p[XX]);
464475
mypolygon.push_back(p[YY]);
465476
}
466-
ATLAS_ASSERT(mypolygon.size() >= 4);
467477

468478
eckit::mpi::Buffer<double> recv_polygons(mpi_size);
469479

‎src/atlas/grid/detail/partitioner/MatchingFunctionSpacePartitionerLonLatPolygon.cc‎

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -49,9 +49,10 @@ void MatchingFunctionSpacePartitionerLonLatPolygon::partition(const Grid& grid,
4949
omp::fill(part, part + grid.size(), 0);
5050
}
5151
else {
52-
const auto& p = partitioned_.polygon();
52+
const auto& poly = partitioned_.polygon();
53+
const auto& proj = partitioned_.projection();
5354

54-
util::PolygonXY poly{p};
55+
util::PolygonXY poly_xy{poly};
5556
{
5657
ATLAS_TRACE("point-in-polygon check for entire grid (" + std::to_string(grid.size()) + " points)");
5758
size_t num_threads = atlas_omp_get_max_threads();
@@ -64,15 +65,17 @@ void MatchingFunctionSpacePartitionerLonLatPolygon::partition(const Grid& grid,
6465
for( size_t chunk=0; chunk < chunks; ++chunk) {
6566
const size_t begin = chunk * size_t(grid.size()) / chunks;
6667
const size_t end = (chunk + 1) * size_t(grid.size()) / chunks;
67-
auto it = grid.xy().begin() + chunk * grid.size() / chunks;
68+
auto lonlat_it = grid.lonlat().begin() + chunk * grid.size() / chunks;
69+
PointXY xy;
6870
for (size_t n = begin; n < end; ++n) {
69-
if (poly.contains(*it)) {
71+
xy = proj.xy(*lonlat_it);
72+
if (poly_xy.contains(xy)) {
7073
part[n] = mpi_rank;
7174
}
7275
else {
7376
part[n] = -1;
7477
}
75-
++it;
78+
++lonlat_it;
7679
}
7780
}
7881
}

‎src/atlas/interpolation/method/structured/RegionalLinear2D.cc‎

Lines changed: 1 addition & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -45,9 +45,6 @@ void RegionalLinear2D::do_setup(const FunctionSpace& source,
4545
source_ = source;
4646
target_ = target;
4747

48-
if (target_.size() == 0) {
49-
return;
50-
}
5148
ASSERT(source_.type() == "StructuredColumns");
5249

5350
// Get grid parameters
@@ -160,7 +157,7 @@ void RegionalLinear2D::do_setup(const FunctionSpace& source,
160157
// Buffer size
161158
targetRecvSize_ = targetRecvPointsList.size();
162159

163-
if (targetRecvSize_ > 0) {
160+
{
164161
// RecvDispls
165162
targetRecvDispls_.push_back(0);
166163
for (size_t jt = 0; jt < comm_.size()-1; ++jt) {
@@ -413,10 +410,6 @@ void RegionalLinear2D::do_execute(const Field& sourceField, Field& targetField,
413410
Metadata&) const {
414411
ATLAS_TRACE("atlas::interpolation::method::RegionalLinear2D::do_execute()");
415412

416-
if (targetField.size() == 0) {
417-
return;
418-
}
419-
420413
// Check number of levels
421414
ASSERT(sourceField.levels() == targetField.levels());
422415
const size_t nz = sourceField.levels() > 0 ? sourceField.levels() : 1;

‎src/atlas/mesh/PartitionPolygon.cc‎

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -238,11 +238,16 @@ void PartitionPolygon::outputPythonScript(const eckit::PathName& filepath, const
238238
return xticks_ss.str();
239239
};
240240
// clang-format off
241-
f << "\n" "ax.set_xlim( " << xmin << "-5, " << xmax << "+5)"
242-
"\n" "ax.set_ylim(-90-5, 90+5)"
243-
"\n" "ax.set_xticks(" << xticks() << ")"
244-
"\n" "ax.set_yticks([-90,-45,0,45,90])"
245-
"\n" "plt.grid()"
241+
if (mesh_.projection().units() == "degrees") {
242+
f << "\n" "ax.set_xlim( " << xmin << "-5, " << xmax << "+5)"
243+
"\n" "ax.set_ylim(-90-5, 90+5)"
244+
"\n" "ax.set_xticks(" << xticks() << ")"
245+
"\n" "ax.set_yticks([-90,-45,0,45,90])";
246+
}
247+
else {
248+
f << "\n" "ax.autoscale()";
249+
}
250+
f << "\n" "plt.grid()"
246251
"\n" "plt.show()";
247252
// clang-format on
248253
}

‎src/atlas/meshgenerator/detail/RegularMeshGenerator.cc‎

Lines changed: 14 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -178,6 +178,7 @@ void RegularMeshGenerator::generate_mesh(const RegularGrid& rg, const grid::Dist
178178
int nparts = options.get<size_t>("nb_parts");
179179
int nx = rg.nx();
180180
int ny = rg.ny();
181+
const int y_numbering = (rg.y().front() < rg.y().back()) ? +1 : -1;
181182

182183
bool periodic_x = options.get<bool>("periodic_x") or rg.periodic();
183184
bool periodic_y = options.get<bool>("periodic_y");
@@ -529,10 +530,19 @@ void RegularMeshGenerator::generate_mesh(const RegularGrid& rg, const grid::Dist
529530
if (!is_ghost_SR[ii]) {
530531
if ((ix_min + ix < nx - 1 || periodic_x) && (iy_min + iy < ny - 1 || periodic_y)) {
531532
// define cell corners (local indices)
532-
quad_nodes[0] = local_idx_SR[ii];
533-
quad_nodes[3] = local_idx_SR[iy * nxl + ix + 1]; // point to the right
534-
quad_nodes[2] = local_idx_SR[(iy + 1) * nxl + ix + 1]; // point above right
535-
quad_nodes[1] = local_idx_SR[(iy + 1) * nxl + ix]; // point above
533+
if (y_numbering < 0) {
534+
quad_nodes[0] = local_idx_SR[ii];
535+
quad_nodes[3] = local_idx_SR[iy * nxl + ix + 1]; // point to the right
536+
quad_nodes[2] = local_idx_SR[(iy + 1) * nxl + ix + 1]; // point above right
537+
quad_nodes[1] = local_idx_SR[(iy + 1) * nxl + ix]; // point above
538+
}
539+
else {
540+
quad_nodes[0] = local_idx_SR[ii];
541+
quad_nodes[1] = local_idx_SR[iy * nxl + ix + 1]; // point to the right
542+
quad_nodes[2] = local_idx_SR[(iy + 1) * nxl + ix + 1]; // point above right
543+
quad_nodes[3] = local_idx_SR[(iy + 1) * nxl + ix]; // point above
544+
545+
}
536546
node_connectivity.set(jcell, quad_nodes);
537547
cells_part(jcell) = mypart;
538548
#if DEBUG_OUTPUT_DETAIL

‎src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -246,7 +246,6 @@ void StructuredMeshGenerator::generate(const Grid& grid, const grid::Distributio
246246
// show distribution
247247
#if DEBUG_OUTPUT
248248
int inode = 0;
249-
const atlas::vector<int>& parts = distribution;
250249
Log::info() << "Partition : " << std::endl;
251250
for (size_t ilat = 0; ilat < rg.ny(); ilat++) {
252251
for (size_t ilon = 0; ilon < rg.nx(ilat); ilon++) {
@@ -360,6 +359,12 @@ Find min and max latitudes used by this part.
360359
}
361360
lat_north = *std::min_element(thread_reduce_lat_north.begin(), thread_reduce_lat_north.end());
362361
lat_south = *std::max_element(thread_reduce_lat_south.begin(), thread_reduce_lat_south.end());
362+
if (lat_north == std::numeric_limits<idx_t>::max()) {
363+
lat_north = -1;
364+
}
365+
if (lat_south == std::numeric_limits<idx_t>::min()) {
366+
lat_south = -1;
367+
}
363368
}
364369
}
365370
}

‎src/atlas/util/Polygon.cc‎

Lines changed: 10 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -258,15 +258,17 @@ std::ostream& operator<<(std::ostream& out, const PolygonCoordinates& pc) {
258258

259259
Polygon::edge_set_t ExplicitPartitionPolygon::compute_edges(idx_t points_size) {
260260
util::Polygon::edge_set_t edges;
261-
auto add_edge = [&](idx_t p1, idx_t p2) {
262-
util::Polygon::edge_t edge = {p1, p2};
263-
edges.insert(edge);
264-
// Log::info() << edge.first << " " << edge.second << std::endl;
265-
};
266-
for (idx_t p = 0; p < points_size - 2; ++p) {
267-
add_edge(p, p + 1);
261+
if (points_size > 0) {
262+
auto add_edge = [&](idx_t p1, idx_t p2) {
263+
util::Polygon::edge_t edge = {p1, p2};
264+
edges.insert(edge);
265+
// Log::info() << edge.first << " " << edge.second << std::endl;
266+
};
267+
for (idx_t p = 0; p < points_size - 2; ++p) {
268+
add_edge(p, p + 1);
269+
}
270+
add_edge(points_size - 2, 0);
268271
}
269-
add_edge(points_size - 2, 0);
270272
return edges;
271273
}
272274

‎src/tests/AtlasTestEnvironment.h‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -224,6 +224,30 @@ bool approx_eq(const Point3& v1, const Point3& v2, const double& t) {
224224
return approx_eq(v1[0], v2[0], t) && approx_eq(v1[1], v2[1], t) && approx_eq(v1[2], v2[2], t);
225225
}
226226

227+
template<typename P1, typename A1, typename P2, typename A2>
228+
bool approx_eq(const std::vector<P1,A1>& v1, const std::vector<P2,A2>& v2) {
229+
if (v1.size() != v2.size()) {
230+
return false;
231+
}
232+
bool _approx_eq = true;
233+
for (size_t j=0; j<v1.size(); ++j) {
234+
_approx_eq &= approx_eq(v1[j], v2[j]);
235+
}
236+
return _approx_eq;
237+
}
238+
239+
template<typename P1, typename A1, typename P2, typename A2>
240+
bool approx_eq(const std::vector<P1,A1>& v1, const std::vector<P2,A2>& v2, const double& t) {
241+
if (v1.size() != v2.size()) {
242+
return false;
243+
}
244+
bool _approx_eq = true;
245+
for (size_t j=0; j<v1.size(); ++j) {
246+
_approx_eq &= approx_eq(v1[j], v2[j], t);
247+
}
248+
return _approx_eq;
249+
}
250+
227251
template <typename T1, typename T2>
228252
std::string expect_message(const std::string& condition, const T1& lhs, const T2& rhs, const eckit::CodeLocation& loc) {
229253
std::stringstream msg;

‎src/tests/interpolation/CMakeLists.txt‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -128,6 +128,14 @@ ecbuild_add_test( TARGET atlas_test_interpolation_structured2D_to_unstructured
128128
ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT}
129129
)
130130

131+
ecbuild_add_test( TARGET atlas_test_interpolation_from_lam
132+
SOURCES test_interpolation_from_lam.cc
133+
LIBS atlas
134+
MPI 4
135+
CONDITION eckit_HAVE_MPI AND MPI_SLOTS GREATER_EQUAL 4
136+
ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT}
137+
)
138+
131139
ecbuild_add_test( TARGET atlas_test_interpolation_structured2D_to_points
132140
SOURCES test_interpolation_structured2D_to_points.cc
133141
LIBS atlas

0 commit comments

Comments
 (0)