fixes bug in routine get_cell_lon_lat_bounds

* average with barycentric coordinates weightening could case an error
  for regular source cells, if one of the edges is on a pole

src/core/interp_method_avg.c

@@ -571,10 +571,10 @@ static void get_cell_lon_lat_bounds(
   int lat_neigh_offset =
     (fabs(coords[0][2] - coords[1][2]) <
      fabs(coords[0][2] - coords[3][2]))?1:3;
-  // get the index of the coord which is closes to the pole
-  // (the one with the lowest absolut z coordinate)
+  // get the index of the coord which is closest to the equator
+  // (the one with the lower absolut z coordinate)
   int closest_to_equator_idx =
-    (fabs(coords[0][2]) < fabs(coords[1][2]))?0:2;
+    (fabs(coords[0][2]) < fabs(coords[2][2]))?0:2;
   int upper_edge_idx = ((coords[0][2]) > coords[2][2])?0:2;
   int lower_edge_idx = upper_edge_idx^2;
   int upper_edge_is_closer_to_pole =
@@ -589,10 +589,10 @@ static void get_cell_lon_lat_bounds(
   XYZtoLL(
     coords[closest_to_equator_idx], &lon[closest_to_equator_edge_ordering],
     &lat[upper_edge_is_closer_to_pole]);
+  double dummy;
   XYZtoLL(
     coords[(closest_to_equator_idx + lat_neigh_offset)%4],
-    &lon[closest_to_equator_edge_ordering^1],
-    &lat[upper_edge_is_closer_to_pole]);
+    &lon[closest_to_equator_edge_ordering^1], &dummy);
   lat[!upper_edge_is_closer_to_pole] =
     M_PI_2 - acos(coords[closest_to_equator_idx ^ 2][2]);
 
@@ -602,6 +602,19 @@ static void get_cell_lon_lat_bounds(
   reorder[2+(cell_ordering^1)] = (upper_edge_idx+lat_neigh_offset)%4;
 
   if (lon[1] < lon[0]) lon[1] += 2.0 * M_PI;
+
+  YAC_ASSERT_F(((lon[0] - lon[1]) != 0.0) && ((lat[0] - lat[1]) != 0.0),
+      "ERROR(get_cell_lon_lat_bounds): internal error\n"
+      "lon[0] = %e lon[1] = %e lat[0] = %e lat[1] = %e\n"
+      "cell coords: ((% .4e; % .4e; % .4e),\n"
+      "              (% .4e; % .4e; % .4e),\n"
+      "              (% .4e; % .4e; % .4e),\n"
+      "              (% .4e; % .4e; % .4e)",
+      lon[0], lon[1], lat[0], lat[1],
+      coords[0][0], coords[0][1], coords[0][2],
+      coords[1][0], coords[1][1], coords[1][2],
+      coords[2][0], coords[2][1], coords[2][2],
+      coords[3][0], coords[3][1], coords[3][2]);
 }
 
 static void get_point_lon_lat(

tests/core/test_interp_method_avg_parallel.c

@@ -3,6 +3,7 @@
 // SPDX-License-Identifier: BSD-3-Clause
 
 #include <stdlib.h>
+#include <string.h>
 
 #include "tests.h"
 #include "test_common.h"
@@ -1669,6 +1670,210 @@ int main(void) {
     yac_basic_grid_delete(grids[0]);
   }
 
+  {// Test 10
+   // barycentric coordinates weighted (cell touches pole and has duplicated
+   // longitude coordinate at pole)
+   // interpolation with two source processes and a single target process
+
+    // the global source grid is a 1x1 grid:
+    //       02--03--03
+    //       |       |
+    //       01  00  02
+    //       |       |
+    //       00--00--01
+    //
+    //---------------
+    // setup
+    //---------------
+
+    int is_tgt = split_comm_size == 1;
+    double src_coordinates_x[2][2] = {{-0.5,0.5}, {0.5,-0.5}};
+    double src_coordinates_y[2][2] = {{89,90}, {90,89}};
+    size_t const src_num_cells[2] = {1,1};
+    size_t src_local_start[2] = {0,0};
+    size_t src_local_count[2] = {1,1};
+    yac_int src_reorder[2][2][4] =
+      {{{0,1,2,3},{2,3,0,1}},{{1,0,3,2},{3,2,1,0}}};
+    double tgt_coordinates_x[5] = {-0.5,-0.25,0.0,0.25,0.5};
+    double tgt_coordinates_y[5] = {89.00,89.25,89.50,89.75,90.00};
+    size_t const tgt_num_cells[2] = {4,4};
+    size_t tgt_local_start[2] = {0,0};
+    size_t tgt_local_count[2] = {4,4};
+
+    for (size_t i = 0; i < 2 * 2; ++i)
+      (&(src_coordinates_x[0][0]))[i] *= YAC_RAD;
+    for (size_t i = 0; i < 2 * 2; ++i)
+      (&(src_coordinates_y[0][0]))[i] *= YAC_RAD;
+    for (size_t i = 0; i < 5; ++i) tgt_coordinates_x[i] *= YAC_RAD;
+    for (size_t i = 0; i < 5; ++i)
+      (&(tgt_coordinates_y[0]))[i] *= YAC_RAD;
+
+    for (int order_x = 0; order_x < 2; ++order_x) {
+      for (int order_y = 0; order_y < 2; ++order_y) {
+
+        struct yac_basic_grid_data grid_data;
+
+        if (is_tgt) {
+          grid_data =
+            yac_generate_basic_grid_data_reg2d(
+              tgt_coordinates_x, tgt_coordinates_y,
+              tgt_num_cells, tgt_local_start, tgt_local_count, 0);
+        } else {
+          grid_data =
+            yac_generate_basic_grid_data_reg2d(
+              src_coordinates_x[order_x],
+              src_coordinates_y[order_y], src_num_cells,
+              src_local_start, src_local_count, 0);
+        }
+
+        if (!is_tgt) {
+
+          yac_int * vertex_ids =
+            xmalloc(grid_data.num_vertices * sizeof(*vertex_ids));
+          for (size_t i = 0; i < grid_data.num_vertices; ++i)
+            vertex_ids[i] = (yac_int)(src_reorder[order_x][order_y][i]);
+          grid_data.vertex_ids = vertex_ids;
+        }
+
+        struct yac_basic_grid * grids[2] =
+          {yac_basic_grid_new(grid_names[is_tgt], grid_data),
+            yac_basic_grid_empty_new(grid_names[is_tgt^1])};
+
+        if (!is_tgt) {
+
+          yac_coordinate_pointer vertex_coordinates =
+            xmalloc(grid_data.num_vertices * sizeof(*vertex_coordinates));
+          memcpy(
+            vertex_coordinates, grid_data.vertex_coordinates,
+            grid_data.num_vertices * sizeof(*vertex_coordinates));
+          // set x and y coordinate of all pole points to 0.0 in order to
+          // reproduce an error that occured in ICON
+          for (size_t i = 0; i < grid_data.num_vertices; ++i) {
+            if (fabs(1.0 - fabs(vertex_coordinates[i][2])) < 1e-9) {
+              vertex_coordinates[i][0] = 0.0;
+              vertex_coordinates[i][1] = 0.0;
+            }
+          }
+          yac_basic_grid_add_coordinates_nocpy(
+            grids[0], YAC_LOC_CORNER, vertex_coordinates);
+        }
+
+        struct yac_dist_grid_pair * grid_pair =
+          yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
+
+        struct yac_interp_field src_fields[] =
+          {{.location = YAC_LOC_CORNER, .coordinates_idx = 0,
+            .masks_idx = SIZE_MAX}};
+        size_t num_src_fields =
+          sizeof(src_fields) / sizeof(src_fields[0]);
+        struct yac_interp_field tgt_field =
+          {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX,
+            .masks_idx = SIZE_MAX};
+
+        struct yac_interp_grid * interp_grid =
+          yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
+                              num_src_fields, src_fields, tgt_field);
+
+        int partial_coverage = 0;
+        struct interp_method * method_stack[] =
+          {yac_interp_method_avg_new(
+              YAC_INTERP_AVG_BARY, partial_coverage),
+            yac_interp_method_fixed_new(-1.0), NULL};
+
+        struct yac_interp_weights * weights =
+          yac_interp_method_do_search(method_stack, interp_grid);
+
+        for (size_t i = 0; i < num_reorder_types; ++i) {
+
+          struct yac_interpolation * interpolation =
+            yac_interp_weights_get_interpolation(
+              weights, reorder_types[i], 1,
+              YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0);
+
+          // check generated interpolation
+          {
+            double * src_field = NULL;
+            double ** src_fields = &src_field;
+            double * tgt_field = NULL;
+            double * ref_tgt_field = NULL;
+            double ref_global_tgt_field[25];
+            double ref_weights[25][4] =
+              {{1.00, 0.00, 0.00, 0.00},
+                {0.75, 0.25, 0.00, 0.00},
+                {0.50, 0.50, 0.00, 0.00},
+                {0.25, 0.75, 0.00, 0.00},
+                {0.00, 1.00, 0.00, 0.00},
+                {0.75, 0.00, 0.25, 0.00},
+                {0.75, 0.00, 0.00, 0.25},
+                {0.50, 0.25, 0.00, 0.25},
+                {0.25, 0.50, 0.00, 0.25},
+                {0.00, 0.75, 0.00, 0.25},
+                {0.50, 0.00, 0.50, 0.00},
+                {0.50, 0.00, 0.25, 0.25},
+                {0.50, 0.00, 0.00, 0.50},
+                {0.25, 0.25, 0.00, 0.50},
+                {0.00, 0.50, 0.00, 0.50},
+                {0.25, 0.00, 0.75, 0.00},
+                {0.25, 0.00, 0.50, 0.25},
+                {0.25, 0.00, 0.25, 0.50},
+                {0.25, 0.00, 0.00, 0.75},
+                {0.00, 0.25, 0.00, 0.75},
+                {0.00, 0.00, 1.00, 0.00},
+                {0.00, 0.00, 0.75, 0.25},
+                {0.00, 0.00, 0.50, 0.50},
+                {0.00, 0.00, 0.25, 0.75},
+                {0.00, 0.00, 0.00, 1.00}};
+
+            for (size_t j = 0; j < 25; ++j) {
+              ref_global_tgt_field[j] =
+                (double)(src_reorder[order_x][order_y][0] + 1) *
+                ref_weights[j][0] +
+                (double)(src_reorder[order_x][order_y][1] + 1) *
+                ref_weights[j][1] +
+                (double)(src_reorder[order_x][order_y][2] + 1) *
+                ref_weights[j][2] +
+                (double)(src_reorder[order_x][order_y][3] + 1) *
+                ref_weights[j][3];
+              if (ref_global_tgt_field[j] == 0.0)
+                ref_global_tgt_field[j] = -1.0;
+            }
+
+            if (is_tgt) {
+              tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
+              ref_tgt_field = xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
+              for (size_t i = 0; i < grid_data.num_vertices; ++i)
+                ref_tgt_field[i] = ref_global_tgt_field[i];
+            } else {
+              src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
+              for (size_t i = 0; i < grid_data.num_vertices; ++i)
+                src_field[i] = (double)(i) + 1.0;
+            }
+
+            yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
+
+            if (is_tgt)
+              for (size_t i = 0; i < grid_data.num_vertices; ++i)
+                if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-3)
+                  PUT_ERR("wrong interpolation result");
+
+            free(src_field);
+            free(tgt_field);
+            free(ref_tgt_field);
+          }
+
+          yac_interpolation_delete(interpolation);
+        }
+
+        yac_interp_weights_delete(weights);
+        yac_interp_method_delete(method_stack);
+        yac_interp_grid_delete(interp_grid);
+        yac_dist_grid_pair_delete(grid_pair);
+        yac_basic_grid_delete(grids[1]);
+        yac_basic_grid_delete(grids[0]);
+      }
+    }
+  }
+
   MPI_Comm_free(&split_comm);
 
   xt_finalize();