A test for the parallel average interpolation method.
#include <stdlib.h>
#include "tests.h"
#include "test_common.h"
#include "dist_grid_utils.h"
#include <mpi.h>
#include <yaxt.h>
#include <netcdf.h>
size_t num_reorder_types = sizeof(reorder_types) / sizeof(reorder_types[0]);
static char const * grid_names[2] = {"src_grid", "tgt_grid"};
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 != 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);
int split_comm_rank, split_comm_size;
MPI_Comm_rank(split_comm, &split_comm_rank);
MPI_Comm_size(split_comm, &split_comm_size);
{
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;
yac_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);
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
for (int partial_coverage = 0; partial_coverage <= 1; ++partial_coverage) {
for (size_t i = 0; i < num_reorder_types; ++i) {
weights, reorder_types[i], 1,
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
double * ref_tgt_field = NULL;
double ref_global_tgt_field[9] =
{2.5, 3.5, 4.5, 6.5, 7.5, 8.5, -1, -1, -1};
if (is_tgt) {
ref_tgt_field[i] = ref_global_tgt_field[grid_data.
vertex_ids[i]];
} else {
}
if (is_tgt)
if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-9)
PUT_ERR("wrong interpolation result");
free(src_field);
free(tgt_field);
free(ref_tgt_field);
}
}
}
}
{
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 global_src_mask[12] = {1, 1, 1, 0,
1, 1, 1, 1,
0, 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;
yac_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);
int * src_mask = NULL;
if (!is_tgt) {
src_mask[i] = global_src_mask[grid_data.
vertex_ids[i]];
}
{{.location =
YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
{.location =
YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
for (size_t i = 0; i < num_reorder_types; ++i) {
weights, reorder_types[i], 1,
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
double * ref_tgt_field = NULL;
double ref_global_tgt_field[9] =
{2.5, 3.5, -1, -1, 7.5, 8.5, -1, -1, -1};
if (is_tgt) {
ref_tgt_field[i] = ref_global_tgt_field[grid_data.
vertex_ids[i]];
} else {
}
if (is_tgt)
if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-9)
PUT_ERR("wrong interpolation result");
free(src_field);
free(tgt_field);
free(ref_tgt_field);
}
}
}
{
int is_tgt = split_comm_size == 1;
double coordinates_x[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
double coordinates_y[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
size_t const num_cells[2][2] = {{3,3}, {3,3}};
size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
size_t local_count[2][2][2] = {{{2,3},{2,3}}, {{3,3}}};
int global_src_mask[16] = {0, 1, 1, 0,
1, 1, 1, 1,
1, 1, 1, 1,
0, 1, 1, 0};
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;
yac_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);
int * src_mask = NULL;
if (!is_tgt) {
src_mask[i] = global_src_mask[grid_data.
vertex_ids[i]];
}
if (is_tgt) {
double cell_coordinates_x[3] = {0.75,1.75,2.75};
double cell_coordinates_y[3] = {0.75,1.75,2.75};
tgt_cell_coordinates =
xmalloc(9 *
sizeof(*tgt_cell_coordinates));
for (int i = 0, k = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j, ++k)
cell_coordinates_x[j], cell_coordinates_y[i],
tgt_cell_coordinates[k]);
}
{{.location =
YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
{.location =
YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX};
for (size_t i = 0; i < num_reorder_types; ++i) {
weights, reorder_types[i], 1,
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
double * ref_tgt_field = NULL;
double ref_global_tgt_field[9] =
{-1.0,
0.14963* 1.0 + 0.20075*( 2.0+ 5.0) + 0.44888* 6.0,
-1.0,
0.14963* 4.0 + 0.20075*( 5.0+ 8.0) + 0.44888* 9.0,
0.14963* 5.0 + 0.20075*( 6.0+ 9.0) + 0.44888*10.0,
0.14963* 6.0 + 0.20075*( 7.0+10.0) + 0.44888*11.0,
-1.0,
0.14963* 9.0 + 0.20075*(10.0+13.0) + 0.44888*14.0,
-1.0};
if (is_tgt) {
ref_tgt_field[i] = ref_global_tgt_field[grid_data.
vertex_ids[i]];
} else {
}
if (is_tgt)
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-3)
PUT_ERR("wrong interpolation result");
free(src_field);
free(tgt_field);
free(ref_tgt_field);
}
}
}
{
int is_tgt = split_comm_size == 1;
double coordinates_x[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
double coordinates_y[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
size_t const num_cells[2][2] = {{3,3}, {3,3}};
size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
size_t local_count[2][2][2] = {{{2,3},{2,3}}, {{3,3}}};
int global_src_mask[16] = {0, 1, 0, 0,
1, 1, 0, 0,
1, 0, 1, 1,
0, 1, 1, 0};
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;
yac_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);
int * src_mask = NULL;
if (!is_tgt) {
src_mask[i] = global_src_mask[grid_data.
vertex_ids[i]];
}
if (is_tgt) {
double cell_coordinates_x[3] = {0.75,1.75,2.75};
double cell_coordinates_y[3] = {0.75,1.75,2.75};
tgt_cell_coordinates =
xmalloc(9 *
sizeof(*tgt_cell_coordinates));
for (int i = 0, k = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j, ++k)
cell_coordinates_x[j], cell_coordinates_y[i],
tgt_cell_coordinates[k]);
}
{{.location =
YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
{.location =
YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX};
for (size_t i = 0; i < num_reorder_types; ++i) {
weights, reorder_types[i], 1,
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
double * ref_tgt_field = NULL;
double inv_dist_a = sqrt(8.0);
double inv_dist_b = sqrt(8.0/5.0);
double inv_dist_c = sqrt(8.0/9.0);
double ref_global_tgt_field[9] =
{(inv_dist_b*( 1.0+ 4.0) + inv_dist_a* 5.0)/(inv_dist_a + 2.0 * inv_dist_b),
(inv_dist_c* 1.0 + inv_dist_b* 5.0)/(inv_dist_b + inv_dist_c),
-1.0,
(inv_dist_c* 4.0 + inv_dist_b*( 5.0+ 8.0))/(2.0 * inv_dist_b + inv_dist_c),
(inv_dist_c* 5.0 + inv_dist_a*10.0)/(inv_dist_a + inv_dist_c),
(inv_dist_b*10.0 + inv_dist_a*11.0)/(inv_dist_a + inv_dist_b),
(inv_dist_c* 8.0 + inv_dist_a*13.0)/(inv_dist_a + inv_dist_c),
(inv_dist_b*(10.0+13.0) + inv_dist_a*14.0)/(inv_dist_a + 2.0 * inv_dist_b),
(inv_dist_c*10.0 + inv_dist_b*(11.0+14.0))/(2.0 * inv_dist_b + inv_dist_c)};
if (is_tgt) {
ref_tgt_field[i] = ref_global_tgt_field[grid_data.
vertex_ids[i]];
} else {
}
if (is_tgt)
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-3)
PUT_ERR("wrong interpolation result");
free(src_field);
free(tgt_field);
free(ref_tgt_field);
}
}
}
{
int is_tgt = split_comm_size == 1;
double coordinates_x[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
double coordinates_y[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
size_t const num_cells[2][2] = {{3,3}, {3,3}};
size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
size_t local_count[2][2][2] = {{{2,3},{2,3}}, {{3,3}}};
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;
yac_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);
if (is_tgt) {
double cell_coordinates_x[3] = {1.0,2.0,3.0};
double cell_coordinates_y[3] = {1.0,2.0,3.0};
tgt_cell_coordinates =
xmalloc(9 *
sizeof(*tgt_cell_coordinates));
for (int i = 0, k = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j, ++k)
cell_coordinates_x[j], cell_coordinates_y[i],
tgt_cell_coordinates[k]);
}
{{.location =
YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
{.location =
YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX};
for (size_t i = 0; i < num_reorder_types; ++i) {
weights, reorder_types[i], 1,
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
double * ref_tgt_field = NULL;
double ref_global_tgt_field[9] =
{5.0, 6.0, 7.0, 9.0, 10.0, 11.0, 13.0, 14.0, 15.0};
if (is_tgt) {
ref_tgt_field[i] = ref_global_tgt_field[grid_data.
vertex_ids[i]];
} else {
}
if (is_tgt)
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-9)
PUT_ERR("wrong interpolation result");
free(src_field);
free(tgt_field);
free(ref_tgt_field);
}
}
}
{
int is_tgt = split_comm_size == 1;
double coordinates_x[8] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
double coordinates_y[8] = {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[2][2][2] = {{{0,0},{0,4}}, {{0,0}}};
size_t local_count[2][2][2] = {{{7,4},{7,3}}, {{7,7}}};
int with_halo = 0;
for (size_t i = 0; i <= num_cells[0]; ++i)
for (size_t i = 0; i <= num_cells[1]; ++i)
yac_generate_basic_grid_data_reg2d(
local_start[is_tgt][split_comm_rank],
local_count[is_tgt][split_comm_rank], with_halo);
if (!is_tgt) {
src_point_coordinates =
for (
size_t i = 0; i < grid_data.
num_cells; ++i) {
double * middle_point = src_point_coordinates[i];
for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
size_t * curr_vertices =
for (size_t j = 0; j < curr_num_vertices; ++j) {
double * curr_vertex_coord =
for (size_t k = 0; k < 3; ++k)
middle_point[k] += curr_vertex_coord[k];
}
}
}
{{.location =
YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX}};
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
{.location =
YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
for (size_t i = 0; i < num_reorder_types; ++i) {
weights, reorder_types[i], 1,
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
double * ref_tgt_field = NULL;
double ref_global_tgt_field[64] =
{-4,-4,-4,-4,-4,-4,-4,-4,
-4, 0+ 1+ 7+ 8, 1+ 2+ 8+ 9, 2+ 3+ 9+10, 3+ 4+10+11, 4+ 5+11+12, 5+ 6+12+13,-4,
-4, 7+ 8+14+15, 8+ 9+15+16, 9+10+16+17,10+11+17+18,11+12+18+19,12+13+19+20,-4,
-4,14+15+21+22,15+16+22+23,16+17+23+24,17+18+24+25,18+19+25+26,19+20+26+27,-4,
-4,21+22+28+29,22+23+29+30,23+24+30+31,24+25+31+32,25+26+32+33,26+27+33+34,-4,
-4,28+29+35+36,29+30+36+37,30+31+37+38,31+32+38+39,32+33+39+40,33+34+40+41,-4,
-4,35+36+42+43,36+37+43+44,37+38+44+45,38+39+45+46,39+40+46+47,40+41+47+48,-4,
-4,-4,-4,-4,-4,-4,-4,-4};
for (size_t i = 0;
i < sizeof(ref_global_tgt_field) / sizeof(ref_global_tgt_field[0]);
++i)
ref_global_tgt_field[i] /= 4.0;
if (is_tgt) {
ref_tgt_field[i] = ref_global_tgt_field[grid_data.
vertex_ids[i]];
} else {
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
src_field[i] = (
double)(grid_data.
cell_ids[i]);
}
if (is_tgt)
if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-9)
PUT_ERR("wrong interpolation result");
free(src_field);
free(tgt_field);
free(ref_tgt_field);
}
}
}
{
int is_tgt = split_comm_size == 1;
double src_coordinates_x[2][2] = {{-0.5,0.5}, {0.5,-0.5}};
double src_coordinates_y[2][2][2] =
{{{-0.5,0.5}, {0.5,-0.5}}, {{89,90}, {90,89}}};
double src_coordinates_x_unstruct[4];
double src_coordinates_y_unstruct[4];
size_t const src_num_cells[2] = {1,1};
size_t src_local_start[2] = {0,0};
size_t src_local_count[2] = {1,1};
int src_cell_mask[4] = {1,1, 0,1};
{{{0,1,2,3},{2,3,0,1}},{{1,0,3,2},{3,2,1,0}}};
double tgt_coordinates_x[5] = {-0.5,-0.25,0.0,0.25,0.5};
double tgt_coordinates_y[2][5] =
{{-0.5,-0.25,0.0,0.25,0.5}, {89.00,89.25,89.50,89.75,90.00}};
size_t const tgt_num_cells[2] = {4,4};
size_t tgt_local_start[2] = {0,0};
size_t tgt_local_count[2] = {4,4};
for (size_t i = 0; i < 2 * 2; ++i)
(&(src_coordinates_x[0][0]))[i] *=
YAC_RAD;
for (size_t i = 0; i < 2 * 2 * 2; ++i)
(&(src_coordinates_y[0][0][0]))[i] *=
YAC_RAD;
for (
size_t i = 0; i < 5; ++i) tgt_coordinates_x[i] *=
YAC_RAD;
for (size_t i = 0; i < 2 * 5; ++i)
(&(tgt_coordinates_y[0][0]))[i] *=
YAC_RAD;
for (int partial_coverage = 0; partial_coverage < 2; ++partial_coverage) {
for (int is_reg = 0; is_reg < 2; ++is_reg) {
for (int at_pole = 0; at_pole <= is_reg; ++at_pole) {
for (int with_mask = 0; with_mask < 2; ++with_mask) {
for (int order_x = 0; order_x < 2; ++order_x) {
for (int order_y = 0; order_y < 2; ++order_y) {
if (is_tgt) {
grid_data =
yac_generate_basic_grid_data_reg2d(
tgt_coordinates_x, tgt_coordinates_y[at_pole],
tgt_num_cells, tgt_local_start, tgt_local_count, 0);
} else {
if (is_reg) {
grid_data =
yac_generate_basic_grid_data_reg2d(
src_coordinates_x[order_x],
src_coordinates_y[at_pole][order_y], src_num_cells,
src_local_start, src_local_count, 0);
} else {
for (int i = 0; i < 4; ++i) {
src_coordinates_x_unstruct[i] =
src_coordinates_x[order_x][i&1];
src_coordinates_y_unstruct[i] =
src_coordinates_y[at_pole][order_y][i>>1];
}
grid_data =
4, 1, (int[]){4},
src_coordinates_x_unstruct,
src_coordinates_y_unstruct, (int[]){0,1,3,2});
}
}
if (!is_tgt) {
}
if (with_mask && !is_tgt) {
mask[i] =
src_cell_mask[src_reorder[order_x][order_y][i]];
}
.masks_idx = (with_mask)?0:SIZE_MAX}};
size_t num_src_fields =
sizeof(src_fields) / sizeof(src_fields[0]);
.masks_idx = SIZE_MAX};
for (size_t i = 0; i < num_reorder_types; ++i) {
weights, reorder_types[i], 1,
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
double * ref_tgt_field = NULL;
double ref_global_tgt_field[25];
double ref_weights[2][2][25][4] =
{{{{1.00, 0.00, 0.00, 0.00},
{0.75, 0.25, 0.00, 0.00},
{0.50, 0.50, 0.00, 0.00},
{0.25, 0.75, 0.00, 0.00},
{0.00, 1.00, 0.00, 0.00},
{0.75, 0.00, 0.25, 0.00},
{0.75, 0.00, 0.00, 0.25},
{0.50, 0.25, 0.00, 0.25},
{0.25, 0.50, 0.00, 0.25},
{0.00, 0.75, 0.00, 0.25},
{0.50, 0.00, 0.50, 0.00},
{0.50, 0.00, 0.25, 0.25},
{0.50, 0.00, 0.00, 0.50},
{0.25, 0.25, 0.00, 0.50},
{0.00, 0.50, 0.00, 0.50},
{0.25, 0.00, 0.75, 0.00},
{0.25, 0.00, 0.50, 0.25},
{0.25, 0.00, 0.25, 0.50},
{0.25, 0.00, 0.00, 0.75},
{0.00, 0.25, 0.00, 0.75},
{0.00, 0.00, 1.00, 0.00},
{0.00, 0.00, 0.75, 0.25},
{0.00, 0.00, 0.50, 0.50},
{0.00, 0.00, 0.25, 0.75},
{0.00, 0.00, 0.00, 1.00}},
{{1.00, 0.00, 0.00, 0.00},
{0.75, 0.25, 0.00, 0.00},
{0.50, 0.50, 0.00, 0.00},
{0.25, 0.75, 0.00, 0.00},
{0.00, 1.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.75, 0.00, 0.00, 0.25},
{0.50, 0.25, 0.00, 0.25},
{0.25, 0.50, 0.00, 0.25},
{0.00, 0.75, 0.00, 0.25},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.50, 0.00, 0.00, 0.50},
{0.25, 0.25, 0.00, 0.50},
{0.00, 0.50, 0.00, 0.50},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.25, 0.00, 0.00, 0.75},
{0.00, 0.25, 0.00, 0.75},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 1.00}}},
{{{1.00, 0.00, 0.00, 0.00},
{0.75, 0.25, 0.00, 0.00},
{0.50, 0.50, 0.00, 0.00},
{0.25, 0.75, 0.00, 0.00},
{0.00, 1.00, 0.00, 0.00},
{0.75, 0.00, 0.25, 0.00},
{0.75, 0.00, 0.00, 0.25},
{0.50, 0.25, 0.00, 0.25},
{0.25, 0.50, 0.00, 0.25},
{0.00, 0.75, 0.00, 0.25},
{0.50, 0.00, 0.50, 0.00},
{0.50, 0.00, 0.25, 0.25},
{0.50, 0.00, 0.00, 0.50},
{0.25, 0.25, 0.00, 0.50},
{0.00, 0.50, 0.00, 0.50},
{0.25, 0.00, 0.75, 0.00},
{0.25, 0.00, 0.50, 0.25},
{0.25, 0.00, 0.25, 0.50},
{0.25, 0.00, 0.00, 0.75},
{0.00, 0.25, 0.00, 0.75},
{0.00, 0.00, 1.00, 0.00},
{0.00, 0.00, 0.75, 0.25},
{0.00, 0.00, 0.50, 0.50},
{0.00, 0.00, 0.25, 0.75},
{0.00, 0.00, 0.00, 1.00}},
{{1.00, 0.00, 0.00, 0.00},
{0.75, 0.25, 0.00, 0.00},
{0.50, 0.50, 0.00, 0.00},
{0.25, 0.75, 0.00, 0.00},
{0.00, 1.00, 0.00, 0.00},
{1.00, 0.00, 0.00, 0.00},
{0.75, 0.00, 0.00, 0.25},
{0.50, 0.25, 0.00, 0.25},
{0.25, 0.50, 0.00, 0.25},
{0.00, 0.75, 0.00, 0.25},
{1.00, 0.00, 0.00, 0.00},
{0.50/0.75, 0.00, 0.00, 0.25/0.75},
{0.50, 0.00, 0.00, 0.50},
{0.25, 0.25, 0.00, 0.50},
{0.00, 0.50, 0.00, 0.50},
{1.00, 0.00, 0.00, 0.00},
{0.50, 0.00, 0.00, 0.50},
{0.25/0.75, 0.00, 0.00, 0.50/0.75},
{0.25, 0.00, 0.00, 0.75},
{0.00, 0.25, 0.00, 0.75},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 1.00},
{0.00, 0.00, 0.00, 1.00},
{0.00, 0.00, 0.00, 1.00},
{0.00, 0.00, 0.00, 1.00}}}};
for (size_t j = 0; j < 25; ++j) {
ref_global_tgt_field[j] =
(double)(src_reorder[order_x][order_y][0] + 1) *
ref_weights[partial_coverage][with_mask][j][0] +
(double)(src_reorder[order_x][order_y][1] + 1) *
ref_weights[partial_coverage][with_mask][j][1] +
(double)(src_reorder[order_x][order_y][2] + 1) *
ref_weights[partial_coverage][with_mask][j][2] +
(double)(src_reorder[order_x][order_y][3] + 1) *
ref_weights[partial_coverage][with_mask][j][3];
if (ref_global_tgt_field[j] == 0.0)
ref_global_tgt_field[j] = -1.0;
}
if (is_tgt) {
ref_tgt_field[i] = ref_global_tgt_field[i];
} else {
src_field[i] = (double)(i) + 1.0;
}
if (is_tgt)
if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-3)
PUT_ERR("wrong interpolation result");
free(src_field);
free(tgt_field);
free(ref_tgt_field);
}
}
}
}
}
}
}
}
}
{
int is_tgt = split_comm_size == 1;
size_t src_num_vertices = 4;
size_t src_num_cells = 2;
int src_num_vertices_per_cell[] = {3, 3};
int src_cell_to_vertex[] = {0,1,3, 0,2,3};
double src_coordinates_x[4] = {-0.5,0.5,-0.5,0.5};
double src_coordinates_y[4] = {-0.5,-0.5,0.5,0.5};
size_t tgt_num_vertices[2] = {5,5};
int tgt_cyclic[2] = {0,0};
double tgt_coordinates_x[5] = {-0.5,-0.25,0.0,0.25,0.5};
double tgt_coordinates_y[5] = {-0.5,-0.25,0.0,0.25,0.5};
int src_cell_mask[4] = {1,1, 0,1};
for (int with_mask = 0; with_mask < 2; ++with_mask) {
if (is_tgt) {
grid_data =
tgt_num_vertices, tgt_cyclic,
tgt_coordinates_x, tgt_coordinates_y);
} else {
grid_data =
src_num_vertices, src_num_cells, src_num_vertices_per_cell,
src_coordinates_x, src_coordinates_y, src_cell_to_vertex);
}
if (with_mask && !is_tgt) {
mask[i] = src_cell_mask[i];
}
.masks_idx = (with_mask)?0:SIZE_MAX}};
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
.masks_idx = SIZE_MAX};
for (int partial_coverage = 0; partial_coverage < 2; ++partial_coverage) {
for (size_t i = 0; i < num_reorder_types; ++i) {
weights, reorder_types[i], 1,
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
double * ref_tgt_field = NULL;
double ref_global_tgt_field[25];
double ref_weights[2][2][25][4] =
{{{{1.00, 0.00, 0.00, 0.00},
{0.75, 0.25, 0.00, 0.00},
{0.50, 0.50, 0.00, 0.00},
{0.25, 0.75, 0.00, 0.00},
{0.00, 1.00, 0.00, 0.00},
{0.75, 0.00, 0.25, 0.00},
{0.75, 0.00, 0.00, 0.25},
{0.50, 0.25, 0.00, 0.25},
{0.25, 0.50, 0.00, 0.25},
{0.00, 0.75, 0.00, 0.25},
{0.50, 0.00, 0.50, 0.00},
{0.50, 0.00, 0.25, 0.25},
{0.50, 0.00, 0.00, 0.50},
{0.25, 0.25, 0.00, 0.50},
{0.00, 0.50, 0.00, 0.50},
{0.25, 0.00, 0.75, 0.00},
{0.25, 0.00, 0.50, 0.25},
{0.25, 0.00, 0.25, 0.50},
{0.25, 0.00, 0.00, 0.75},
{0.00, 0.25, 0.00, 0.75},
{0.00, 0.00, 1.00, 0.00},
{0.00, 0.00, 0.75, 0.25},
{0.00, 0.00, 0.50, 0.50},
{0.00, 0.00, 0.25, 0.75},
{0.00, 0.00, 0.00, 1.00}},
{{1.00, 0.00, 0.00, 0.00},
{0.75, 0.25, 0.00, 0.00},
{0.50, 0.50, 0.00, 0.00},
{0.25, 0.75, 0.00, 0.00},
{0.00, 1.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.75, 0.00, 0.00, 0.25},
{0.50, 0.25, 0.00, 0.25},
{0.25, 0.50, 0.00, 0.25},
{0.00, 0.75, 0.00, 0.25},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.50, 0.00, 0.00, 0.50},
{0.25, 0.25, 0.00, 0.50},
{0.00, 0.50, 0.00, 0.50},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.25, 0.00, 0.00, 0.75},
{0.00, 0.25, 0.00, 0.75},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 1.00}}},
{{{1.00, 0.00, 0.00, 0.00},
{0.75, 0.25, 0.00, 0.00},
{0.50, 0.50, 0.00, 0.00},
{0.25, 0.75, 0.00, 0.00},
{0.00, 1.00, 0.00, 0.00},
{0.75, 0.00, 0.25, 0.00},
{0.75, 0.00, 0.00, 0.25},
{0.50, 0.25, 0.00, 0.25},
{0.25, 0.50, 0.00, 0.25},
{0.00, 0.75, 0.00, 0.25},
{0.50, 0.00, 0.50, 0.00},
{0.50, 0.00, 0.25, 0.25},
{0.50, 0.00, 0.00, 0.50},
{0.25, 0.25, 0.00, 0.50},
{0.00, 0.50, 0.00, 0.50},
{0.25, 0.00, 0.75, 0.00},
{0.25, 0.00, 0.50, 0.25},
{0.25, 0.00, 0.25, 0.50},
{0.25, 0.00, 0.00, 0.75},
{0.00, 0.25, 0.00, 0.75},
{0.00, 0.00, 1.00, 0.00},
{0.00, 0.00, 0.75, 0.25},
{0.00, 0.00, 0.50, 0.50},
{0.00, 0.00, 0.25, 0.75},
{0.00, 0.00, 0.00, 1.00}},
{{1.00, 0.00, 0.00, 0.00},
{0.75, 0.25, 0.00, 0.00},
{0.50, 0.50, 0.00, 0.00},
{0.25, 0.75, 0.00, 0.00},
{0.00, 1.00, 0.00, 0.00},
{1.00, 0.00, 0.00, 0.00},
{0.75, 0.00, 0.00, 0.25},
{0.50, 0.25, 0.00, 0.25},
{0.25, 0.50, 0.00, 0.25},
{0.00, 0.75, 0.00, 0.25},
{1.00, 0.00, 0.00, 0.00},
{0.50/0.75, 0.00, 0.00, 0.25/0.75},
{0.50, 0.00, 0.00, 0.50},
{0.25, 0.25, 0.00, 0.50},
{0.00, 0.50, 0.00, 0.50},
{1.00, 0.00, 0.00, 0.00},
{0.50, 0.00, 0.00, 0.50},
{0.25/0.75, 0.00, 0.00, 0.50/0.75},
{0.25, 0.00, 0.00, 0.75},
{0.00, 0.25, 0.00, 0.75},
{0.00, 0.00, 0.00, 0.00},
{0.00, 0.00, 0.00, 1.00},
{0.00, 0.00, 0.00, 1.00},
{0.00, 0.00, 0.00, 1.00},
{0.00, 0.00, 0.00, 1.00}}}};
for (size_t j = 0; j < 25; ++j) {
ref_global_tgt_field[j] =
1.0 * ref_weights[partial_coverage][with_mask][j][0] +
2.0 * ref_weights[partial_coverage][with_mask][j][1] +
3.0 * ref_weights[partial_coverage][with_mask][j][2] +
4.0 * ref_weights[partial_coverage][with_mask][j][3];
if (ref_global_tgt_field[j] == 0.0)
ref_global_tgt_field[j] = -1.0;
}
if (is_tgt) {
ref_tgt_field[i] = ref_global_tgt_field[i];
} else {
src_field[i] = (double)(i) + 1.0;
}
if (is_tgt)
if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-3)
PUT_ERR("wrong interpolation result");
free(src_field);
free(tgt_field);
free(ref_tgt_field);
}
}
}
}
}
{
int is_tgt = split_comm_size == 1;
double coordinates_x[2][7] = {{0.0,0.1,0.2}, {0.0,0.05,0.1,0.15,0.2}};
double coordinates_y[2][7] = {{0.0,0.1,0.2}, {0.0,0.05,0.1,0.15,0.2}};
size_t const num_cells[2][2] = {{2,2}, {4,4}};
size_t local_start[2][2][2] = {{{0,0},{0,1}}, {{0,0}}};
size_t local_count[2][2][2] = {{{2,1},{2,1}}, {{4,4}}};
int global_mask[2][5*5] = {{1,1,1,
0,1,1,
0,0,1},
{1,1,1,1,1,
1,1,1,1,0,
1,1,1,0,0,
1,1,0,0,0,
1,0,0,0,0}};
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;
yac_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);
int field_mask[5*5];
field_mask[i] = global_mask[is_tgt][grid_data.
vertex_ids[i]];
for (int with_src_mask = 0; with_src_mask < 2; ++with_src_mask) {
for (int with_tgt_mask = 0; with_tgt_mask < 2; ++with_tgt_mask) {
.coordinates_idx = SIZE_MAX,
.masks_idx = (with_src_mask)?0:SIZE_MAX}};
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
.coordinates_idx = SIZE_MAX,
.masks_idx = (with_tgt_mask)?0:SIZE_MAX};
size_t const num_weight_types =
sizeof(weight_types) / sizeof(weight_types[0]);
for (size_t weight_types_idx = 0; weight_types_idx < num_weight_types;
++weight_types_idx) {
for (int partial_coverage = 0; partial_coverage < 2;
++partial_coverage) {
weight_types[weight_types_idx], partial_coverage),
for (size_t i = 0; i < num_reorder_types; ++i) {
weights, reorder_types[i], 1,
{
double src_field_data[3*3];
double * src_field = src_field_data;
double ** src_fields = &src_field;
double tgt_field_data[5*5] =
{-2.0,-2.0,-2.0,-2.0,-2.0,
-2.0,-2.0,-2.0,-2.0,-2.0,
-2.0,-2.0,-2.0,-2.0,-2.0,
-2.0,-2.0,-2.0,-2.0,-2.0,
-2.0,-2.0,-2.0,-2.0,-2.0};
double * tgt_field = tgt_field_data;
if (!is_tgt)
interpolation, &src_fields, &tgt_field);
if (is_tgt) {
size_t src_per_tgt[5*5][4] =
{{0,1,3,4},{0,1,3,4},{0,1,3,4},{1,2,4,5},{1,2,4,5},
{0,1,3,4},{0,1,3,4},{0,1,3,4},{1,2,4,5},{1,2,4,5},
{0,1,3,4},{0,1,3,4},{0,1,3,4},{1,2,4,5},{1,2,4,5},
{3,4,6,7},{3,4,6,7},{3,4,6,7},{4,5,7,8},{4,5,7,8},
{3,4,6,7},{3,4,6,7},{3,4,6,7},{4,5,7,8},{4,5,7,8}};
double A = 0.345491502812526;
double B = 0.154508497187473;
double C = 0.25;
double D = 0.5;
double weights[3][5*5][4] =
{{{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},
{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},
{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},
{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},
{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C}},
{{1,0,0,0},{A,A,B,B},{0,1,0,0},{A,A,B,B},{0,1,0,0},
{A,B,A,B},{C,C,C,C},{B,A,B,A},{C,C,C,C},{B,A,B,A},
{0,0,1,0},{B,B,A,A},{0,0,0,1},{B,B,A,A},{0,0,0,1},
{A,B,A,B},{C,C,C,C},{B,A,B,A},{C,C,C,C},{B,A,B,A},
{0,0,1,0},{B,B,A,A},{0,0,0,1},{B,B,A,A},{0,0,0,1}},
{{1,0,0,0},{D,D,0,0},{0,1,0,0},{D,D,0,0},{0,1,0,0},
{D,0,D,0},{D,0,0,D},{0,D,0,D},{D,0,0,D},{0,D,0,D},
{0,0,1,0},{0,0,D,D},{0,0,0,1},{0,0,D,D},{0,0,0,1},
{D,0,D,0},{D,0,0,D},{0,D,0,D},{D,0,0,D},{0,D,0,D},
{0,0,1,0},{0,0,D,D},{0,0,0,1},{0,0,D,D},{0,0,0,1}}};
double src_values[4];
double w[4];
double ref_tgt_value = 0.0;
int contains_masked_src = 0;
int contains_unmasked_src = 0;
for (int j = 0; j < 4; ++j) {
size_t src_idx = src_per_tgt[tgt_idx][j];
w[j] = weights[weight_types_idx][tgt_idx][j];
if (with_src_mask &&
!global_mask[0][src_idx] &&
(w[j] != 0.0)) {
src_values[j] = -1.0;
w[j] = 0.0;
contains_masked_src = 1;
} else {
src_values[j] = (double)src_idx;
contains_unmasked_src |= w[j] != 0.0;
}
}
if (with_tgt_mask && !global_mask[1][tgt_idx]) {
ref_tgt_value = -2.0;
} else if (contains_masked_src) {
if (partial_coverage && contains_unmasked_src) {
for (int j = 0; j < 4; ++j)
ref_tgt_value += w[j] * src_values[j];
} else {
ref_tgt_value = -1.0;
}
} else {
for (int j = 0; j < 4; ++j)
ref_tgt_value += w[j] * src_values[j];
}
if (with_src_mask && partial_coverage &&
(weight_types[weight_types_idx] ==
(!with_tgt_mask || global_mask[1][tgt_idx])) {
switch (tgt_idx) {
case(10):
ref_tgt_value = 1.73879612503;
break;
case(20):
case(22):
ref_tgt_value = 4.0;
default:
break;
}
}
if (fabs(ref_tgt_value - tgt_field[i]) > 1e-3)
PUT_ERR("wrong interpolation result");
}
}
}
}
}
}
}
}
}
MPI_Comm_free(&split_comm);
xt_finalize();
MPI_Finalize();
return TEST_EXIT_CODE;
}
struct yac_basic_grid * yac_basic_grid_new(char const *name, struct yac_basic_grid_data grid_data)
size_t yac_basic_grid_add_mask(struct yac_basic_grid *grid, enum yac_location location, int const *mask, size_t count, char const *mask_name)
size_t yac_basic_grid_add_coordinates_nocpy(struct yac_basic_grid *grid, enum yac_location location, yac_coordinate_pointer coordinates)
struct yac_basic_grid * yac_basic_grid_empty_new(char const *name)
void yac_basic_grid_delete(struct yac_basic_grid *grid)
size_t yac_basic_grid_add_mask_nocpy(struct yac_basic_grid *grid, enum yac_location location, int const *mask, char const *mask_name)
struct yac_basic_grid_data yac_generate_basic_grid_data_unstruct_deg(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell, double *x_vertices, double *y_vertices, int *cell_to_vertex)
struct yac_basic_grid_data yac_generate_basic_grid_data_unstruct(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell, double *x_vertices, double *y_vertices, int *cell_to_vertex)
struct yac_basic_grid_data yac_generate_basic_grid_data_reg_2d_deg(size_t nbr_vertices[2], int cyclic[2], double *lon_vertices, double *lat_vertices)
void yac_correct_weights(size_t nSourceCells, double *weight)
correct interpolation weights
int main(int argc, char **argv)
void yac_dist_grid_pair_delete(struct yac_dist_grid_pair *grid_pair)
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)
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
static void normalise_vector(double v[])
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)
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_avg_new(enum yac_interp_avg_weight_type weight_type, int partial_coverage)
yac_interp_avg_weight_type
@ YAC_INTERP_AVG_ARITHMETIC
struct interp_method * yac_interp_method_fixed_new(double value)
void yac_interp_weights_delete(struct yac_interp_weights *weights)
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, char const *yaxt_exchanger_name)
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_coordinate_pointer vertex_coordinates
size_t * cell_to_vertex_offsets
int * num_vertices_per_cell
enum yac_location location
struct yac_interp_field tgt_field
struct yac_dist_grid_pair * grid_pair
struct yac_interp_field src_fields[]
struct yac_basic_grid ** grids
double(* yac_coordinate_pointer)[3]