A test for the parallel source point mapping 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_weight_types = sizeof(weight_types) / sizeof(weight_types[0]);
};
size_t num_scale_types = sizeof(scale_types) / sizeof(scale_types[0]);
};
size_t num_reorder_types = sizeof(reorder_types) / sizeof(reorder_types[0]);
static char const * grid_names[2] = {"grid1", "grid2"};
static double compute_scale(
double src_cell_area, double tgt_cell_area);
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 != 6) {
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 < 1, 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[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][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{7,7}}};
int global_mask[7*7] = {
0,0,1,1,1,1,1,
0,1,1,2,2,2,2,
1,1,2,2,2,2,2,
1,2,2,1,1,1,2,
1,2,2,1,0,1,2,
1,2,2,1,1,1,2,
1,2,2,2,2,2,2};
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);
{
int valid_mask_value = (is_tgt)?2:1;
for (
size_t i = 0; i < grid_data.
num_cells; ++i) {
double * middle_point = 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];
}
mask[i] = global_mask[grid_data.
cell_ids[i]] == valid_mask_value;
}
}
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
for (size_t i = 0; i < num_reorder_types; ++i) {
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
if (is_tgt) {
for (
size_t i = 0; i < grid_data.
num_cells; ++i) tgt_field[i] = -1;
} else {
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
src_field[i] = (
double)(grid_data.
cell_ids[i]);
}
double ref_tgt_field[7*7] = {
-1,-1,-1,-1,-1,-1,-1,
-1,-1,-1,2+3+9,4,5,6,
-1,-1,8+15,-2,25,-2,-2,
-1,14+21,24,-1,-1,-1,26,
-1,28,31,-1,-1,-1,33,
-1,35,38,-1,-1,-1,40,
-1,42,-2,-2,39,-2,-2};
if (is_tgt)
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
if (ref_tgt_field[grid_data.
cell_ids[i]] != tgt_field[i])
PUT_ERR("wrong interpolation result");
free(tgt_field);
free(src_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][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{7,7}}};
int global_mask[7*7] = {
0,0,1,1,1,1,1,
0,1,1,2,2,2,2,
1,1,2,2,2,2,2,
1,2,2,1,1,1,2,
1,2,2,1,0,1,2,
1,2,2,1,1,1,2,
1,2,2,2,2,2,2};
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);
{
int valid_mask_value = (is_tgt)?2:1;
for (
size_t i = 0; i < grid_data.
num_cells; ++i) {
double * middle_point = 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];
}
mask[i] = global_mask[grid_data.
cell_ids[i]] == valid_mask_value;
}
}
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
for (size_t i = 0; i < num_reorder_types; ++i) {
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
if (is_tgt) {
for (
size_t i = 0; i < grid_data.
num_cells; ++i) tgt_field[i] = -1;
} else {
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
src_field[i] = (
double)(grid_data.
cell_ids[i]);
}
double ref_tgt_field[7*7] = {
-1,-1,-1,-1,-1,-1,-1,
-1,-1,-1, 0, 0, 0, 0,
-1,-1, 0, 0, 0, 0, 0,
-1, 0, 0,-1,-1,-1, 0,
-1, 0, 0,-1,-1,-1, 0,
-1, 0, 0,-1,-1,-1, 0,
-1, 0, 0, 0, 0, 0, 0};
size_t coast_point[] = {2,3,4,5,6,
8,9,
14,15,
21,24,25,26,
28,31,33,
35,38,39,40,
42};
size_t num_coast_points = sizeof(coast_point)/sizeof(coast_point[0]);
size_t num_tgt_per_coast[] = {3,3,4,4,3,
3,3,
3,3,
3,4,4,3,
4,4,3,
4,4,3,3,
3};
size_t tgts[] = {10,11,17, 10,11,17, 10,11,12,18, 11,12,13,19, 12,13,20,
16,17,23, 10,11,17,
22,23,29, 16,17,23,
22,23,29, 16,22,23,30, 11,17,18,19, 20,27,34,
22,29,30,36, 23,29,30,37, 27,34,41,
29,36,37,43, 30,36,37,44, 45,46,47, 34,41,48,
36,43,44};
for (size_t i = 0, k = 0; i < num_coast_points; ++i) {
size_t curr_num_tgt = num_tgt_per_coast[i];
double curr_data = (double)(coast_point[i]) / (double)curr_num_tgt;
for (size_t j = 0; j < curr_num_tgt; ++j, ++k)
ref_tgt_field[tgts[k]] += curr_data;
}
if (is_tgt)
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
if (fabs(ref_tgt_field[grid_data.
cell_ids[i]] -
tgt_field[i]) > 1e-6)
PUT_ERR("wrong interpolation result");
free(tgt_field);
free(src_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][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{7,7}}};
int global_mask[7*7] = {
0,0,0,0,0,0,1,
1,0,0,0,0,1,2,
2,1,0,0,1,2,2,
2,2,1,1,2,2,2,
2,2,1,1,2,2,2,
2,1,0,0,1,2,2,
1,0,0,0,0,1,2};
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);
{
int valid_mask_value = (is_tgt)?2:1;
for (
size_t i = 0; i < grid_data.
num_cells; ++i) {
double * middle_point = 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];
}
mask[i] = global_mask[grid_data.
cell_ids[i]] == valid_mask_value;
}
}
double grid_cell_areas[7*7];
for (int i = 0, k = 0; i < 7; ++i)
for (int j = 0; j < 7; ++j, ++k)
grid_cell_areas[k] =
fabs(
(coordinates_x[j+1] - coordinates_x[j+0]) *
(sin(coordinates_y[i+0]) - sin(coordinates_y[i+1])));
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
double const radius[] = {1.0, 2.0};
enum {NUM_RADIUS = sizeof(radius) / sizeof(radius[0])};
for (size_t weight_type_idx = 0; weight_type_idx < num_weight_types;
++weight_type_idx) {
for (size_t scale_type_idx = 0; scale_type_idx < num_scale_types;
++scale_type_idx) {
for (size_t src_sphere_radius_idx = 0;
src_sphere_radius_idx < NUM_RADIUS; ++src_sphere_radius_idx) {
for (size_t tgt_sphere_radius_idx = 0;
tgt_sphere_radius_idx < NUM_RADIUS; ++tgt_sphere_radius_idx) {
weight_types[weight_type_idx], scale_types[scale_type_idx],
radius[src_sphere_radius_idx], radius[tgt_sphere_radius_idx]),
for (size_t i = 0; i < num_reorder_types; ++i) {
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
if (is_tgt) {
for (
size_t i = 0; i < grid_data.
num_cells; ++i) tgt_field[i] = -1;
} else {
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
src_field[i] = (
double)(grid_data.
cell_ids[i]);
}
if (is_tgt) {
double ref_tgt_field[7*7] = {
-1,-1,-1,-1,-1,-1,-1,
-1,-1,-1,-1,-1,-1, 0,
0,-1,-1,-1,-1, 0, 0,
0, 0,-1,-1, 0, 0, 0,
0, 0,-1,-1, 0, 0, 0,
0,-1,-1,-1,-1, 0, 0,
-1,-1,-1,-1,-1,-1, 0};
size_t coast_point[] = {7,15,23,30,36,42,
6,12,18,24,31,39,47};
size_t num_coast_points = sizeof(coast_point)/sizeof(coast_point[0]);
size_t num_tgt_per_coast[] = {6,6,6,6,6,6,
9,9,11,12,12,11,9};
size_t tgts[] = {14,21,22,28,29,35,
14,21,22,28,29,35,
14,21,22,28,29,35,
14,21,22,28,29,35,
14,21,22,28,29,35,
14,21,22,28,29,35,
13,19,20,25,26,27,32,33,34,
13,19,20,25,26,27,32,33,34,
13,19,20,25,26,27,32,33,34,40,41,
13,19,20,25,26,27,32,33,34,40,41,48,
13,19,20,25,26,27,32,33,34,40,41,48,
19,20,25,26,27,32,33,34,40,41,48,
25,26,27,32,33,34,40,41,48};
switch (weight_types[weight_type_idx]) {
default:
for (size_t i = 0, k = 0; i < num_coast_points; ++i) {
size_t curr_num_tgt = num_tgt_per_coast[i];
double curr_data = (double)(coast_point[i]) / (double)curr_num_tgt;
for (size_t j = 0; j < curr_num_tgt; ++j, ++k)
ref_tgt_field[tgts[k]] +=
curr_data *
compute_scale(
scale_types[scale_type_idx],
grid_cell_areas[coast_point[i]] *
radius[src_sphere_radius_idx] *
radius[src_sphere_radius_idx],
grid_cell_areas[tgts[k]] *
radius[tgt_sphere_radius_idx] *
radius[tgt_sphere_radius_idx]);
}
break;
}
size_t * curr_tgts = tgts;
for (size_t i = 0, k = 0; i < num_coast_points; ++i) {
size_t curr_num_tgt = num_tgt_per_coast[i];
size_t curr_coast_point = coast_point[i];
double curr_src_data = (double)(curr_coast_point);
double inv_distances[curr_num_tgt];
double inv_distances_sum = 0.0;
for (size_t j = 0; j < curr_num_tgt; ++j)
inv_distances_sum +=
((inv_distances[j] =
point_coordinates[curr_coast_point],
point_coordinates[curr_tgts[j]])));
for (size_t j = 0; j < curr_num_tgt; ++j, ++k)
ref_tgt_field[curr_tgts[j]] +=
curr_src_data
* (inv_distances[j] / inv_distances_sum) *
compute_scale(
scale_types[scale_type_idx],
grid_cell_areas[coast_point[i]] *
radius[src_sphere_radius_idx] *
radius[src_sphere_radius_idx],
grid_cell_areas[tgts[k]] *
radius[tgt_sphere_radius_idx] *
radius[tgt_sphere_radius_idx]);
curr_tgts += curr_num_tgt;
}
}
}
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
if (fabs(ref_tgt_field[grid_data.
cell_ids[i]] -
tgt_field[i]) > 1.0e-6)
PUT_ERR("wrong interpolation result");
double src_sum = 0.0;
for (size_t i = 0; i < num_coast_points; ++i) {
double curr_src_data = (double)(coast_point[i]);
if (
(scale_types[scale_type_idx] ==
(scale_types[scale_type_idx] ==
curr_src_data *=
grid_cell_areas[coast_point[i]] *
radius[src_sphere_radius_idx] *
radius[src_sphere_radius_idx];
src_sum += curr_src_data;
}
double tgt_sum = 0.0;
for (
size_t i = 0; i < grid_data.
num_cells; ++i) {
if (tgt_field[i] != -1.0) {
double curr_tgt_data = tgt_field[i];
if (
(scale_types[scale_type_idx] ==
(scale_types[scale_type_idx] ==
curr_tgt_data *=
grid_cell_areas[grid_data.
cell_ids[i]] *
radius[tgt_sphere_radius_idx] *
radius[tgt_sphere_radius_idx];
tgt_sum += curr_tgt_data;
}
}
if (fabs(src_sum - tgt_sum) > 1.0e-6)
PUT_ERR("wrong interpolation result (not conservative)");
}
free(tgt_field);
free(src_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][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{7,7}}};
int global_mask[7*7] = {
1,0,0,0,0,0,0,
0,1,0,0,0,0,0,
0,0,1,0,0,0,0,
0,0,0,1,0,0,1,
0,0,0,0,1,1,2,
0,0,0,0,1,2,2,
0,0,0,1,2,2,2};
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);
{
int valid_mask_value = (is_tgt)?2:1;
for (
size_t i = 0; i < grid_data.
num_cells; ++i) {
double * middle_point = 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];
}
mask[i] = global_mask[grid_data.
cell_ids[i]] == valid_mask_value;
}
}
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
for (size_t i = 0; i < num_reorder_types; ++i) {
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
if (is_tgt) {
for (
size_t i = 0; i < grid_data.
num_cells; ++i) tgt_field[i] = -1;
} else {
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
src_field[i] = (
double)(grid_data.
cell_ids[i]);
}
if (is_tgt) {
double ref_tgt_field[7*7] = {
-1,-1,-1,-1,-1,-1,-1,
-1,-1,-1,-1,-1,-1,-1,
-1,-1,-1,-1,-1,-1,-1,
-1,-1,-1,-1,-1,-1,-1,
-1,-1,-1,-1,-1,-1, 0,
-1,-1,-1,-1,-1, 0,-2,
-1,-1,-1,-1, 0,-2,-2};
size_t coast_point[] = {0,8,16,24,27,32,33,39,45};
size_t num_coast_points = sizeof(coast_point)/sizeof(coast_point[0]);
size_t num_tgt_per_coast[] = {0,0,0,1,1,1,1,1,1};
size_t tgts[] = {40,34,40,34,40,46};
for (size_t i = 0, k = 0; i < num_coast_points; ++i) {
size_t curr_num_tgt = num_tgt_per_coast[i];
if (curr_num_tgt == 0) continue;
double curr_data = (double)(coast_point[i]) / (double)curr_num_tgt;
for (size_t j = 0; j < curr_num_tgt; ++j, ++k)
ref_tgt_field[tgts[k]] += curr_data;
}
for (
size_t i = 0; i < grid_data.
num_cells; ++i)
if (fabs(ref_tgt_field[grid_data.
cell_ids[i]] -
tgt_field[i]) > 1e-6)
PUT_ERR("wrong interpolation result");
}
free(tgt_field);
free(src_field);
}
}
}
MPI_Comm_free(&split_comm);
xt_finalize();
MPI_Finalize();
return TEST_EXIT_CODE;
}
static double compute_scale(
double src_cell_area, double tgt_cell_area) {
double scale = 1.0;
scale *= src_cell_area;
scale /= tgt_cell_area;
return scale;
}
struct yac_basic_grid * yac_basic_grid_new(char const *name, struct yac_basic_grid_data grid_data)
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)
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 normalise_vector(double v[])
static double get_vector_angle(double const a[3], double const b[3])
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_fixed_new(double value)
struct interp_method * yac_interp_method_spmap_new(double spread_distance, double max_search_distance, enum yac_interp_spmap_weight_type weight_type, enum yac_interp_spmap_scale_type scale_type, double src_sphere_radius, double tgt_sphere_radius)
#define YAC_INTERP_SPMAP_MAX_SEARCH_DISTANCE_DEFAULT
#define YAC_INTERP_SPMAP_SCALE_DEFAULT
yac_interp_spmap_scale_type
@ YAC_INTERP_SPMAP_NONE
weights are not scaled
@ YAC_INTERP_SPMAP_INVTGTAREA
@ YAC_INTERP_SPMAP_SRCAREA
@ YAC_INTERP_SPMAP_FRACAREA
#define YAC_INTERP_SPMAP_TGT_SPHERE_RADIUS_DEFAULT
#define YAC_INTERP_SPMAP_SPREAD_DISTANCE_DEFAULT
yac_interp_spmap_weight_type
#define YAC_INTERP_SPMAP_SRC_SPHERE_RADIUS_DEFAULT
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_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[]
int main(int argc, char **argv)
struct yac_basic_grid ** grids
double(* yac_coordinate_pointer)[3]
