test_interp_method_parallel.c

Go to the documentation of this file.

// Copyright (c) 2024 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#include <stdlib.h>
#include <unistd.h>
#include <mpi.h>
#include <yaxt.h>
#include <string.h>
 
#include "tests.h"
#include "test_common.h"
#include "geometry.h"
#include "read_icon_grid.h"
#include "interpolation/methods/interp_method.h"
#include "interpolation/methods/interp_method_avg.h"
#include "interpolation/methods/interp_method_callback.h"
#include "interpolation/methods/interp_method_check.h"
#include "interpolation/methods/interp_method_conserv.h"
#include "interpolation/methods/interp_method_creep.h"
#include "interpolation/methods/interp_method_file.h"
#include "interpolation/methods/interp_method_fixed.h"
#include "interpolation/methods/interp_method_hcsbb.h"
#include "interpolation/methods/interp_method_ncc.h"
#include "interpolation/methods/interp_method_nnn.h"
#include "interpolation/methods/interp_method_spmap.h"
#include "weight_file_common.h"
#include "instrument.h"
#include "yac_mpi.h"
 
static void utest_compute_weights_callback(
  double const tgt_coords[3], int src_cell_id, size_t src_cell_idx,
  int const ** global_results_points, double ** result_weights,
  size_t * result_count, void * user_data);
 
typedef struct {
  int start_calls[YAC_INSTR_REGION_COUNT];
  int stop_calls[YAC_INSTR_REGION_COUNT];
  int last_start_depth[YAC_INSTR_REGION_COUNT];
  int bad_handle_calls;
} utest_instr_counters;
 
static void utest_instrumentation_enable(utest_instr_counters * instr_ctx);
static void utest_instrumentation_disable_and_check(
  utest_instr_counters * instr_ctx, int const * ref_start_calls);
static void utest_init_ref_start_calls(
  int * ref_start_calls, int default_start_calls);
 
int main(int argc, char** argv) {
 
  MPI_Init(NULL, NULL);
 
  xt_initialize(MPI_COMM_WORLD);
 
  int comm_rank, comm_size;
  MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
 
  utest_set_even_io_rank_list(MPI_COMM_WORLD);
 
  if (argc != 2) {
    PUT_ERR("ERROR: missing grid file directory");
    xt_finalize();
    MPI_Finalize();
    return TEST_EXIT_CODE;
  }
 
  char * filenames[2];
  char * grid_filenames[] ={"icon_grid_R02B01.nc", "icon_grid_R03B01.nc"};
  for (int i = 0; i < 2; ++i)
    filenames[i] =
      strcat(
        strcpy(
          malloc(strlen(argv[1]) + strlen(grid_filenames[i]) + 2), argv[1]),
        grid_filenames[i]);
 
  struct yac_basic_grid_data grid_data[2];
 
  for (int i = 0; i < 2; ++i)
    grid_data[i] =
      yac_read_icon_basic_grid_data_parallel(
        filenames[i], MPI_COMM_WORLD);
 
  struct yac_basic_grid * grids[2] =
    {yac_basic_grid_new(filenames[0], grid_data[0]),
     yac_basic_grid_new(filenames[1], grid_data[1])};
 
  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_CELL, .coordinates_idx = SIZE_MAX, .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_CELL, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
 
  struct yac_interp_grid * interp_grid =
    yac_interp_grid_new(grid_pair, filenames[0], filenames[1],
                        num_src_fields, src_fields, tgt_field);
 
  struct interp_method * method_stack[] =
    {yac_interp_method_fixed_new(-1.0), NULL};
 
  // the search over the whole method stack and the do_search routine of each
  // individual method are instrumented; the callbacks are registered only
  // around the search itself, so no unrelated region can contribute to the
  // counts checked afterwards
  utest_instr_counters instr_ctx;
  utest_instrumentation_enable(&instr_ctx);
 
  struct yac_interp_weights * weights =
    yac_interp_method_do_search(method_stack, interp_grid);
 
  { // only the fixed value method is in the stack
    int ref_start_calls[YAC_INSTR_REGION_COUNT];
    utest_init_ref_start_calls(ref_start_calls, 0);
    ref_start_calls[YAC_INSTR_INTERP_SEARCH_FIXED] = 1;
 
    utest_instrumentation_disable_and_check(&instr_ctx, ref_start_calls);
  }
 
  struct yac_interpolation * interpolation_src =
    yac_interp_weights_get_interpolation(
      weights, YAC_MAPPING_ON_SRC, 1,
      YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
  struct yac_interpolation * interpolation_tgt =
    yac_interp_weights_get_interpolation(
      weights, YAC_MAPPING_ON_TGT, 1,
      YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
 
  yac_interpolation_delete(interpolation_src);
  yac_interpolation_delete(interpolation_tgt);
 
  yac_interp_weights_delete(weights);
 
  yac_interp_method_delete(method_stack);
 
  yac_interp_grid_delete(interp_grid);
  yac_dist_grid_pair_delete(grid_pair);
 
  for (int i = 0; i < 2; ++i) {
    yac_basic_grid_delete(grids[i]);
    free(filenames[i]);
  }
 
 { // test interpolation_complete argument
 
    char const * weight_file_empty =
      "test_interp_method_parallel_weights_empty.nc";
    char const * weight_file_fixed =
      "test_interp_method_parallel_weights_fixed.nc";
 
    char const * src_grid_name = "src_grid";
    char const * tgt_grid_name = "tgt_grid";
 
    // create weight file (one empty and one with a fixed value)
    if (comm_rank == 0) {
 
      { // empty file
        int * tgt_indices = NULL;
        int * src_indices = NULL;
        double * weights = NULL;
        size_t num_links = 0;
        enum yac_location src_locations[] = {YAC_LOC_CELL};
        enum yac_location tgt_location = YAC_LOC_CELL;
        enum {
          NUM_SRC_FIELDS = sizeof(src_locations) / sizeof(src_locations[0])};
        int num_links_per_field[NUM_SRC_FIELDS] = {num_links};
        int * tgt_id_fixed = NULL;
        size_t num_fixed_tgt = 0;
        double * fixed_values = NULL;
        int * num_tgt_per_fixed_value = NULL;
        size_t num_fixed_values = 0;
 
        utest_write_weight_file(
          weight_file_empty, src_indices, tgt_indices, weights, num_links,
          src_locations, NUM_SRC_FIELDS, num_links_per_field, tgt_id_fixed,
          num_fixed_tgt, fixed_values, num_tgt_per_fixed_value,
          num_fixed_values, tgt_location, src_grid_name, tgt_grid_name);
      }
 
      { // file with fixed value
        int * tgt_indices = NULL;
        int * src_indices = NULL;
        double * weights = NULL;
        size_t num_links = 0;
        enum yac_location src_locations[] = {YAC_LOC_CELL};
        enum yac_location tgt_location = YAC_LOC_CELL;
        enum {
          NUM_SRC_FIELDS = sizeof(src_locations) / sizeof(src_locations[0])};
        int num_links_per_field[NUM_SRC_FIELDS] = {num_links};
        int tgt_id_fixed[] = {0, 1, 2, 3};
        size_t num_fixed_tgt = 4;
        double fixed_values[] = {999.0};
        int num_tgt_per_fixed_value[] = {4};
        size_t num_fixed_values = 1;
 
        utest_write_weight_file(
          weight_file_fixed, src_indices, tgt_indices, weights, num_links,
          src_locations, NUM_SRC_FIELDS, num_links_per_field, tgt_id_fixed,
          num_fixed_tgt, fixed_values, num_tgt_per_fixed_value,
          num_fixed_values, tgt_location, src_grid_name, tgt_grid_name);
      }
    }
    double coordinates_x[] = {0.0, 1.0, 2.0};
    double coordinates_y[] = {0.0, 1.0, 2.0};
    size_t num_vertices[2] = {3,3};
    int cyclic[2] = {0,0};
 
    struct yac_basic_grid * src_grid =
      yac_basic_grid_reg_2d_deg_new(
        src_grid_name, num_vertices, cyclic, coordinates_x, coordinates_y);
    struct yac_basic_grid * tgt_grid =
      yac_basic_grid_reg_2d_deg_new(
        tgt_grid_name, num_vertices, cyclic, coordinates_x, coordinates_y);
 
    double cell_coords[4][3];
    double cell_coords_x[] = {0.5, 1.5};
    double cell_coords_y[] = {0.5, 1.5};
    for (int i = 0, k = 0; i < 2; ++i)
      for (int j = 0; j < 2; ++j, ++k)
        LLtoXYZ_deg(cell_coords_x[j], cell_coords_y[i], cell_coords[k]);
    yac_basic_grid_add_coordinates(src_grid, YAC_LOC_CELL, cell_coords, 4);
    yac_basic_grid_add_coordinates(tgt_grid, YAC_LOC_CELL, cell_coords, 4);
 
    struct yac_dist_grid_pair * grid_pair =
      yac_dist_grid_pair_new(src_grid, tgt_grid, MPI_COMM_WORLD);
 
    struct yac_interp_field src_fields[] =
      {{.location = YAC_LOC_CELL, .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_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX};
 
    struct yac_interp_grid * interp_grid =
      yac_interp_grid_new(grid_pair, src_grid_name, tgt_grid_name,
                          num_src_fields, src_fields, tgt_field);
 
    //----------------------------------------
    // test generation of interpolation method
    //----------------------------------------
 
    struct interp_method * method_stack[] = {
      yac_interp_method_file_new( // this method stops the stack
        weight_file_empty, YAC_INTERP_FILE_MISSING_ERROR,
        YAC_INTERP_FILE_SUCCESS_STOP),
      yac_interp_method_avg_new(
        (enum yac_interp_avg_weight_type)
          YAC_INTERP_AVG_WEIGHT_TYPE_DEFAULT,
        YAC_INTERP_AVG_PARTIAL_COVERAGE_DEFAULT),
      yac_interp_method_callback_new(utest_compute_weights_callback, NULL),
      yac_interp_method_check_new( NULL, NULL, NULL, NULL),
      yac_interp_method_conserv_new(
        YAC_INTERP_CONSERV_ORDER_DEFAULT,
        YAC_INTERP_CONSERV_ENFORCED_CONSERV_DEFAULT,
        YAC_INTERP_CONSERV_PARTIAL_COVERAGE_DEFAULT,
        (enum yac_interp_method_conserv_normalisation)
          YAC_INTERP_CONSERV_NORMALISATION_DEFAULT),
      yac_interp_method_creep_new(YAC_INTERP_CREEP_DISTANCE_DEFAULT),
      yac_interp_method_file_new(
        weight_file_fixed, YAC_INTERP_FILE_MISSING_ERROR,
        YAC_INTERP_FILE_SUCCESS_CONT),
      yac_interp_method_fixed_new(YAC_INTERP_FIXED_VALUE_DEFAULT),
      yac_interp_method_hcsbb_new(),
      yac_interp_method_ncc_new(
        (enum yac_interp_ncc_weight_type)
          YAC_INTERP_NCC_WEIGHT_TYPE_DEFAULT,
        YAC_INTERP_NCC_PARTIAL_COVERAGE_DEFAULT),
      yac_interp_method_nnn_new(
        (struct yac_interp_method_nnn_config) {
          .type =
            (enum yac_interp_nnn_weight_type)YAC_INTERP_NNN_WEIGHTED_DEFAULT,
          .n = YAC_INTERP_NNN_N_DEFAULT,
          .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT,
          .data.gauss_scale = YAC_INTERP_NNN_GAUSS_SCALE_DEFAULT}),
      yac_interp_method_spmap_new(
        YAC_INTERP_SPMAP_DEFAULT_CONFIG, YAC_INTERP_SPMAP_OVERWRITE_DEFAULT),
      NULL};
 
    //-----------------
    // generate weights
    //-----------------
 
    utest_instr_counters instr_ctx;
    utest_instrumentation_enable(&instr_ctx);
 
    struct yac_interp_weights * weights =
      yac_interp_method_do_search(method_stack, interp_grid);
 
    if (yac_interp_weights_get_interp_count(weights) != 0)
      PUT_ERR(
        "error in handling of interpolation_complete argument by do_search");
 
    { // the first method in the stack stops the remaining stack, hence all
      // following methods return immediately without entering their region;
      // the check method is the only exception, because it ignores
      // interpolation_complete
      int ref_start_calls[YAC_INSTR_REGION_COUNT];
      utest_init_ref_start_calls(ref_start_calls, 0);
      ref_start_calls[YAC_INSTR_INTERP_SEARCH_FILE]  = 1;
      ref_start_calls[YAC_INSTR_INTERP_SEARCH_CHECK] = 1;
 
      utest_instrumentation_disable_and_check(&instr_ctx, ref_start_calls);
    }
 
    //-----------------------------------------------------------------
    // check the instrumentation for a method that actually searches
    //-----------------------------------------------------------------
 
    { // the distributed grid searches used by a method are instrumented as
      // well and have to be nested inside the region of that method
      struct interp_method * nnn_method_stack[] = {
        yac_interp_method_nnn_new(
          (struct yac_interp_method_nnn_config) {
            .type = (enum yac_interp_nnn_weight_type)YAC_INTERP_NNN_AVG,
            .n = 2,
            .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT,
            .data.gauss_scale = YAC_INTERP_NNN_GAUSS_SCALE_DEFAULT}),
        NULL};
 
      utest_instr_counters nnn_instr_ctx;
      utest_instrumentation_enable(&nnn_instr_ctx);
 
      struct yac_interp_weights * nnn_weights =
        yac_interp_method_do_search(nnn_method_stack, interp_grid);
 
      int ref_start_calls[YAC_INSTR_REGION_COUNT];
      utest_init_ref_start_calls(ref_start_calls, 0);
      ref_start_calls[YAC_INSTR_INTERP_SEARCH_NNN] = 1;
 
      utest_instrumentation_disable_and_check(&nnn_instr_ctx, ref_start_calls);
 
      int nnn_search_calls =
        nnn_instr_ctx.start_calls[YAC_INSTR_DIST_GRID_PAIR_NNN_SEARCH];
      int nnn_search_depth =
        nnn_instr_ctx.last_start_depth[YAC_INSTR_DIST_GRID_PAIR_NNN_SEARCH];
 
      if (nnn_search_calls != 1) {
        PUT_ERR_F(
          "error in instrumentation of yac_dist_grid_pair_do_nnn_search: "
          "region was entered %d times (expected 1)\n", nnn_search_calls);
      } else if (nnn_search_depth != 2) {
        PUT_ERR_F(
          "error in instrumentation of yac_dist_grid_pair_do_nnn_search: "
          "region was entered at depth %d (expected 2)\n", nnn_search_depth);
      }
 
      yac_interp_weights_delete(nnn_weights);
      yac_interp_method_delete(nnn_method_stack);
    }
 
    //---------------
    // cleanup
    //---------------
 
    yac_interp_method_delete(method_stack);
    yac_interp_weights_delete(weights);
    yac_interp_grid_delete(interp_grid);
    yac_dist_grid_pair_delete(grid_pair);
    yac_basic_grid_delete(tgt_grid);
    yac_basic_grid_delete(src_grid);
 
    if (comm_rank == 0) {
      unlink(weight_file_empty);
      unlink(weight_file_fixed);
    }
  }
 
  xt_finalize();
 
  MPI_Finalize();
 
  return TEST_EXIT_CODE;
}
 
static void utest_compute_weights_callback(
  double const tgt_coords[3], int src_cell_id, size_t src_cell_idx,
  int const ** global_results_points, double ** result_weights,
  size_t * result_count, void * user_data) {
 
  UNUSED(tgt_coords);
  UNUSED(src_cell_id);
  UNUSED(src_cell_idx);
  UNUSED(global_results_points);
  UNUSED(result_weights);
  UNUSED(result_count);
  UNUSED(user_data);
 
  PUT_ERR("ERROR: in handling interpolation_complete in do_search_callback");
}
 
static int utest_instr_handles[YAC_INSTR_REGION_COUNT];
 
static void * utest_instrument_create_cb(int id, void * user_data) {
 
  UNUSED(user_data);
 
  return (void *)&utest_instr_handles[id];
}
 
static void utest_instrument_start_cb(
  int id, void * handle, int depth, void * user_data) {
 
  utest_instr_counters * instr_ctx = (utest_instr_counters *)user_data;
 
  if (handle != (void *)&utest_instr_handles[id])
    instr_ctx->bad_handle_calls++;
  instr_ctx->start_calls[id]++;
  instr_ctx->last_start_depth[id] = depth;
}
 
static void utest_instrument_stop_cb(
  int id, void * handle, int depth, void * user_data) {
 
  UNUSED(depth);
 
  utest_instr_counters * instr_ctx = (utest_instr_counters *)user_data;
 
  if (handle != (void *)&utest_instr_handles[id])
    instr_ctx->bad_handle_calls++;
  instr_ctx->stop_calls[id]++;
}
 
static void utest_instrumentation_enable(utest_instr_counters * instr_ctx) {
 
  memset(instr_ctx, 0, sizeof(*instr_ctx));
  for (size_t i = 0; i < YAC_INSTR_REGION_COUNT; ++i)
    instr_ctx->last_start_depth[i] = -1;
 
  yac_instrument_set_callbacks(
    utest_instrument_create_cb, utest_instrument_start_cb,
    utest_instrument_stop_cb, NULL, instr_ctx);
}
 
static void utest_init_ref_start_calls(
  int * ref_start_calls, int default_start_calls) {
 
  // only the interpolation search regions (a contiguous block in
  // YAC_INSTRUMENT_REGION_TABLE) are checked; the regions entered by the
  // methods themselves (e.g. the distributed grid searches) are not, because
  // how often they are entered depends on the domain decomposition
  for (size_t i = 0; i < YAC_INSTR_REGION_COUNT; ++i)
    ref_start_calls[i] = -1;
  for (int i = YAC_INSTR_INTERP_SEARCH;
       i <= YAC_INSTR_INTERP_SEARCH_SPMAP; ++i)
    ref_start_calls[i] = default_start_calls;
 
  // the search over the method stack is entered exactly once per search
  ref_start_calls[YAC_INSTR_INTERP_SEARCH] = 1;
}
 
static void utest_instrumentation_disable_and_check(
  utest_instr_counters * instr_ctx, int const * ref_start_calls) {
 
  yac_instrument_clear_callbacks();
 
  if (instr_ctx->bad_handle_calls != 0)
    PUT_ERR(
      "error in instrumentation of do_search "
      "(handle does not match the one created for the region)");
 
  for (size_t i = 0; i < YAC_INSTR_REGION_COUNT; ++i) {
 
    char const * region_name =
      yac_instrument_region_name((yac_instrument_region_id)i);
 
    if (instr_ctx->start_calls[i] != instr_ctx->stop_calls[i])
      PUT_ERR_F(
        "error in instrumentation of do_search: region \"%s\" was started "
        "%d times, but stopped %d times\n",
        region_name, instr_ctx->start_calls[i], instr_ctx->stop_calls[i]);
 
    if (ref_start_calls[i] == -1) continue;
 
    if (instr_ctx->start_calls[i] != ref_start_calls[i])
      PUT_ERR_F(
        "error in instrumentation of do_search: region \"%s\" was started "
        "%d times (expected %d)\n",
        region_name, instr_ctx->start_calls[i], ref_start_calls[i]);
 
    // the search over the method stack is the outermost region; the regions
    // of the individual methods are nested directly inside of it
    if (instr_ctx->start_calls[i] > 0) {
      int ref_depth = (i == (size_t)YAC_INSTR_INTERP_SEARCH)?0:1;
      if (instr_ctx->last_start_depth[i] != ref_depth)
        PUT_ERR_F(
          "error in instrumentation of do_search: region \"%s\" was started "
          "at depth %d (expected %d)\n",
          region_name, instr_ctx->last_start_depth[i], ref_depth);
    }
  }
}