interp_method_spmap.c

Go to the documentation of this file.

// Copyright (c) 2024 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#ifdef HAVE_CONFIG_H
// Get the definition of the 'restrict' keyword.
#include "config.h"
#endif
 
#include <float.h>
#include <string.h>
 
#include "interp_method_internal.h"
#include "interp_method_spmap.h"
#include "ensure_array_size.h"
#include "area.h"
#include "io_utils.h"
#include "yac_mpi_internal.h"
 
static size_t do_search_spmap(struct interp_method * method,
                              struct yac_interp_grid * interp_grid,
                              size_t * tgt_points, size_t count,
                              struct yac_interp_weights * weights,
                              int * interpolation_complete);
static void delete_spmap(struct interp_method * method);
 
static struct interp_method_vtable
  interp_method_spmap_vtable = {
    .do_search = do_search_spmap,
    .delete = delete_spmap
};
 
struct interp_method_spmap {
 
  struct interp_method_vtable * vtable;
  double spread_distance;
  double max_search_distance;
  enum yac_interp_spmap_weight_type weight_type;
  struct yac_spmap_scale_config scale_config;
};
 
static inline int compare_size_t(const void * a, const void * b) {
 
  size_t const * a_ = a, * b_ = b;
 
  return (*a_ > *b_) - (*b_ > *a_);
}
 
static int lists_overlap(
  size_t * list_a, size_t count_a, size_t * list_b, size_t count_b) {
 
  if ((count_a == 0) || (count_b == 0)) return 0;
 
  size_t i = 0, j = 0;
  size_t curr_a = SIZE_MAX, curr_b = list_b[0];
 
  do {
    while ((i < count_a) && (((curr_a = list_a[i++])) < curr_b));
    if (curr_a == curr_b) return 1;
    while ((j < count_b) && (((curr_b = list_b[j++])) < curr_a));
    if (curr_a == curr_b) return 1;
  } while ((i < count_a) || (j < count_b));
 
  return 0;
}
 
static void merge_lists(
  size_t * list, size_t * list_size, size_t * insert, size_t insert_size) {
 
  if (insert_size == 0) return;
 
  size_t new_list_size = *list_size;
  size_t old_list_size = *list_size;
 
  for (size_t i = 0, j = 0; i < insert_size; ++i) {
    size_t curr_insert = insert[i];
    while ((j < old_list_size) && (list[j] < curr_insert)) ++j;
    if ((j >= old_list_size) || (list[j] != curr_insert))
      list[new_list_size++] = curr_insert;
  }
 
  if (new_list_size != old_list_size) {
    qsort(list, new_list_size, sizeof(*list), compare_size_t);
    *list_size = new_list_size;
  }
}
 
static void check_spread_distance(
  yac_const_coordinate_pointer tgt_field_coords, double max_distance,
  size_t tgt_start_point, size_t * tgt_points, size_t * count) {
 
  double const * start_coord = tgt_field_coords[tgt_start_point];
  size_t new_count = 0;
  for (size_t i = 0, old_count = *count; i < old_count; ++i) {
    if (get_vector_angle(
          start_coord, tgt_field_coords[tgt_points[i]]) <= max_distance) {
      if (new_count != i) tgt_points[new_count] = tgt_points[i];
      ++new_count;
    }
  }
 
  *count = new_count;
}
 
static void remove_disconnected_points(
  struct yac_interp_grid * interp_grid, size_t tgt_start_point,
  size_t * from_tgt_points, size_t * to_tgt_points,
  size_t * count, int * flag, size_t * temp_cell_edges,
  size_t ** edges_buffer, size_t * edges_buffer_array_size) {
 
  struct yac_const_basic_grid_data * tgt_grid_data =
    yac_interp_grid_get_basic_grid_data_tgt(interp_grid);
  const_size_t_pointer cell_to_edge = tgt_grid_data->cell_to_edge;
  const_size_t_pointer cell_to_edge_offsets =
    tgt_grid_data->cell_to_edge_offsets;
  const int * num_vertices_per_cell =
    tgt_grid_data->num_vertices_per_cell;
 
  size_t old_count = *count;
  memset(flag, 0, old_count * sizeof(*flag));
 
  for (size_t i = 0; i < old_count; ++i) {
    if (from_tgt_points[i] == tgt_start_point) {
      flag[i] = 1;
      break;
    }
  }
 
  size_t * edges = *edges_buffer;
  size_t num_edges = 0;
  size_t edges_array_size = *edges_buffer_array_size;
 
  const_size_t_pointer curr_edges =
    cell_to_edge + cell_to_edge_offsets[tgt_start_point];
  size_t curr_num_edges = num_vertices_per_cell[tgt_start_point];
 
  ENSURE_ARRAY_SIZE(edges, edges_array_size, curr_num_edges);
  num_edges = curr_num_edges;
  memcpy(edges, curr_edges, num_edges * sizeof(*edges));
  qsort(edges, num_edges, sizeof(*edges), compare_size_t);
 
  int change_flag = 0;
  do {
 
    change_flag = 0;
 
    for (size_t i = 0; i < old_count; ++i) {
 
      if (flag[i]) continue;
 
      const_size_t_pointer curr_edges =
        tgt_grid_data->cell_to_edge + cell_to_edge_offsets[from_tgt_points[i]];
      curr_num_edges = num_vertices_per_cell[from_tgt_points[i]];
      memcpy(
        temp_cell_edges, curr_edges, curr_num_edges * sizeof(*temp_cell_edges));
      qsort(temp_cell_edges, curr_num_edges, sizeof(*edges), compare_size_t);
      ENSURE_ARRAY_SIZE(edges, edges_array_size, num_edges + curr_num_edges);
 
      if (lists_overlap(edges, num_edges, temp_cell_edges, curr_num_edges)) {
        merge_lists(edges, &num_edges, temp_cell_edges, curr_num_edges);
        flag[i] = 1;
        change_flag = 1;
      }
    }
  } while (change_flag);
 
  *edges_buffer = edges;
  *edges_buffer_array_size = edges_array_size;
 
  size_t new_count = 0;
  for (size_t i = 0; i < old_count; ++i)
    if (flag[i]) to_tgt_points[new_count++] = from_tgt_points[i];
 
  *count = new_count;
}
 
static size_t check_tgt_result_points(
  struct yac_interp_grid * interp_grid, double spread_distance,
  size_t num_src_points, size_t const * const tgt_result_points,
  size_t * num_tgt_per_src, size_t * spread_tgt_result_points) {
 
  size_t max_num_tgt_per_src = 0;
  for (size_t i = 0; i < num_src_points; ++i)
    if (num_tgt_per_src[i] > max_num_tgt_per_src)
      max_num_tgt_per_src = num_tgt_per_src[i];
 
  int * flag = xmalloc(max_num_tgt_per_src * sizeof(*flag));
 
  size_t new_offset = 0;
  size_t * cell_edge_buffer = NULL;
  size_t cell_edge_buffer_array_size = 0;
  size_t * edge_buffer = NULL;
  size_t edge_buffer_array_size = 0;
  int max_num_vertice_per_tgt = 0;
  const int * num_vertices_per_tgt =
    yac_interp_grid_get_basic_grid_data_tgt(interp_grid)->
      num_vertices_per_cell;
 
  yac_const_coordinate_pointer tgt_field_coords =
    yac_interp_grid_get_tgt_field_coords(interp_grid);
 
  // for all source points
  for (size_t i = 0, old_offset = 0; i < num_src_points; ++i) {
 
    size_t * old_results = spread_tgt_result_points + old_offset;
    size_t * new_results = spread_tgt_result_points + new_offset;
    old_offset += num_tgt_per_src[i];
 
    // remove all target points, which exceed the spread distance from
    // the original tgt
    check_spread_distance(
      tgt_field_coords, spread_distance, tgt_result_points[i],
      old_results, num_tgt_per_src + i);
 
    // check buffer sizes required by routine remove_disconnected_points
    for (size_t j = 0, curr_num_tgt_per_src = num_tgt_per_src[i];
        j < curr_num_tgt_per_src; ++j)
      if (num_vertices_per_tgt[old_results[j]] > max_num_vertice_per_tgt)
        max_num_vertice_per_tgt = num_vertices_per_tgt[old_results[j]];
 
    ENSURE_ARRAY_SIZE(
      cell_edge_buffer, cell_edge_buffer_array_size,
      max_num_vertice_per_tgt);
 
    // remove all tgts that are not directly connected to the original tgt
    remove_disconnected_points(
      interp_grid, tgt_result_points[i],
      old_results, new_results, num_tgt_per_src + i, flag,
      cell_edge_buffer, &edge_buffer, &edge_buffer_array_size);
 
    new_offset += num_tgt_per_src[i];
  }
 
  free(edge_buffer);
  free(cell_edge_buffer);
  free(flag);
 
  return new_offset;
}
 
static double * compute_weights(
  struct yac_interp_grid * interp_grid,
  enum yac_interp_spmap_weight_type weight_type, size_t num_src_points,
  size_t const * const src_points, size_t const * const num_tgt_per_src,
  size_t total_num_tgt, size_t const * const tgt_result_points) {
 
  double * weights = xmalloc(total_num_tgt * sizeof(*weights));
  YAC_ASSERT(
    (weight_type == YAC_INTERP_SPMAP_AVG) ||
    (weight_type == YAC_INTERP_SPMAP_DIST),
    "ERROR(do_search_spmap): invalid weight_type")
  switch (weight_type) {
    case (YAC_INTERP_SPMAP_AVG): {
      for (size_t i = 0, offset = 0; i < num_src_points; ++i) {
        size_t curr_num_tgt = num_tgt_per_src[i];
        if (curr_num_tgt == 0) continue;
        if (curr_num_tgt > 1) {
          double curr_weight_data = 1.0 / (double)(curr_num_tgt);
          for (size_t j = 0; j < curr_num_tgt; ++j, ++offset) {
            weights[offset] = curr_weight_data;
          }
        } else {
          weights[offset] = 1.0;
          ++offset;
        }
      }
      break;
    }
    default:
    case (YAC_INTERP_SPMAP_DIST): {
 
      yac_const_coordinate_pointer src_field_coords =
        yac_interp_grid_get_src_field_coords(interp_grid, 0);
      yac_const_coordinate_pointer tgt_field_coords =
        yac_interp_grid_get_tgt_field_coords(interp_grid);
 
      for (size_t i = 0, offset = 0; i < num_src_points; ++i) {
 
        size_t curr_num_tgt = num_tgt_per_src[i];
 
        if (curr_num_tgt == 0) continue;
 
        double * curr_weights = weights + offset;
 
        if (curr_num_tgt > 1) {
 
          size_t const * const  curr_result_points =
            tgt_result_points + offset;
          double const * curr_src_coord = src_field_coords[src_points[offset]];
          offset += curr_num_tgt;
 
          int match_flag = 0;
 
          for (size_t j = 0; j < curr_num_tgt; ++j) {
 
            double distance =
              get_vector_angle(
                (double*)curr_src_coord,
                (double*)tgt_field_coords[curr_result_points[j]]);
 
            if (distance < yac_angle_tol) {
              for (size_t k = 0; k < curr_num_tgt; ++k) curr_weights[k] = 0.0;
              curr_weights[j] = 1.0;
              match_flag = 1;
              break;
            }
            curr_weights[j] = 1.0 / distance;
          }
 
          if (!match_flag) {
 
            // compute scaling factor for the weights
            double inv_distance_sum = 0.0;
            for (size_t j = 0; j < curr_num_tgt; ++j)
              inv_distance_sum += curr_weights[j];
            double scale = 1.0 / inv_distance_sum;
 
            for (size_t j = 0; j < curr_num_tgt; ++j) curr_weights[j] *= scale;
          }
        } else {
          *curr_weights = 1.0;
          ++offset;
        }
      }
      break;
    }
  };
 
  return weights;
}
 
static void compute_cell_areas(
  struct yac_const_basic_grid_data * basic_grid_data,
  int * required_cell_areas, double area_scale, char const * type,
  double * cell_areas) {
 
  struct yac_grid_cell grid_cell;
  yac_init_grid_cell(&grid_cell);
 
  size_t num_cells = basic_grid_data->count[YAC_LOC_CELL];
 
  for (size_t i = 0; i < num_cells; ++i) {
    if (required_cell_areas[i]) {
      yac_const_basic_grid_data_get_grid_cell(basic_grid_data, i, &grid_cell);
      double cell_area = yac_huiliers_area(grid_cell) * area_scale;
      YAC_ASSERT_F(
        cell_area > YAC_AREA_TOL,
        "ERROR(get_cell_areas): "
        "area of %s cell (global id %"XT_INT_FMT") is close to zero (%e)",
        type, basic_grid_data->ids[YAC_LOC_CELL][i], cell_area);
      cell_areas[i] = cell_area;
    } else {
      cell_areas[i] = 0.0;
    }
  }
 
  yac_free_grid_cell(&grid_cell);
}
 
static void dist_read_cell_areas(
  char const * filename, char const * varname, MPI_Comm comm,
  int * io_ranks, int io_rank_idx, int num_io_ranks,
  double ** dist_cell_areas, size_t * dist_cell_areas_global_size) {
 
#ifndef YAC_NETCDF_ENABLED
 
  UNUSED(filename);
  UNUSED(varname);
  UNUSED(comm);
  UNUSED(io_ranks);
  UNUSED(io_rank_idx);
  UNUSED(num_io_ranks);
 
  *dist_cell_areas = NULL;
  *dist_cell_areas_global_size = 0;
 
  die(
    "ERROR(interp_method_spmap::dist_read_cell_areas): "
    "YAC is built without the NetCDF support");
#else
 
  if ((io_rank_idx >= 0) && (io_rank_idx < num_io_ranks)) {
 
    // open file
    int ncid;
    yac_nc_open(filename, NC_NOWRITE, &ncid);
 
    // get variable id
    int varid;
    yac_nc_inq_varid(ncid, varname, &varid);
 
    // get dimension ids
    int ndims;
    int dimids[NC_MAX_VAR_DIMS];
    YAC_HANDLE_ERROR(
      nc_inq_var(ncid, varid, NULL, NULL, &ndims, dimids, NULL));
 
    YAC_ASSERT_F(
      (ndims == 1) || (ndims == 2),
      "ERROR(dist_read_cell_areas): "
      "invalid number of dimensions for cell area variable \"%s\" from "
      "file \"%s\" (is %d, but should be either 1 or 2)",
      varname, filename, ndims);
 
    // get size of dimensions
    size_t dimlen[2];
    *dist_cell_areas_global_size = 1;
    for (int i = 0; i < ndims; ++i) {
      YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dimids[i], &dimlen[i]));
      YAC_ASSERT_F(
        dimlen[i] > 0,
        "ERROR(dist_read_cell_areas): "
        "invalid dimension size for cell area variable \"%s\" from "
        "file \"%s\" (is %zu, should by > 0)",
        varname, filename, dimlen[i]);
      *dist_cell_areas_global_size *= dimlen[i];
    }
 
    // compute start/count
    // (in 2D case we have to round a little bit and then adjust the
    //  data afterwards)
    size_t start[2], count[2], offset, read_size;
    size_t local_start =
      (size_t)(
        ((long long)*dist_cell_areas_global_size * (long long)io_rank_idx) /
        (long long)num_io_ranks);
    size_t local_count =
      (size_t)(
        ((long long)*dist_cell_areas_global_size *
          (long long)(io_rank_idx+1)) / (long long)num_io_ranks) - local_start;
    if (ndims == 1) {
      start[0] = local_start;
      count[0] = local_count;
      offset = 0;
      read_size = local_count;
    } else {
      start[0] = local_start / dimlen[1];
      count[0] =
        (local_start + local_count + dimlen[1] - 1) / dimlen[1] - start[0];
      start[1] = 0;
      count[1] = dimlen[1];
      offset = local_start - start[0] * dimlen[1];
      read_size = count[0] * count[1];
    }
 
    // read data
    *dist_cell_areas = xmalloc(read_size * sizeof(**dist_cell_areas));
    YAC_HANDLE_ERROR(
      nc_get_vara_double(ncid, varid, start, count, *dist_cell_areas));
 
    // adjust data if necessary
    if (ndims == 2)
      memmove(
        *dist_cell_areas, *dist_cell_areas + offset,
        local_count * sizeof(**dist_cell_areas));
 
    // close file
    YAC_HANDLE_ERROR(nc_close(ncid));
 
  } else {
    *dist_cell_areas = xmalloc(1*sizeof(**dist_cell_areas));
    *dist_cell_areas_global_size = 0;
  }
 
  yac_mpi_call(
    MPI_Bcast(
      dist_cell_areas_global_size, 1, YAC_MPI_SIZE_T, io_ranks[0], comm),
    comm);
 
#endif // YAC_NETCDF_ENABLED
}
 
static void read_cell_areas(
  struct yac_const_basic_grid_data * basic_grid_data,
  struct yac_spmap_cell_area_file_config file_config, MPI_Comm comm,
  int * required_cell_areas, char const * type,
  double * cell_areas) {
 
  char const * routine = "read_cell_areas";
 
  size_t num_cells = basic_grid_data->count[YAC_LOC_CELL];
 
  // get io configuration
  int local_is_io, * io_ranks, num_io_ranks;
  yac_get_io_ranks(comm, &local_is_io, &io_ranks, &num_io_ranks);
 
  int comm_rank, comm_size;
  yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);
  yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
 
  int io_rank_idx = 0;
  while ((io_rank_idx < num_io_ranks) &&
          (comm_rank != io_ranks[io_rank_idx]))
    ++io_rank_idx;
  YAC_ASSERT_F(
    !local_is_io || (io_rank_idx < num_io_ranks),
    "ERROR(%s): unable to determine io_rank_idx", routine);
 
  double * dist_cell_areas;
  size_t dist_cell_areas_global_size;
 
  // read the data on the io processes
  dist_read_cell_areas(
    file_config.filename, file_config.varname, comm,
    io_ranks, io_rank_idx, num_io_ranks,
    &dist_cell_areas, &dist_cell_areas_global_size);
 
  // count the number of locally required cell areas
  size_t num_required_cell_areas = 0;
  for (size_t i = 0; i < num_cells; ++i)
    if (required_cell_areas[i]) ++num_required_cell_areas;
 
  size_t * global_idx = xmalloc(num_required_cell_areas * sizeof(*global_idx));
  size_t * reorder_idx =
    xmalloc(num_required_cell_areas * sizeof(*reorder_idx));
 
  size_t * sendcounts, * recvcounts, * sdispls, * rdispls;
  yac_get_comm_buffers(
    1, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
 
  // determine global indices of locally required cell areas
  yac_int const * global_cell_ids = basic_grid_data->ids[YAC_LOC_CELL];
  yac_int min_global_id = file_config.min_global_id;
  for (size_t i = 0, j = 0; i < num_cells; ++i) {
    if (required_cell_areas[i]) {
      YAC_ASSERT_F(
        global_cell_ids[i] >= min_global_id,
        "ERROR(%s): %s global id %" XT_INT_FMT " is smaller than "
        "the minimum global id provided by the user (%" XT_INT_FMT ")",
        routine, type, global_cell_ids[i], min_global_id);
      size_t idx = (size_t)(global_cell_ids[i] - min_global_id);
      YAC_ASSERT_F(
        idx < dist_cell_areas_global_size,
        "ERROR(%s): %s global id %" XT_INT_FMT " exceeds "
        "available size of array \"%s\" in file \"%s\" "
        "(min_global_id %" XT_INT_FMT ")", routine, type, global_cell_ids[i],
        file_config.varname, file_config.filename, file_config.min_global_id);
      global_idx[j] = idx;
      reorder_idx[j] = i;
      int dist_rank_idx =
        ((long long)idx * (long long)num_io_ranks +
         (long long)num_io_ranks - 1) / (long long)dist_cell_areas_global_size;
      sendcounts[io_ranks[dist_rank_idx]]++;
      ++j;
    }
  }
  free(io_ranks);
 
  yac_generate_alltoallv_args(
    1, sendcounts, recvcounts, sdispls, rdispls, comm);
 
  // sort required data by their global index
  yac_quicksort_index_size_t_size_t(
    global_idx, num_required_cell_areas, reorder_idx);
 
  // send required points to io processes
  size_t request_count = recvcounts[comm_size - 1] + rdispls[comm_size - 1];
  size_t * request_global_idx =
    xmalloc(request_count * sizeof(*request_global_idx));
  yac_alltoallv_p2p(
    global_idx, sendcounts, sdispls+1,
    request_global_idx, recvcounts, rdispls,
    sizeof(*global_idx), YAC_MPI_SIZE_T, comm, routine, __LINE__);
  free(global_idx);
 
  // pack requested cell areas
  double * requested_cell_areas =
    xmalloc(request_count * sizeof(*requested_cell_areas));
  size_t global_idx_offset =
    ((long long)io_rank_idx * (long long)dist_cell_areas_global_size) /
    (long long)num_io_ranks;
  for (size_t i = 0; i < request_count; ++i)
    requested_cell_areas[i] =
      dist_cell_areas[request_global_idx[i] - global_idx_offset];
  free(request_global_idx);
  free(dist_cell_areas);
 
  // return cell areas
  double * temp_cell_areas =
    xmalloc(num_required_cell_areas * sizeof(*temp_cell_areas));
  yac_alltoallv_p2p(
    requested_cell_areas, recvcounts, rdispls,
    temp_cell_areas, sendcounts, sdispls+1,
    sizeof(*requested_cell_areas), MPI_DOUBLE, comm, routine, __LINE__);
  free(requested_cell_areas);
 
  yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
 
  // unpack cell areas
  for (size_t i = 0; i < num_required_cell_areas; ++i)
    cell_areas[reorder_idx[i]] = temp_cell_areas[i];
 
  free(temp_cell_areas);
  free(reorder_idx);
}
 
static double * get_cell_areas(
  struct yac_const_basic_grid_data * basic_grid_data,
  char const * type, struct yac_spmap_cell_area_config cell_area_config,
  MPI_Comm comm, size_t const * required_points, size_t num_required_points) {
 
  size_t num_cells = basic_grid_data->count[YAC_LOC_CELL];
 
  // determine which cell areas are actually required
  int * required_cell_areas = xcalloc(num_cells, sizeof(*required_cell_areas));
  for (size_t i = 0; i < num_required_points; ++i)
    required_cell_areas[required_points[i]] = 1;
 
  // compute and scale cell areas
  double * cell_areas = xmalloc(num_cells * sizeof(*cell_areas));
 
  YAC_ASSERT_F(
    (cell_area_config.cell_area_provider == YAC_INTERP_SPMAP_CELL_AREA_YAC) ||
    (cell_area_config.cell_area_provider == YAC_INTERP_SPMAP_CELL_AREA_FILE),
    "ERROR(get_cell_areas): unsupported %s cell area origin", type);
 
  switch (cell_area_config.cell_area_provider) {
 
    default:
    case (YAC_INTERP_SPMAP_CELL_AREA_YAC): {
 
      double area_scale =
        cell_area_config.sphere_radius * cell_area_config.sphere_radius;
      compute_cell_areas(
        basic_grid_data, required_cell_areas, area_scale, type,
        cell_areas);
      break;
    }
    case (YAC_INTERP_SPMAP_CELL_AREA_FILE): {
 
      read_cell_areas(
        basic_grid_data, cell_area_config.file_config, comm,
        required_cell_areas, type, cell_areas);
    }
  };
 
  free(required_cell_areas);
 
  return cell_areas;
}
 
static void scale_weights(
  struct yac_interp_grid * interp_grid,
  struct yac_spmap_scale_config scale_config,
  size_t num_src_points, size_t const * src_points,
  size_t const * num_tgt_per_src, size_t const * tgt_points,
  size_t total_num_weights, double * weights) {
 
  YAC_ASSERT(
    (scale_config.type == YAC_INTERP_SPMAP_NONE) ||
    (scale_config.type == YAC_INTERP_SPMAP_SRCAREA) ||
    (scale_config.type == YAC_INTERP_SPMAP_INVTGTAREA) ||
    (scale_config.type == YAC_INTERP_SPMAP_FRACAREA),
    "ERROR(scale_weights): invalid scale_type")
 
  // if there is no scaling
  if (scale_config.type == YAC_INTERP_SPMAP_NONE) return;
 
  // get cell areas if they are required
  double const * src_cell_areas =
    ((scale_config.type == YAC_INTERP_SPMAP_SRCAREA) ||
     (scale_config.type == YAC_INTERP_SPMAP_FRACAREA))?
      get_cell_areas(
        yac_interp_grid_get_basic_grid_data_src(interp_grid),
        "source", scale_config.src, yac_interp_grid_get_MPI_Comm(interp_grid),
        src_points, total_num_weights):NULL;
  double const * tgt_cell_areas =
    ((scale_config.type == YAC_INTERP_SPMAP_INVTGTAREA) ||
     (scale_config.type == YAC_INTERP_SPMAP_FRACAREA))?
      get_cell_areas(
        yac_interp_grid_get_basic_grid_data_tgt(interp_grid),
        "target", scale_config.tgt, yac_interp_grid_get_MPI_Comm(interp_grid),
        tgt_points, total_num_weights):NULL;
 
#define SCALE_WEIGHTS( \
  SRC_CELL_AREA, TGT_CELL_AREA) \
{ \
  for (size_t i = 0, offset = 0; i < num_src_points; ++i) { \
    size_t curr_num_tgt = num_tgt_per_src[i]; \
    for (size_t j = 0; j < curr_num_tgt; ++j, ++offset) { \
      weights[offset] *= SRC_CELL_AREA / TGT_CELL_AREA; \
    } \
  } \
}
 
  switch (scale_config.type) {
    case(YAC_INTERP_SPMAP_SRCAREA):
      SCALE_WEIGHTS(src_cell_areas[src_points[offset]], 1.0)
      break;
    case(YAC_INTERP_SPMAP_INVTGTAREA):
      SCALE_WEIGHTS(1.0, tgt_cell_areas[tgt_points[offset]])
      break;
    default:
    case(YAC_INTERP_SPMAP_FRACAREA):
      SCALE_WEIGHTS(
        src_cell_areas[src_points[offset]], tgt_cell_areas[tgt_points[offset]])
      break;
  }
 
  free((void*)src_cell_areas);
  free((void*)tgt_cell_areas);
}
 
static void spread_src_data(
  struct yac_interp_grid * interp_grid, double spread_distance,
  enum yac_interp_spmap_weight_type weight_type,
  struct yac_spmap_scale_config scale_config,
  size_t num_src_points, size_t ** src_points_,
  size_t ** tgt_result_points_, double ** weights_,
  size_t * total_num_weights_) {
 
  // shortcut in case there is only a single "1.0" weight for all source points
  if ((spread_distance <= 0) && (scale_config.type == YAC_INTERP_SPMAP_NONE)) {
 
    *weights_ = NULL;
    *total_num_weights_ = num_src_points;
    return;
  }
 
  size_t * src_points = *src_points_;
  size_t * tgt_result_points = *tgt_result_points_;
 
  size_t * num_tgt_per_src =
    xmalloc(num_src_points * sizeof(*num_tgt_per_src));
  size_t total_num_weights;
 
  // search for additional target points if spread distance is bigger than 0.0
  if (spread_distance > 0.0) {
 
    struct sin_cos_angle inc_angle =
      sin_cos_angle_new(sin(spread_distance), cos(spread_distance));
    yac_const_coordinate_pointer tgt_field_coords =
      yac_interp_grid_get_tgt_field_coords(interp_grid);
 
    struct bounding_circle * search_bnd_circles =
      xmalloc(num_src_points * sizeof(*search_bnd_circles));
    for (size_t i = 0; i < num_src_points; ++i) {
      memcpy(
        search_bnd_circles[i].base_vector,
        tgt_field_coords[tgt_result_points[i]], sizeof(*tgt_field_coords));
      search_bnd_circles[i].inc_angle = inc_angle;
      search_bnd_circles[i].sq_crd = DBL_MAX;
    }
    size_t * spread_tgt_result_points = NULL;
 
    // do bounding circle search around found tgt points
    yac_interp_grid_do_bnd_circle_search_tgt(
      interp_grid, search_bnd_circles, num_src_points,
      &spread_tgt_result_points, num_tgt_per_src);
    free(search_bnd_circles);
 
    // remove target points which exceed the spread distance and only keep
    // targets that are connected to the original target point or other
    // target that have already been selected
    total_num_weights =
      check_tgt_result_points(
        interp_grid, spread_distance, num_src_points, tgt_result_points,
        num_tgt_per_src, spread_tgt_result_points);
    free((void*)tgt_result_points);
    tgt_result_points =
      xrealloc(
        spread_tgt_result_points,
        total_num_weights * sizeof(*spread_tgt_result_points));
 
    // adjust src_points (one source per target)
    size_t * new_src_points =
      xmalloc(total_num_weights * sizeof(*new_src_points));
    for (size_t i = 0, offset = 0; i < num_src_points; ++i)
      for (size_t j = 0, curr_src_point = src_points[i];
            j < num_tgt_per_src[i]; ++j, ++offset)
        new_src_points[offset] = curr_src_point;
    free((void*)src_points);
    src_points = new_src_points;
 
  } else {
 
    for (size_t i = 0; i < num_src_points; ++i) num_tgt_per_src[i] = 1;
    total_num_weights = num_src_points;
  }
 
  // compute weights
  double * weights =
    compute_weights(
      interp_grid, weight_type, num_src_points, src_points,
      num_tgt_per_src, total_num_weights, tgt_result_points);
 
  // scale weights
  scale_weights(
    interp_grid, scale_config, num_src_points, src_points, num_tgt_per_src,
    tgt_result_points, total_num_weights, weights);
 
  free(num_tgt_per_src);
 
  // set return values
  *tgt_result_points_ = tgt_result_points;
  *src_points_ = src_points;
  *weights_ = weights;
  *total_num_weights_ = total_num_weights;
}
 
static size_t do_search_spmap (struct interp_method * method,
                               struct yac_interp_grid * interp_grid,
                               size_t * tgt_points, size_t count,
                               struct yac_interp_weights * weights,
                               int * interpolation_complete) {
 
  if (*interpolation_complete) return 0;
 
  YAC_ASSERT(
    yac_interp_grid_get_num_src_fields(interp_grid) == 1,
    "ERROR(do_search_spmap): invalid number of source fields")
  YAC_ASSERT(
    yac_interp_grid_get_src_field_location(interp_grid, 0) == YAC_LOC_CELL,
    "ERROR(do_search_spmap): "
    "invalid source field location (has to be YAC_LOC_CELL)")
  YAC_ASSERT(
    yac_interp_grid_get_tgt_field_location(interp_grid) == YAC_LOC_CELL,
    "ERROR(do_search_spmap): "
    "invalid target field location (has to be YAC_LOC_CELL)")
 
  // get coordinates of all source points
  size_t * src_points;
  size_t num_src_points;
  yac_interp_grid_get_src_points(
    interp_grid, 0, &src_points, &num_src_points);
  yac_coordinate_pointer src_coords = xmalloc(num_src_points * sizeof(*src_coords));
  yac_interp_grid_get_src_coordinates(
    interp_grid, src_points, num_src_points, 0, src_coords);
 
  // search for matching tgt points
  size_t * tgt_result_points =
    xmalloc(num_src_points * sizeof(*tgt_result_points));
  yac_interp_grid_do_nnn_search_tgt(
    interp_grid, src_coords, num_src_points, 1, tgt_result_points,
    ((struct interp_method_spmap*)method)->max_search_distance);
 
  free(src_coords);
 
  // remove source points for which matching target point was found
  {
    size_t new_num_src_points = 0;
    for (size_t i = 0; i < num_src_points; ++i) {
      if (tgt_result_points[i] != SIZE_MAX) {
        if (i != new_num_src_points) {
          src_points[new_num_src_points] = src_points[i];
          tgt_result_points[new_num_src_points] =
            tgt_result_points[i];
        }
        ++new_num_src_points;
      }
    }
    num_src_points = new_num_src_points;
  }
 
  // spread the data from each source point to multiple target points
  // (weight_data is set to NULL if no spreading was applied)
  double * weight_data;
  size_t total_num_weights;
  spread_src_data(
    interp_grid, ((struct interp_method_spmap*)method)->spread_distance,
    ((struct interp_method_spmap*)method)->weight_type,
    ((struct interp_method_spmap*)method)->scale_config, num_src_points,
    &src_points, &tgt_result_points, &weight_data, &total_num_weights);
 
  // relocate source-target-point-pairs to dist owners of the respective
  // target points
  size_t result_count = total_num_weights;
  int to_tgt_owner = 1;
  yac_interp_grid_relocate_src_tgt_pairs(
    interp_grid, to_tgt_owner,
    0, &src_points, &tgt_result_points, &weight_data, &result_count);
  total_num_weights = result_count;
 
  // sort source-target-point-pairs by target points
  yac_quicksort_index_size_t_size_t_double(
    tgt_result_points, result_count, src_points, weight_data);
 
  // generate num_src_per_tgt and compact tgt_result_points
  size_t * num_src_per_tgt = xmalloc(result_count * sizeof(*num_src_per_tgt));
  size_t num_unique_tgt_result_points = 0;
  for (size_t i = 0; i < result_count;) {
    size_t prev_i = i;
    size_t curr_tgt = tgt_result_points[i];
    while ((i < result_count) && (curr_tgt == tgt_result_points[i])) ++i;
    num_src_per_tgt[num_unique_tgt_result_points] = i - prev_i;
    tgt_result_points[num_unique_tgt_result_points] = curr_tgt;
    ++num_unique_tgt_result_points;
  }
  result_count = num_unique_tgt_result_points;
  num_src_per_tgt =
    xrealloc(num_src_per_tgt, result_count * sizeof(*num_src_per_tgt));
  tgt_result_points =
    xrealloc(tgt_result_points, result_count * sizeof(*tgt_result_points));
 
  // match tgt_result_points with available target points
  qsort(tgt_points, count, sizeof(*tgt_points), compare_size_t);
  int * reorder_flag = xmalloc(count * sizeof(*tgt_points));
  {
    size_t j = 0;
    for (size_t i = 0; i < result_count; ++i) {
      size_t curr_result_tgt = tgt_result_points[i];
      while ((j < count) && (tgt_points[j] < curr_result_tgt))
        reorder_flag[j++] = 1;
      YAC_ASSERT(
        (j < count) && (curr_result_tgt == tgt_points[j]),
        "ERROR(do_search_spmap): "
        "required target points already in use or not available")
      reorder_flag[j++] = 0;
    }
    for (; j < count; ++j) reorder_flag[j] = 1;
  }
 
  // sort used target points to the beginning of the array
  yac_quicksort_index_int_size_t(reorder_flag, count, tgt_points);
  free(reorder_flag);
 
  struct remote_point * srcs =
    yac_interp_grid_get_src_remote_points(
      interp_grid, 0, src_points, total_num_weights);
  struct remote_points tgts = {
    .data =
      yac_interp_grid_get_tgt_remote_points(
        interp_grid, tgt_result_points, result_count),
    .count = result_count};
  free(tgt_result_points);
 
  // store results
  if (weight_data == NULL)
    yac_interp_weights_add_sum(
      weights, &tgts, num_src_per_tgt, srcs);
  else
    yac_interp_weights_add_wsum(
      weights, &tgts, num_src_per_tgt, srcs, weight_data);
 
  free(weight_data);
  free(src_points);
  free(tgts.data);
  free(srcs);
  free(num_src_per_tgt);
 
  return result_count;
}
 
static void check_cell_area_config(
  struct yac_spmap_cell_area_config * cell_area_config, char const * type) {
 
  YAC_ASSERT_F(
    (cell_area_config->cell_area_provider == YAC_INTERP_SPMAP_CELL_AREA_YAC) ||
    (cell_area_config->cell_area_provider == YAC_INTERP_SPMAP_CELL_AREA_FILE),
    "ERROR(check_cell_area_config): invalid %s cell_area_provider type", type);
 
  switch (cell_area_config->cell_area_provider) {
    default:
    case(YAC_INTERP_SPMAP_CELL_AREA_YAC): {
      YAC_ASSERT_F(
        cell_area_config->sphere_radius > 0.0,
        "ERROR(check_cell_area_config): invalid %s sphere_radius "
        "(has to be >= 0.0", type);
      cell_area_config->file_config.filename = NULL;
      cell_area_config->file_config.varname = NULL;
      cell_area_config->file_config.min_global_id = XT_INT_MAX;
      break;
    }
    case(YAC_INTERP_SPMAP_CELL_AREA_FILE): {
#ifndef YAC_NETCDF_ENABLED
      die(
        "ERROR(interp_method_spmap::check_cell_area_config): "
        "YAC is built without the NetCDF support");
#endif
      YAC_ASSERT_F(
        (cell_area_config->file_config.filename != NULL) &&
        (strlen(cell_area_config->file_config.filename) > 0),
        "ERROR(check_cell_area_config): invalid filename for %s areas", type);
      YAC_ASSERT_F(
        (cell_area_config->file_config.varname!= NULL) &&
        (strlen(cell_area_config->file_config.varname) > 0),
        "ERROR(check_cell_area_config): invalid varname for %s areas", type);
      cell_area_config->sphere_radius = 0.0;
      cell_area_config->file_config.filename =
        strdup(cell_area_config->file_config.filename);
      cell_area_config->file_config.varname =
        strdup(cell_area_config->file_config.varname);
      break;
    }
  }
}
 
struct interp_method * yac_interp_method_spmap_new(
  double spread_distance, double max_search_distance,
  enum yac_interp_spmap_weight_type weight_type,
  struct yac_spmap_scale_config scale_config) {
 
  struct interp_method_spmap * method = xmalloc(1 * sizeof(*method));
 
  method->vtable = &interp_method_spmap_vtable;
  method->spread_distance = spread_distance;
  method->max_search_distance =
    (max_search_distance == 0.0)?M_PI:max_search_distance;
  method->weight_type = weight_type;
  method->scale_config = scale_config;
 
  YAC_ASSERT(
    (spread_distance >= 0.0) && (spread_distance <= M_PI_2),
    "ERROR(yac_interp_method_spmap_new): invalid spread_distance "
    "(has to be >= 0 and <= PI/2")
 
  YAC_ASSERT(
    (max_search_distance >= 0.0) && (max_search_distance <= M_PI),
    "ERROR(yac_interp_method_spmap_new): invalid max_search_distance "
    "(has to be >= 0 and <= PI")
 
  check_cell_area_config(&(method->scale_config.src), "source");
  check_cell_area_config(&(method->scale_config.tgt), "target");
 
  return (struct interp_method*)method;
}
 
static void delete_spmap(struct interp_method * method) {
 
  struct interp_method_spmap * method_spmap =
    (struct interp_method_spmap*)(method);
 
  free((void*)method_spmap->scale_config.src.file_config.filename);
  free((void*)method_spmap->scale_config.src.file_config.varname);
  free((void*)method_spmap->scale_config.tgt.file_config.filename);
  free((void*)method_spmap->scale_config.tgt.file_config.varname);
  free(method);
}