test_interp_operator_direct.c

Go to the documentation of this file.

// Copyright (c) 2025 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
 
#include <mpi.h>
#include <yaxt.h>
 
#include "tests.h"
#include "interpolation/interpolation.h"
#include "yac_mpi.h"
#include "test_common.h"
#include "interpolation/operators/interp_operator_direct.h"
#include "collection_selection_internal.h"
 
static void utest_init_data(
  struct yac_collection_selection * sel, double **** src_fields_collection,
  double **** src_fields_frac_mask_collection, double *** tgt_field_collection);
static void utest_free_data(
  struct yac_collection_selection * sel, double *** src_fields_collection,
  double *** src_fields_frac_mask_collection, double ** tgt_field_collection);
static void utest_check_direct(
  Xt_redist redist, struct yac_collection_selection * sel,
  double *** src_fields_collection, double ** tgt_field_collection);
 
enum {
  MAX_COLLECTION_SIZE = 4,
  MAX_NUM_SRC_FIELDS = 1,
  NUM_SRC_POINTS = 2,
  NUM_TGT_POINTS = 2,
  NUM_PROCS = 2,
};
#define TGT_UNSET_VALUE (-1.0)
 
int main(void) {
 
  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);
 
  if (comm_size != 2) {
    PUT_ERR("ERROR: test requires 2 processes");
    xt_finalize();
    MPI_Finalize();
    return TEST_EXIT_CODE;
  }
 
  // generate yaxt index list containing global ids of all locally available
  // source points (each source field contains the same global ids)
  // each process holds the following global ids:
  //  i + comm_rank * NUM_SRC_POINTS
  //  with i in [0..NUM_SRC_POINTS]
  Xt_int src_global_ids[NUM_SRC_POINTS];
  for (size_t src_idx = 0; src_idx < NUM_SRC_POINTS; ++src_idx)
    src_global_ids[src_idx] =
      (Xt_int)(src_idx + comm_rank * NUM_SRC_POINTS);
  Xt_idxlist src_idxlist =
    xt_idxvec_new(src_global_ids, NUM_SRC_POINTS);
 
  // collection selection configurations
  struct {
    size_t N;
    size_t * indices;
  } collection_selection_configs[] =
    {{.N = 3, .indices = (size_t[]){0, 2, 3}},
     {.N = 4, .indices = (size_t[]){0, 1, 2, 3}},
     {.N = 2, .indices = NULL}};
  enum {
    NUM_COLL_SEL_CONFIGS =
      sizeof(collection_selection_configs)/
      sizeof(collection_selection_configs[0]) };
 
  // loop over collection selections
  for (size_t sel_idx = 0; sel_idx < NUM_COLL_SEL_CONFIGS; ++sel_idx) {
 
    struct yac_collection_selection * sel =
      yac_collection_selection_new(
        collection_selection_configs[sel_idx].N,
        collection_selection_configs[sel_idx].indices);
 
    // interpolation configurations
    // We assume that the global ids of the source source points of each field
    // are:
    // * for rank 0: [0, 1]
    // * for rank 1: [2, 3]
    struct {
      struct {
 
        //----------------------------------------------------------------------
        // required source points
        //----------------------------------------------------------------------
        struct {
          int num_indices;  // number of source points to be received
          Xt_int * indices; // global ids of required source points
        } required_src_points;
 
        //----------------------------------------------------------------------
        // reference data
        //----------------------------------------------------------------------
        int is_source;
        int is_target;
 
      } interp_data[NUM_PROCS];
    } interp_configs[] =
       // no redistribution
       // (processes are neither source nor target)
      {{.interp_data =
          {{.required_src_points = {.num_indices = 0, .indices = NULL},
            .is_source = 0,
            .is_target = 0},
           {.required_src_points = {.num_indices = 0, .indices = NULL},
            .is_source = 0,
            .is_target = 0}}},
       // rank 0 sends all its data to rank 1
       // (rank 0 is source and rank 1 is target)
       {.interp_data =
          {{.required_src_points = {.num_indices = 0, .indices = NULL},
            .is_source = 1,
            .is_target = 0},
           {.required_src_points = {.num_indices = 2,
                                    .indices = (Xt_int[]){0, 1}},
            .is_source = 0,
            .is_target = 1}}},
       // target data on each rank consists of its local source data
       // (all processes are both source and target)
       {.interp_data =
          {{.required_src_points = {.num_indices = 2,
                                    .indices = (Xt_int[]){0, 1}},
            .is_source = 1,
            .is_target = 1},
           {.required_src_points = {.num_indices = 2,
                                    .indices = (Xt_int[]){2, 3}},
            .is_source = 1,
            .is_target = 1}}},
       // target data on each rank consists of remote and local source data
       // (all processes are both source and target)
       {.interp_data =
          {{.required_src_points = {.num_indices = 2,
                                    .indices = (Xt_int[]){0, 2}},
            .is_source = 1,
            .is_target = 1},
           {.required_src_points = {.num_indices = 2,
                                    .indices = (Xt_int[]){1, 3}},
            .is_source = 1,
            .is_target = 1}}}};
    enum {
      NUM_INTERP_CONFIGS = sizeof(interp_configs) / sizeof(interp_configs[0])};
 
    // loop over all interpolation configurations
    for (size_t interp_config_idx = 0; interp_config_idx < NUM_INTERP_CONFIGS;
         ++interp_config_idx) {
 
      Xt_idxlist required_source_idxlist =
        (interp_configs[interp_config_idx].
           interp_data[comm_rank].required_src_points.num_indices > 0)?
          xt_idxvec_new(
            interp_configs[interp_config_idx].
              interp_data[comm_rank].required_src_points.indices,
            interp_configs[interp_config_idx].
              interp_data[comm_rank].required_src_points.num_indices):
          xt_idxempty_new();
 
      Xt_xmap xmap =
        xt_xmap_all2all_new(
          src_idxlist, required_source_idxlist, MPI_COMM_WORLD);
      Xt_redist redist = xt_redist_p2p_new(xmap, MPI_DOUBLE);
 
      xt_xmap_delete(xmap);
      xt_idxlist_delete(required_source_idxlist);
 
      // create direct interpolation operator
      struct yac_interp_operator * interp =
        yac_interp_operator_direct_new(sel, redist);
 
      for (int use_copy = 0; use_copy <= 1; ++use_copy) {
 
        if (use_copy) {
          struct yac_interp_operator * interp_copy =
            yac_interp_operator_copy(interp);
          yac_interp_operator_delete(interp);
          interp = interp_copy;
        }
 
        for (int with_frac_mask = 0; with_frac_mask <= 1; ++with_frac_mask) {
 
          // if without fractional masking, set fallback value to
          // YAC_FRAC_MASK_NO_VALUE, which deacitavtes it
          double frac_mask_fallback_value =
            with_frac_mask?-10.0:YAC_FRAC_MASK_NO_VALUE;
 
          { // check is_source and is_target
 
            if (yac_interp_operator_is_source(interp) !=
                interp_configs[interp_config_idx].
                  interp_data[comm_rank].is_source)
              PUT_ERR("ERROR in is_source");
 
            if (yac_interp_operator_is_target(interp) !=
                interp_configs[interp_config_idx].
                  interp_data[comm_rank].is_target)
              PUT_ERR("ERROR in is_target");
          }
 
          { // synchronous execute
 
            // initialise source and target fields
            double *** src_fields_collection;
            double *** src_fields_frac_mask_collection;
            double ** tgt_field_collection;
            utest_init_data(
              sel, &src_fields_collection, &src_fields_frac_mask_collection,
              &tgt_field_collection);
 
            yac_interp_operator_execute(
              interp, src_fields_collection, src_fields_frac_mask_collection,
              tgt_field_collection, frac_mask_fallback_value, 1.0, 0.0);
 
            // there should be no open put or get operation
            if (yac_interp_operator_execute_put_test(interp) !=
                YAC_INTERP_COMPLETE)
              PUT_ERR("ERROR in execute_put_test");
            if (yac_interp_operator_execute_get_test(interp) !=
                YAC_INTERP_COMPLETE)
              PUT_ERR("ERROR in execute_get_test");
 
            // check results
            utest_check_direct(
              redist, sel, src_fields_collection, tgt_field_collection);
 
            // free source and target fields
            utest_free_data(
              sel, src_fields_collection, src_fields_frac_mask_collection,
              tgt_field_collection);
          }
 
          { // independent put + get
 
            // initialise source and target fields
            double *** src_fields_collection;
            double *** src_fields_frac_mask_collection;
            double ** tgt_field_collection;
 
            utest_init_data(
              sel, &src_fields_collection, &src_fields_frac_mask_collection,
              &tgt_field_collection);
            yac_interp_operator_execute_put(
              interp, src_fields_collection, src_fields_frac_mask_collection,
              1, frac_mask_fallback_value, 1.0, 0.0);
            yac_interp_operator_execute_get(
              interp, tgt_field_collection,
              frac_mask_fallback_value, 1.0, 0.0);
 
            // there should be no open put or get operation
            if (yac_interp_operator_execute_put_test(interp) !=
                YAC_INTERP_COMPLETE)
              PUT_ERR("ERROR in execute_put_test");
            if (yac_interp_operator_execute_get_test(interp) !=
                YAC_INTERP_COMPLETE)
              PUT_ERR("ERROR in execute_get_test");
 
            // check data
            utest_check_direct(
              redist, sel, src_fields_collection, tgt_field_collection);
 
            // free source and target fields
            utest_free_data(
              sel, src_fields_collection, src_fields_frac_mask_collection,
              tgt_field_collection);
          }
 
          { // independent put + async get + wait
 
            // initialise source and target fields
            double *** src_fields_collection;
            double *** src_fields_frac_mask_collection;
            double ** tgt_field_collection;
 
            utest_init_data(
              sel, &src_fields_collection, &src_fields_frac_mask_collection,
              &tgt_field_collection);
 
            yac_interp_operator_execute_put(
              interp, src_fields_collection, src_fields_frac_mask_collection,
              1, frac_mask_fallback_value, 1.0, 0.0);
            yac_interp_operator_execute_get_async(
              interp, tgt_field_collection,
              frac_mask_fallback_value, 1.0, 0.0);
            yac_interp_operator_execute_wait(interp);
 
            // there should be no open put or get operation
            if (yac_interp_operator_execute_put_test(interp) !=
                YAC_INTERP_COMPLETE)
              PUT_ERR("ERROR in execute_put_test");
            if (yac_interp_operator_execute_get_test(interp) !=
                YAC_INTERP_COMPLETE)
              PUT_ERR("ERROR in execute_get_test");
 
            // check data
            utest_check_direct(
              redist, sel, src_fields_collection, tgt_field_collection);
 
            // free source and target fields
            utest_free_data(
              sel, src_fields_collection, src_fields_frac_mask_collection,
              tgt_field_collection);
          }
 
          { // independent put + async get + test_get-loop
 
            // initialise source and target fields
            double *** src_fields_collection;
            double *** src_fields_frac_mask_collection;
            double ** tgt_field_collection;
            utest_init_data(
              sel, &src_fields_collection, &src_fields_frac_mask_collection,
              &tgt_field_collection);
 
            yac_interp_operator_execute_put(
              interp, src_fields_collection,
              src_fields_frac_mask_collection,
              1, frac_mask_fallback_value, 1.0, 0.0);
            yac_interp_operator_execute_get_async(
              interp, tgt_field_collection,
              frac_mask_fallback_value, 1.0, 0.0);
            while(
              yac_interp_operator_execute_get_test(interp) ==
              YAC_INTERP_INCOMPLETE);
 
            // there should be no open put or get operation
            if (yac_interp_operator_execute_put_test(interp) !=
                YAC_INTERP_COMPLETE)
              PUT_ERR("ERROR in execute_put_test");
            if (yac_interp_operator_execute_get_test(interp) !=
                YAC_INTERP_COMPLETE)
              PUT_ERR("ERROR in execute_get_test");
 
            // check data
            utest_check_direct(
              redist, sel, src_fields_collection, tgt_field_collection);
 
            // free source and target fields
            utest_free_data(
              sel, src_fields_collection, src_fields_frac_mask_collection,
              tgt_field_collection);
          }
 
          { // independent async get + put + wait
 
            // initialise source and target fields
            double *** src_fields_collection;
            double *** src_fields_frac_mask_collection;
            double ** tgt_field_collection;
            utest_init_data(
              sel, &src_fields_collection, &src_fields_frac_mask_collection,
              &tgt_field_collection);
 
            yac_interp_operator_execute_get_async(
              interp, tgt_field_collection,
              frac_mask_fallback_value, 1.0, 0.0);
            yac_interp_operator_execute_put(
              interp, src_fields_collection, src_fields_frac_mask_collection,
              1, frac_mask_fallback_value, 1.0, 0.0);
            yac_interp_operator_execute_wait(interp);
 
            // there should be no open put or get operation
            if (yac_interp_operator_execute_put_test(interp) !=
                YAC_INTERP_COMPLETE)
              PUT_ERR("ERROR in execute_put_test");
            if (yac_interp_operator_execute_get_test(interp) !=
                YAC_INTERP_COMPLETE)
              PUT_ERR("ERROR in execute_get_test");
 
            // check data
            utest_check_direct(
              redist, sel, src_fields_collection, tgt_field_collection);
 
            // free source and target fields
            utest_free_data(
              sel, src_fields_collection, src_fields_frac_mask_collection,
              tgt_field_collection);
          }
 
        } // with_frac_mask
 
      } // use_copy
 
      yac_interp_operator_delete(interp);
      xt_redist_delete(redist);
 
    } // redist configs
 
    yac_collection_selection_delete(sel);
 
  } // collection selection configs
 
  xt_idxlist_delete(src_idxlist);
 
  xt_finalize();
  MPI_Finalize();
 
  return TEST_EXIT_CODE;
}
 
static void utest_init_data(
  struct yac_collection_selection * sel, double **** src_fields_collection,
  double **** src_fields_frac_mask_collection,
  double *** tgt_field_collection) {
 
  size_t collection_size = yac_collection_selection_get_collection_size(sel);
 
  int rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
  *src_fields_collection =
    xmalloc(MAX_COLLECTION_SIZE * sizeof(**src_fields_collection));
  *src_fields_frac_mask_collection =
    xmalloc(MAX_COLLECTION_SIZE * sizeof(**src_fields_frac_mask_collection));
 
  for (size_t c = 0; c < MAX_COLLECTION_SIZE; ++c) {
 
    (*src_fields_collection)[c] =
      xmalloc(MAX_NUM_SRC_FIELDS * sizeof(***src_fields_collection));
    (*src_fields_frac_mask_collection)[c] =
      xmalloc(MAX_NUM_SRC_FIELDS * sizeof(***src_fields_frac_mask_collection));
 
    for (size_t f = 0; f < MAX_NUM_SRC_FIELDS; ++f) {
 
      (*src_fields_collection)[c][f] =
        xmalloc(NUM_SRC_POINTS * sizeof(***src_fields_collection));
      (*src_fields_frac_mask_collection)[c][f] =
        xmalloc(NUM_SRC_POINTS * sizeof(***src_fields_frac_mask_collection));
 
      for (size_t i = 0; i < NUM_SRC_POINTS; ++i) {
        (*src_fields_collection)[c][f][i] =
          (double)(c*1000 + f*100 + rank*10 + i);
        (*src_fields_frac_mask_collection)[c][f][i] = 1.0;
      }
    }
  }
 
  *tgt_field_collection =
    xmalloc(collection_size * sizeof(**tgt_field_collection));
 
  for (size_t c = 0; c < collection_size; ++c) {
 
    (*tgt_field_collection)[c] =
      xmalloc(NUM_TGT_POINTS * sizeof(**tgt_field_collection));
 
    for (size_t i = 0; i < NUM_TGT_POINTS; ++i) {
      (*tgt_field_collection)[c][i] = TGT_UNSET_VALUE;
    }
  }
}
 
static void utest_free_data(
  struct yac_collection_selection * sel, double *** src_fields_collection,
  double *** src_fields_frac_mask_collection, double ** tgt_field_collection) {
 
  size_t collection_size = yac_collection_selection_get_collection_size(sel);
 
  for (size_t c = 0; c < MAX_COLLECTION_SIZE; ++c) {
 
    for (size_t f = 0; f < MAX_NUM_SRC_FIELDS; ++f) {
 
      free(src_fields_collection[c][f]);
      free(src_fields_frac_mask_collection[c][f]);
    }
 
    free(src_fields_collection[c]);
    free(src_fields_frac_mask_collection[c]);
  }
 
  free(src_fields_collection);
  free(src_fields_frac_mask_collection);
 
  for (size_t c = 0; c < collection_size; ++c) {
 
    free(tgt_field_collection[c]);
  }
 
  free(tgt_field_collection);
}
 
static void utest_check_direct(
  Xt_redist redist, struct yac_collection_selection * sel,
  double *** src_fields_collection, double ** tgt_field_collection) {
 
  size_t collection_size =
    yac_collection_selection_get_collection_size(sel);
  size_t const * collection_indices =
    yac_collection_selection_get_indices(sel);
 
  // initialise reference target field with dummy values
  // and generate array containg pointers to the field associated to
  // each collection
  double ** ref_tgt_field_collection =
    xmalloc(collection_size * sizeof(*ref_tgt_field_collection));
  for (size_t c = 0; c < collection_size; ++c) {
    ref_tgt_field_collection[c] =
      xmalloc(NUM_TGT_POINTS * sizeof(**ref_tgt_field_collection));
    for (size_t i = 0; i < NUM_TGT_POINTS; ++i) {
      ref_tgt_field_collection[c][i] = TGT_UNSET_VALUE;
    }
  }
 
  // for each collection entry
  for (size_t c = 0; c < collection_size; ++c) {
 
    // get collection index (if non-contiguous) otherwise use c
    size_t c_idx = (collection_indices == NULL)?c:collection_indices[c];
 
    // receive source points
    xt_redist_s_exchange1(
      redist, src_fields_collection[c_idx][0], ref_tgt_field_collection[c]);
  }
 
  for (size_t c = 0; c < collection_size; ++c) {
    for (size_t i = 0; i < NUM_TGT_POINTS; ++i) {
      if (tgt_field_collection[c][i] != ref_tgt_field_collection[c][i])
        PUT_ERR("wrong data in direct interpolation");
    }
  }
 
  for (size_t c = 0; c < collection_size; ++c) {
    free(ref_tgt_field_collection[c]);
  }
  free(ref_tgt_field_collection);
}