YetAnotherCoupler 3.2.0_a
Loading...
Searching...
No Matches
test_interp_method_nnn_parallel2.c

A test for the parallel nearest neighbour interpolation method.

// Copyright (c) 2024 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "tests.h"
#include "interp_method.h"
#include "dist_grid_utils.h"
#include "yac_mpi.h"
#include "geometry.h"
#include <mpi.h>
#include <yaxt.h>
#include <netcdf.h>
char const src_grid_name[] = "src_grid";
char const tgt_grid_name[] = "tgt_grid";
static void target_main(MPI_Comm global_comm, MPI_Comm target_comm);
static void source_main(MPI_Comm global_comm, MPI_Comm source_comm);
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);
if (comm_size != 5) {
PUT_ERR("ERROR: wrong number of processes");
xt_finalize();
MPI_Finalize();
return TEST_EXIT_CODE;
}
int tgt_flag = comm_rank < 1;
MPI_Comm split_comm;
MPI_Comm_split(MPI_COMM_WORLD, tgt_flag, 0, &split_comm);
if (tgt_flag) target_main(MPI_COMM_WORLD, split_comm);
else source_main(MPI_COMM_WORLD, split_comm);
MPI_Comm_free(&split_comm);
xt_finalize();
MPI_Finalize();
return TEST_EXIT_CODE;
}
static void source_main(MPI_Comm global_comm,
MPI_Comm source_comm) {
// 1 and 5 nearest neighbours per target point with fixed value fallback
// corner and cell ids for a 7 x 7 grid (x == target point position)
// 56-----57-----58-----59-----60-----61-----62-----63
// | x| x| x| x| x| x| x|
// | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
// | | | | | | | |
// 48-----49-----50-----51-----52-----53-----54-----55
// | x| x| x| x| x| x| x|
// | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
// | | | | | | | |
// 40-----41-----42-----43-----44-----45-----46-----47
// | x| x| x| x| x| x| x|
// | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
// | | | | | | | |
// 32-----33-----34-----35-----36-----37-----38-----39
// | x| x| x| x| x| x| x|
// | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
// | | | | | | | |
// 24-----25-----26-----27-----28-----29-----30-----31
// | x| x| x| x| x| x| x|
// | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
// | | | | | | | |
// 16-----17-----18-----19-----20-----21-----22-----23
// | x| x| x| x| x| x| x|
// | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
// | | | | | | | |
// 08-----09-----10-----11-----12-----13-----14-----15
// | x| x| x| x| x| x| x|
// | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
// | | | | | | | |
// 00-----01-----02-----03-----04-----05-----06-----07
//
// the grid is distributed among the processes as follows:
// (index == process)
//
// 3---3---3---3---3---3---3---3
// | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
// 3---3---3---3---3---3---3---3
// | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
// 3---3---3---1---2---2---3---3
// | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
// 1---1---1---2---2---2---2---2
// | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
// 1---1---1---1---2---2---2---2
// | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
// 1---1---1---0---0---0---2---2
// | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
// 0---0---0---0---0---0---0---0
// | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
// 0---0---0---0---0---0---0---0
//
// the source mask looks as follows (# == masked point)
//
// +---+---#---#---#---#---#---#
// | | | | | | | |
// +---+---#---#---#---#---#---#
// | | | | | | | |
// #---#---#---#---#---#---#---#
// | | | | | | | |
// #---#---#---#---#---#---#---#
// | | | | | | | |
// #---#---#---#---#---#---#---#
// | | | | | | | |
// #---#---#---#---#---#---#---#
// | | | | | | | |
// +---#---#---#---#---#---#---+
// | | | | | | | |
// +---+---#---#---#---#---+---+
//
// the target mask looks as follows (# == masked point)
//
// +---+---+---+---+---+---+---+
// | | | | | | | |
// +---+---+---+---+---+---+---+
// | | | | | | | |
// +---+---#---#---#---#---+---+
// | | | | | | | |
// +---+---#---#---#---#---+---+
// | | | | | | | |
// +---+---#---#---#---#---+---+
// | | | | | | | |
// +---+---#---#---#---#---+---+
// | | | | | | | |
// +---+---+---+---+---+---+---+
// | | | | | | | |
// +---+---+---+---+---+---+---+
//--------------------
// set up global search
//--------------------
int my_source_rank;
MPI_Comm_rank(source_comm, &my_source_rank);
double coordinates_x[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
double coordinates_y[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
size_t const num_cells[2] = {7,7};
size_t local_start[4][2] = {{0,0},{0,2},{3,2},{0,5}};
size_t local_count[4][2] = {{7,2},{3,3},{4,3},{7,2}};
int with_halo = 0;
int global_corner_mask[8][8] = {
{1,1,0,0,0,0,1,1},
{1,0,0,0,0,0,0,1},
{0,0,0,0,0,0,0,0},
{0,0,0,0,0,0,0,0},
{0,0,0,0,0,0,0,0},
{0,0,0,0,0,0,0,0},
{1,1,0,0,0,0,0,0},
{1,1,0,0,0,0,0,0}};
for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
struct yac_basic_grid_data grid_data =
yac_generate_basic_grid_data_reg2d(
coordinates_x, coordinates_y, num_cells,
local_start[my_source_rank], local_count[my_source_rank], with_halo);
int * src_corner_mask =
xmalloc(grid_data.num_vertices * sizeof(*src_corner_mask));
for (size_t i = 0; i < grid_data.num_vertices; ++i)
src_corner_mask[i] =
((int*)(&(global_corner_mask[0][0])))[grid_data.vertex_ids[i]];
struct yac_basic_grid * src_grid =
yac_basic_grid_new(src_grid_name, grid_data);
src_grid, YAC_LOC_CORNER, src_corner_mask, NULL);
struct yac_basic_grid * tgt_grid = yac_basic_grid_empty_new(tgt_grid_name);
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_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_CELL, .coordinates_idx = 0, .masks_idx = 0};
struct yac_interp_grid * interp_grid =
struct interp_method * method_stack[2] = {
(struct yac_nnn_config){.type = YAC_INTERP_NNN_AVG, .n = 4}), NULL};
struct yac_interp_weights * weights =
yac_interp_method_do_search(method_stack, interp_grid);
yac_interp_method_delete(method_stack);
enum yac_interp_weights_reorder_type reorder_type[2] =
for (size_t i = 0; i < 2; ++i) {
struct yac_interpolation * interpolation =
weights, reorder_type[i], 1, YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0);
// check generated interpolation
{
double * src_field = NULL;
double ** src_fields = &src_field;
src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
for (size_t i = 0; i < grid_data.num_vertices; ++i)
src_field[i] = (double)(grid_data.vertex_ids[i]);
yac_interpolation_execute(interpolation, &src_fields, NULL);
free(src_field);
}
yac_interpolation_delete(interpolation);
}
yac_interp_grid_delete(interp_grid);
}
static void target_main(MPI_Comm global_comm,
MPI_Comm target_comm) {
double coordinates_x[] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
double coordinates_y[] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
double cell_coordinates_x[] = {0.75,1.75,2.75,3.75,4.75,5.75,6.75};
double cell_coordinates_y[] = {0.75,1.75,2.75,3.75,4.75,5.75,6.75};
yac_coordinate_pointer cell_coordinates = xmalloc(49 * sizeof(*cell_coordinates));
size_t const num_cells[2] = {7,7};
size_t local_start[2] = {0,0};
size_t local_count[2] = {7,7};
int with_halo = 0;
int global_cell_mask[7][7] = {
{1,1,1,1,1,1,1},
{1,1,1,1,1,1,1},
{1,1,0,0,0,1,1},
{1,1,0,0,0,1,1},
{1,1,0,0,0,1,1},
{1,1,1,1,1,1,1},
{1,1,1,1,1,1,1}};
for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
for (size_t i = 0, k = 0; i < num_cells[1]; ++i)
for (size_t j = 0; j < num_cells[0]; ++j, ++k)
cell_coordinates_x[j], cell_coordinates_y[i], cell_coordinates[k]);
struct yac_basic_grid_data grid_data =
yac_generate_basic_grid_data_reg2d(
coordinates_x, coordinates_y, num_cells,
local_start, local_count, with_halo);
int * tgt_cell_mask =
xmalloc(grid_data.num_cells * sizeof(*tgt_cell_mask));
for (size_t i = 0; i < grid_data.num_cells; ++i)
tgt_cell_mask[i] =
((int*)(&(global_cell_mask[0][0])))[grid_data.cell_ids[i]];
struct yac_basic_grid * tgt_grid =
yac_basic_grid_new(tgt_grid_name, grid_data);
tgt_grid, YAC_LOC_CELL, tgt_cell_mask, NULL);
tgt_grid, YAC_LOC_CELL, cell_coordinates);
struct yac_basic_grid * src_grid = yac_basic_grid_empty_new(src_grid_name);
struct yac_dist_grid_pair * grid_pair =
yac_dist_grid_pair_new(tgt_grid, src_grid, 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_CELL, .coordinates_idx = 0, .masks_idx = 0};
struct yac_interp_grid * interp_grid =
struct interp_method * method_stack[2] = {
(struct yac_nnn_config){.type = YAC_INTERP_NNN_AVG, .n = 4}), NULL};
struct yac_interp_weights * weights =
yac_interp_method_do_search(method_stack, interp_grid);
yac_interp_method_delete(method_stack);
enum yac_interp_weights_reorder_type reorder_type[2] =
for (size_t i = 0; i < 2; ++i) {
struct yac_interpolation * interpolation =
weights, reorder_type[i], 1, YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0);
// check generated interpolation
{
double * tgt_field = NULL;
double ref_tgt_field[49] =
{(0+1+8+49)/4.0,(0+1+8+6)/4.0,(0+1+8+6)/4.0,(1+6+7+15)/4.0,
(1+6+7+15)/4.0,(1+6+7+15)/4.0,(1+6+7+15)/4.0,
(0+1+8+49)/4.0,(0+1+8+49)/4.0,(0+1+8+6)/4.0,(1+6+7+15)/4.0,
(1+6+7+15)/4.0,(1+6+7+15)/4.0,(1+6+7+15)/4.0,
(8+1+0+49)/4.0,(8+1+0+49)/4.0,-1,-1,
-1,(1+6+7+15)/4.0,(1+6+7+15)/4.0,
(48+49+8+57)/4.0,(48+49+8+57)/4.0,-1,-1,
-1,(6+7+15+49)/4.0,(6+7+15+49)/4.0,
(48+49+56+57)/4.0,(48+49+56+57)/4.0,-1,-1,
-1,(6+7+15+49)/4.0,(6+7+15+49)/4.0,
(48+49+56+57)/4.0,(48+49+56+57)/4.0,(48+49+56+57)/4.0,(48+49+56+57)/4.0,
(48+49+56+57)/4.0,(48+49+15+57)/4.0,(6+7+15+49)/4.0,
(48+49+56+57)/4.0,(48+49+56+57)/4.0,(48+49+56+57)/4.0,(48+49+56+57)/4.0,
(48+49+56+57)/4.0,(48+49+56+57)/4.0,(49+57+56+15)/4.0};
tgt_field = xmalloc(grid_data.num_cells * sizeof(*tgt_field));
for (size_t j = 0; j < grid_data.num_cells; ++j) tgt_field[j] = -1.0;
yac_interpolation_execute(interpolation, NULL, &tgt_field);
for (size_t j = 0; j < grid_data.num_cells; ++j)
if (fabs(ref_tgt_field[j] - tgt_field[j]) > 1e-9)
PUT_ERR("wrong interpolation result");
free(tgt_field);
}
yac_interpolation_delete(interpolation);
}
yac_interp_grid_delete(interp_grid);
}
struct yac_basic_grid * yac_basic_grid_new(char const *name, struct yac_basic_grid_data grid_data)
Definition basic_grid.c:45
size_t yac_basic_grid_add_coordinates_nocpy(struct yac_basic_grid *grid, enum yac_location location, yac_coordinate_pointer coordinates)
Definition basic_grid.c:203
struct yac_basic_grid * yac_basic_grid_empty_new(char const *name)
Definition basic_grid.c:58
void yac_basic_grid_delete(struct yac_basic_grid *grid)
Definition basic_grid.c:65
size_t yac_basic_grid_add_mask_nocpy(struct yac_basic_grid *grid, enum yac_location location, int const *mask, char const *mask_name)
Definition basic_grid.c:253
void yac_dist_grid_pair_delete(struct yac_dist_grid_pair *grid_pair)
Definition dist_grid.c:2957
struct yac_dist_grid_pair * yac_dist_grid_pair_new(struct yac_basic_grid *grid_a, struct yac_basic_grid *grid_b, MPI_Comm comm)
Definition dist_grid.c:2714
#define YAC_RAD
Definition geometry.h:30
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:304
void yac_interp_grid_delete(struct yac_interp_grid *interp_grid)
struct yac_interp_grid * yac_interp_grid_new(struct yac_dist_grid_pair *grid_pair, char const *src_grid_name, char const *tgt_grid_name, size_t num_src_fields, struct yac_interp_field const *src_fields, struct yac_interp_field const tgt_field)
Definition interp_grid.c:31
void yac_interp_method_delete(struct interp_method **method)
struct yac_interp_weights * yac_interp_method_do_search(struct interp_method **method, struct yac_interp_grid *interp_grid)
struct interp_method * yac_interp_method_nnn_new(struct yac_nnn_config config)
@ YAC_INTERP_NNN_AVG
average of n source points
struct yac_interpolation * yac_interp_weights_get_interpolation(struct yac_interp_weights *weights, enum yac_interp_weights_reorder_type reorder, size_t collection_size, double frac_mask_fallback_value, double scaling_factor, double scaling_summand)
void yac_interp_weights_delete(struct yac_interp_weights *weights)
yac_interp_weights_reorder_type
@ YAC_MAPPING_ON_TGT
weights will be applied at target processes
@ YAC_MAPPING_ON_SRC
weights will be appied at source processes
void yac_interpolation_execute(struct yac_interpolation *interp, double ***src_fields, double **tgt_field)
void yac_interpolation_delete(struct yac_interpolation *interp)
double const YAC_FRAC_MASK_NO_VALUE
@ YAC_LOC_CORNER
Definition location.h:15
@ YAC_LOC_CELL
Definition location.h:14
#define xmalloc(size)
Definition ppm_xfuncs.h:66
enum yac_location location
Definition basic_grid.h:18
struct yac_interp_field tgt_field
Definition interp_grid.c:26
size_t num_src_fields
Definition interp_grid.c:27
struct yac_dist_grid_pair * grid_pair
Definition interp_grid.c:25
char * src_grid_name
Definition interp_grid.c:23
struct yac_interp_field src_fields[]
Definition interp_grid.c:28
char * tgt_grid_name
Definition interp_grid.c:24
int main(int argc, char **argv)
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19