test_interp_method_fixed_parallel.c

Go to the documentation of this file.

// Copyright (c) 2026 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#include <stdlib.h>
#include <math.h>
#include <string.h>
 
#include "tests.h"
#include "test_common.h"
#include "geometry.h"
#include "interpolation/methods/interp_method.h"
#include "interpolation/methods/interp_method_avg.h"
#include "interpolation/methods/interp_method_fixed.h"
#include "dist_grid_utils.h"
#include "yac_mpi.h"
 
#include <mpi.h>
#include <yaxt.h>
 
int split_comm_rank, split_comm_size;
 
enum yac_interp_weights_reorder_type reorder_types[] =
  {YAC_MAPPING_ON_SRC, YAC_MAPPING_ON_TGT};
size_t num_reorder_types = sizeof(reorder_types) / sizeof(reorder_types[0]);
static char const * grid_names[2] = {"src_grid", "tgt_grid"};
 
 
static void utest_fixed_value(double fixed_value)
{
  int is_tgt = split_comm_size == 1;
  double coordinates_x[2][4] = {{0.0,1.0,2.0,3.0}, {0.5,1.5,2.5}};
  double coordinates_y[2][3] = {{0.0,1.0,2.0}, {0.5,1.5,2.5}};
  size_t const num_cells[2][2] = {{3,2}, {2,2}};
  size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
  size_t local_count[2][2][2] = {{{2,2},{2,2}}, {{2,2}}};
  int with_halo = 0;
  for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
    coordinates_x[is_tgt][i] *= YAC_RAD;
  for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
    coordinates_y[is_tgt][i] *= YAC_RAD;
 
  struct yac_basic_grid_data grid_data =
    utest_generate_basic_grid_data_reg2d(
      coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
      local_start[is_tgt][split_comm_rank],
      local_count[is_tgt][split_comm_rank], with_halo);
 
  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])};
 
  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 = 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_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);
 
  struct interp_method * method_stack[] =
    {yac_interp_method_fixed_new(fixed_value), 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, NULL, 1, 1);
 
    {
      double * src_field = NULL;
      double ** src_fields_ptr = &src_field;
      double * tgt_field_buf = NULL;
 
      if (is_tgt) {
        tgt_field_buf =
          xmalloc(grid_data.num_vertices * sizeof(*tgt_field_buf));
        for (size_t j = 0; j < grid_data.num_vertices; ++j)
          tgt_field_buf[j] = -42.0;
      } else {
        src_field =
          xmalloc(grid_data.num_vertices * sizeof(*src_field));
        for (size_t j = 0; j < grid_data.num_vertices; ++j)
          src_field[j] = -666.0;
      }
 
      yac_interpolation_execute(interpolation, &src_fields_ptr, &tgt_field_buf);
 
      if (is_tgt)
        for (size_t j = 0; j < grid_data.num_vertices; ++j)
          if (memcmp(&tgt_field_buf[j], &fixed_value, sizeof(double)) != 0)
            PUT_ERR_F("wrong fixed interpolation result (expected %g, got %g)",
                      fixed_value, tgt_field_buf[j]);
 
      free(src_field);
      free(tgt_field_buf);
    }
 
    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]);
}
 
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);
  MPI_Barrier(MPI_COMM_WORLD);
 
  // Guard: this test requires exactly 3 MPI processes (2 src + 1 tgt)
  if (comm_size != 3) {
    PUT_ERR("ERROR: wrong number of processes");
    xt_finalize();
    MPI_Finalize();
    return TEST_EXIT_CODE;
  }
 
  MPI_Comm split_comm;
  MPI_Comm_split(
    MPI_COMM_WORLD, comm_rank < 2, 0, &split_comm);
 
  MPI_Comm_rank(split_comm, &split_comm_rank);
  MPI_Comm_size(split_comm, &split_comm_size);
 
  // Test 1: Basic fixed value = -1.0
  // fixed-value interpolation with two source processes and one target process
  // ALL target points must receive exactly -1.0
  // tests both YAC_MAPPING_ON_SRC and YAC_MAPPING_ON_TGT reorder types
  utest_fixed_value(-1.0);
 
  // Test 2: Fixed value = 0.0
  // Every target point must receive exactly 0.0
  utest_fixed_value(0.0);
 
  // Test 3: Fixed value = NaN
  // Every target point must receive NaN (verified via memcmp on bit pattern)
  utest_fixed_value(NAN);
 
  { // Test 4: Fixed as fallback in a method stack (avg first, then fixed)
    // Source grid has a mask that leaves the top-right corner uncovered.
    // Covered target points get avg-interpolated values; uncovered ones get -2.0.
    //
    // source mask (vertex mask, 4x3 vertices = 12 total):
    //   1-------1-------1-------1
    //   |       |       |       |
    //   |       |       |       |
    //   |       |       |       |
    //   1-------1-------1-------1
    //   |       |       |       |
    //   |       |       |       |
    //   |       |       |       |
    //   1-------1-------1-------0    <- vertex (3,0) masked out
    //
    // Source vertex ids: row 0: 0,1,2,3 / row 1: 4,5,6,7 / row 2: 8,9,10,11
    // Masked vertex: global id 3 (bottom-right corner).
    // Target corners near masked vertex will fall back to fixed = -2.0.
 
    int is_tgt = split_comm_size == 1;
    double coordinates_x[2][4] = {{0.0,1.0,2.0,3.0}, {0.5,1.5,2.5}};
    double coordinates_y[2][3] = {{0.0,1.0,2.0}, {0.5,1.5,2.5}};
    size_t const num_cells[2][2] = {{3,2}, {2,2}};
    size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
    size_t local_count[2][2][2] = {{{2,2},{2,2}}, {{2,2}}};
    // mask: all vertices valid except global id 3 (bottom-right of src grid)
    int global_src_mask[12] = {1, 1, 1, 0,
                               1, 1, 1, 1,
                               1, 1, 1, 1};
    int with_halo = 0;
    for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
      coordinates_x[is_tgt][i] *= YAC_RAD;
    for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
      coordinates_y[is_tgt][i] *= YAC_RAD;
 
    struct yac_basic_grid_data grid_data =
      utest_generate_basic_grid_data_reg2d(
        coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
        local_start[is_tgt][split_comm_rank],
        local_count[is_tgt][split_comm_rank], with_halo);
 
    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])};
 
    int * src_mask = NULL;
    if (!is_tgt) {
      src_mask = xmalloc(grid_data.num_vertices * sizeof(*src_mask));
      for (size_t i = 0; i < grid_data.num_vertices; ++i)
        src_mask[i] = global_src_mask[grid_data.vertex_ids[i]];
      yac_basic_grid_add_mask_nocpy(grids[0], YAC_LOC_CORNER, src_mask, NULL);
    }
 
    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 = SIZE_MAX, .masks_idx = 0}};
    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);
 
    // avg first (no partial coverage), then fixed = -2.0 as fallback
    struct interp_method * method_stack[] =
      {yac_interp_method_avg_new(YAC_INTERP_AVG_ARITHMETIC, 0),
       yac_interp_method_fixed_new(-2.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, NULL, 1, 1);
 
      {
        double * src_field = NULL;
        double ** src_fields_ptr = &src_field;
        double * tgt_field_buf = NULL;
        // Reference target values:
        // Global target vertex ids 0..8 in a 3x3 arrangement.
        // src vertex values = their global ids (0..11).
        // Vertices near the masked src corner (global id 3) cannot be
        // fully covered by avg => they fall back to -2.0.
        // Covered tgt vertices get arithmetic avg of surrounding src corners.
        // ref_global_tgt_field[9] matches the avg test with partial_coverage=0
        // except where masking forces fallback to -2.0.
        double ref_global_tgt_field[9] =
          {2.5, 3.5, -2.0, // <- last tgt vertex is in src cell with masked vertex
           6.5, 7.5, 8.5,
           -2.0, -2.0, -2.0}; // <- out of src domain
 
        if (is_tgt) {
          tgt_field_buf =
            xmalloc(grid_data.num_vertices * sizeof(*tgt_field_buf));
          for(size_t j = 0; j<grid_data.num_vertices; ++j)
            tgt_field_buf[j] = -42;
          double * ref_tgt_field =
            xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
          for (size_t j = 0; j < grid_data.num_vertices; ++j)
            ref_tgt_field[j] = ref_global_tgt_field[grid_data.vertex_ids[j]];
 
          yac_interpolation_execute(interpolation, &src_fields_ptr, &tgt_field_buf);
 
          for (size_t j = 0; j < grid_data.num_vertices; ++j)
            if (fabs(ref_tgt_field[j] - tgt_field_buf[j]) > 1e-9)
              PUT_ERR("Test 4: wrong interpolation result (avg+fixed fallback)");
 
          free(ref_tgt_field);
          free(tgt_field_buf);
        } else {
          src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
          for (size_t j = 0; j < grid_data.num_vertices; ++j)
            src_field[j] = (double)(grid_data.vertex_ids[j]);
 
          yac_interpolation_execute(interpolation, &src_fields_ptr, &tgt_field_buf);
 
          free(src_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]);
  }
 
  { // Test 5: collection_size > 1 with fixed value = 42.0
    // All target points in all 3 collections must receive 42.0
 
    int is_tgt = split_comm_size == 1;
    double coordinates_x[2][4] = {{0.0,1.0,2.0,3.0}, {0.5,1.5,2.5}};
    double coordinates_y[2][3] = {{0.0,1.0,2.0}, {0.5,1.5,2.5}};
    size_t const num_cells[2][2] = {{3,2}, {2,2}};
    size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
    size_t local_count[2][2][2] = {{{2,2},{2,2}}, {{2,2}}};
    int with_halo = 0;
    size_t const collection_size = 3;
    for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
      coordinates_x[is_tgt][i] *= YAC_RAD;
    for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
      coordinates_y[is_tgt][i] *= YAC_RAD;
 
    struct yac_basic_grid_data grid_data =
      utest_generate_basic_grid_data_reg2d(
        coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
        local_start[is_tgt][split_comm_rank],
        local_count[is_tgt][split_comm_rank], with_halo);
 
    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])};
 
    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 = 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_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);
 
    struct interp_method * method_stack[] =
      {yac_interp_method_fixed_new(42.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], collection_size,
            YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL,
            1, 1);
 
      {
        // src_data layout: [collection_idx][field_idx][local_idx]
        double *** src_data =
          xmalloc(collection_size * sizeof(*src_data));
        // tgt_data layout: [collection_idx][local_idx]
        double ** tgt_data =
          xmalloc(collection_size * sizeof(*tgt_data));
 
        for (size_t c = 0; c < collection_size; ++c) {
          // one src field per collection
          src_data[c] = xmalloc(1 * sizeof(**src_data));
          src_data[c][0] = NULL;
          tgt_data[c] = NULL;
        }
 
        if (is_tgt)
          for (size_t c = 0; c < collection_size; ++c){
            tgt_data[c] =
              xmalloc(grid_data.num_vertices * sizeof(*tgt_data[c]));
            for(size_t j = 0; j<grid_data.num_vertices; ++j)
              tgt_data[c][j] = -42;
          }
        else
          for (size_t c = 0; c < collection_size; ++c){
            src_data[c][0] =
              xmalloc(grid_data.num_vertices * sizeof(*src_data[c][0]));
            for(size_t j = 0; j<grid_data.num_vertices; ++j)
              src_data[c][0][j] = -666;
          }
 
        yac_interpolation_execute(interpolation, src_data, tgt_data);
 
        if (is_tgt)
          for (size_t c = 0; c < collection_size; ++c)
            for (size_t j = 0; j < grid_data.num_vertices; ++j)
              if (tgt_data[c][j] != 42.0)
                PUT_ERR("Test 5: wrong fixed interpolation result "
                        "(expected 42.0 in all collections)");
 
        for (size_t c = 0; c < collection_size; ++c) {
          free(src_data[c][0]);
          free(src_data[c]);
          free(tgt_data[c]);
        }
        free(src_data);
        free(tgt_data);
      }
 
      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();
  MPI_Finalize();
 
  return TEST_EXIT_CODE;
}