test_interp_operator_direct_mf.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_mf.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_mf(
  Xt_redist * redists, size_t num_src_fields,
  struct yac_collection_selection * sel,
  double *** src_fields_collection, double ** tgt_field_collection);
 
enum {
  MAX_COLLECTION_SIZE = 4,
  MAX_NUM_SRC_FIELDS = 2,
  NUM_SRC_POINTS = 2,
  NUM_TGT_POINTS = 4,
  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 != NUM_PROCS) {
    PUT_ERR("ERROR: test requires 2 processes");
    xt_finalize();
    MPI_Finalize();
    return TEST_EXIT_CODE;
  }
 
  // collection selection configurations
  struct {
    size_t N;
    size_t * indices;
  } collection_selection_configs[] =
    {{.N = 3, .indices = (size_t[]){0, 3, 2}},
     {.N = 4, .indices = (size_t[]){0, 1, 2, 3}},
     {.N = MAX_COLLECTION_SIZE, .indices = NULL}};
  enum {
    NUM_COLL_SEL_CONFIGS =
      sizeof(collection_selection_configs) / sizeof(collection_selection_configs[0])
  };
 
  // 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);
 
  // 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[MAX_NUM_SRC_FIELDS];
 
        //----------------------------------------------------------------------
        // reference data
        //----------------------------------------------------------------------
        int is_source;
        int is_target;
        int num_indices;
 
      } interp_data[NUM_PROCS];
    } interp_configs[] =
 
       // no redistribution
       // (processes are neither source nor target)
      {{.interp_data =
          {{.required_src_points = {{.num_indices = 0},
                                    {.num_indices = 0}},
            .is_source = 0,
            .is_target = 0},
           {.required_src_points = {{.num_indices = 0},
                                    {.num_indices = 0}},
            .is_source = 0,
            .is_target = 0}}},
       // sending all data of rank 0 to rank 1
       // (rank 0 is source and rank 1 is target)
       {.interp_data =
          {{.required_src_points = {{.num_indices = 0},
                                    {.num_indices = 0}},
            .is_source = 1,
            .is_target = 0},
           {.required_src_points = {{.num_indices = 2,
                                     .indices = (Xt_int[]){0, 1}},
                                    {.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}},
                                    {.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}},
                                    {.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}},
                                    {.num_indices = 2,
                                     .indices = (Xt_int[]){1, 3}}},
            .is_source = 1,
            .is_target = 1},
           {.required_src_points = {{.num_indices = 2,
                                     .indices = (Xt_int[]){0, 2}},
                                    {.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) {
 
      // create an array of redistributions for each source field
      Xt_redist redists[MAX_NUM_SRC_FIELDS];
      for (size_t f = 0; f < MAX_NUM_SRC_FIELDS; ++f) {
 
        Xt_idxlist required_source_idxlist =
          (interp_configs[interp_config_idx].
            interp_data[comm_rank].required_src_points[f].num_indices > 0)?
            xt_idxvec_new(
              interp_configs[interp_config_idx].
                interp_data[comm_rank].required_src_points[f].indices,
              interp_configs[interp_config_idx].
                interp_data[comm_rank].required_src_points[f].num_indices):
            xt_idxempty_new();
        Xt_xmap xmap =
          xt_xmap_all2all_new(
            src_idxlist, required_source_idxlist, MPI_COMM_WORLD);
        redists[f] = xt_redist_p2p_new(xmap, MPI_DOUBLE);
        xt_xmap_delete(xmap);
        xt_idxlist_delete(required_source_idxlist);
      }
 
      struct yac_interp_operator * interp =
        yac_interp_operator_direct_mf_new(sel, redists, MAX_NUM_SRC_FIELDS);
 
      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;
        }
 
        { // 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");
        }
 
        {// test 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,
            YAC_FRAC_MASK_NO_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_mf(
            redists, MAX_NUM_SRC_FIELDS, 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);
        }
 
        { // test 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,
            YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0);
          yac_interp_operator_execute_get(
            interp, tgt_field_collection,
            YAC_FRAC_MASK_NO_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_mf(
            redists, MAX_NUM_SRC_FIELDS, sel,
            src_fields_collection, tgt_field_collection);
 
          // 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");
 
          // free source and target fields
          utest_free_data(
            sel, src_fields_collection, src_fields_frac_mask_collection,
            tgt_field_collection);
        }
 
        { // test 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,
            YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0);
          yac_interp_operator_execute_get_async(
            interp, tgt_field_collection,
            YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0);
          yac_interp_operator_execute_wait(interp);
 
          // check data
          utest_check_direct_mf(
            redists, MAX_NUM_SRC_FIELDS, 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);
        }
 
        { // test 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,
            YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0);
          yac_interp_operator_execute_get_async(
            interp, tgt_field_collection,
            YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0);
          while(
            yac_interp_operator_execute_get_test(interp) ==
            YAC_INTERP_INCOMPLETE);
 
          // check data
          utest_check_direct_mf(
            redists, MAX_NUM_SRC_FIELDS, 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);
        }
 
        { // test 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,
            YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0);
          yac_interp_operator_execute_put(
            interp, src_fields_collection,
            src_fields_frac_mask_collection, 1,
            YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0);
          yac_interp_operator_execute_wait(interp);
 
          // check data
          utest_check_direct_mf(
            redists, MAX_NUM_SRC_FIELDS, 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);
        }
 
      } // use_copy
 
      yac_interp_operator_delete(interp);
 
      for (size_t f = 0; f < MAX_NUM_SRC_FIELDS; ++f)
        xt_redist_delete(redists[f]);
 
    } // NUM_INTERP_CONFIGS
 
    yac_collection_selection_delete(sel);
  } // NUM_COLL_SEL_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_mf(
  Xt_redist * redists, size_t num_src_fields,
  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];
 
    // for each source field
    for (size_t f = 0; f < num_src_fields; ++f) {
 
      // receive source points
      xt_redist_s_exchange1(
        redists[f], src_fields_collection[c_idx][f],
        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 multi-field interpolation");
      }
    }
  }
 
  // clean up
  for (size_t c = 0; c < collection_size; ++c) {
    free(ref_tgt_field_collection[c]);
  }
  free(ref_tgt_field_collection);
}