YAC 3.12.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_interp_method_avg_parallel.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#include <stdlib.h>
6#include <string.h>
7
8#include "tests.h"
9#include "test_common.h"
10#include "geometry.h"
11#include "interpolation/methods/interp_method.h"
12#include "interpolation/methods/interp_method_avg.h"
13#include "interpolation/methods/interp_method_fixed.h"
14#include "dist_grid_utils.h"
15#include "clipping.h"
16#include "yac_mpi.h"
17
18#include <mpi.h>
19#include <yaxt.h>
20#include <netcdf.h>
21
29size_t num_reorder_types = sizeof(reorder_types) / sizeof(reorder_types[0]);
30static char const * grid_names[2] = {"src_grid", "tgt_grid"};
31
32int main(void) {
33
34 MPI_Init(NULL, NULL);
35
36 xt_initialize(MPI_COMM_WORLD);
37
38 int comm_rank, comm_size;
39 MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
40 MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
41 MPI_Barrier(MPI_COMM_WORLD);
42
43 if (comm_size != 3) {
44 PUT_ERR("ERROR: wrong number of processes");
45 xt_finalize();
46 MPI_Finalize();
47 return TEST_EXIT_CODE;
48 }
49
50 MPI_Comm split_comm;
51 MPI_Comm_split(
52 MPI_COMM_WORLD, comm_rank < 2, 0, &split_comm);
53
54 int split_comm_rank, split_comm_size;
55 MPI_Comm_rank(split_comm, &split_comm_rank);
56 MPI_Comm_size(split_comm, &split_comm_size);
57
58 {// Test 1a
59 // linear point interpolation (fixed is backup) with two source processes and a
60 // single target process the source fields have two struct points per grid
61 // (the test checks how the interpolation methods handle NULL target coordinates )
62
63 // the global grid is a 4x2 grid:
64 // 08--14--09--15--10--16--11
65 // | | | |
66 // 08 03 10 04 12 05 13
67 // | | | |
68 // 04--07--05--09--06--11--07
69 // | | | |
70 // 01 00 03 01 05 02 06
71 // | | | |
72 // 00--00--01--02--02--04--03
73 //
74 //---------------
75 // setup
76 //---------------
77
78 int is_tgt = split_comm_size == 1;
79 double coordinates_x[2][4] = {{0.0,1.0,2.0,3.0}, {0.5,1.5,2.5}};
80 double coordinates_y[2][3] = {{0.0,1.0,2.0}, {0.5,1.5,2.5}};
81 size_t const num_cells[2][2] = {{3,2}, {2,2}};
82 size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
83 size_t local_count[2][2][2] = {{{2,2},{2,2}}, {{2,2}}};
84 int with_halo = 0;
85 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
86 coordinates_x[is_tgt][i] *= YAC_RAD;
87 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
88 coordinates_y[is_tgt][i] *= YAC_RAD;
89
90 struct yac_basic_grid_data grid_data =
91 yac_generate_basic_grid_data_reg2d(
92 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
93 local_start[is_tgt][split_comm_rank],
94 local_count[is_tgt][split_comm_rank], with_halo);
95
96 struct yac_basic_grid * grids[2] =
97 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
98 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
99
100 struct yac_dist_grid_pair * grid_pair =
101 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
102
103 struct yac_interp_field src_fields[] =
104 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
105 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
106 struct yac_interp_field tgt_field =
107 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
108
109 struct yac_interp_grid * interp_grid =
110 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
111 num_src_fields, src_fields, tgt_field);
112
113 for (int partial_coverage = 0; partial_coverage <= 1; ++partial_coverage) {
114
115 struct interp_method * method_stack[] =
116 {yac_interp_method_avg_new(
117 YAC_INTERP_AVG_ARITHMETIC, partial_coverage),
118 yac_interp_method_fixed_new(-1.0), NULL};
119
120 struct yac_interp_weights * weights =
121 yac_interp_method_do_search(method_stack, interp_grid);
122
123 for (size_t i = 0; i < num_reorder_types; ++i) {
124
125 struct yac_interpolation * interpolation =
126 yac_interp_weights_get_interpolation(
127 weights, reorder_types[i], 1,
128 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
129
130 // check generated interpolation
131 {
132 double * src_field = NULL;
133 double ** src_fields = &src_field;
134 double * tgt_field = NULL;
135 double * ref_tgt_field = NULL;
136 double ref_global_tgt_field[9] =
137 {2.5, 3.5, 4.5, 6.5, 7.5, 8.5, -1, -1, -1};
138
139 if (is_tgt) {
140 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
141 ref_tgt_field = xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
142 for (size_t i = 0; i < grid_data.num_vertices; ++i)
143 ref_tgt_field[i] = ref_global_tgt_field[grid_data.vertex_ids[i]];
144 } else {
145 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
146 for (size_t i = 0; i < grid_data.num_vertices; ++i)
147 src_field[i] = (double)(grid_data.vertex_ids[i]);
148 }
149
150 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
151
152 if (is_tgt)
153 for (size_t i = 0; i < grid_data.num_vertices; ++i)
154 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-9)
155 PUT_ERR("wrong interpolation result");
156
157 free(src_field);
158 free(tgt_field);
159 free(ref_tgt_field);
160 }
161
162 yac_interpolation_delete(interpolation);
163 }
164
165 yac_interp_weights_delete(weights);
166 yac_interp_method_delete(method_stack);
167 } // partial_coverage
168 yac_interp_grid_delete(interp_grid);
169 yac_dist_grid_pair_delete(grid_pair);
170 yac_basic_grid_delete(grids[1]);
171 yac_basic_grid_delete(grids[0]);
172 }
173
174 {// Test 2a
175 // linear point interpolation (fixed is backup) with two source processes and a
176 // single target process the source fields have two struct points per grid
177 // (the test check how the interpolation methods handle NULL target coordinates )
178
179 // the global grid is a 4x2 grid:
180 // 08--14--09--15--10--16--11
181 // | | | |
182 // 08 03 10 04 12 05 13
183 // | | | |
184 // 04--07--05--09--06--11--07
185 // | | | |
186 // 01 00 03 01 05 02 06
187 // | | | |
188 // 00--00--01--02--02--04--03
189 //
190 // source mask:
191 // 0-------1-------1-------1
192 // | | | |
193 // | | | |
194 // | | | |
195 // 1-------1-------1-------1
196 // | | | |
197 // | | | |
198 // | | | |
199 // 1-------1-------1-------0
200 //
201 //---------------
202 // setup
203 //---------------
204
205 int is_tgt = split_comm_size == 1;
206 double coordinates_x[2][4] = {{0.0,1.0,2.0,3.0}, {0.5,1.5,2.5}};
207 double coordinates_y[2][3] = {{0.0,1.0,2.0}, {0.5,1.5,2.5}};
208 size_t const num_cells[2][2] = {{3,2}, {2,2}};
209 size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
210 size_t local_count[2][2][2] = {{{2,2},{2,2}}, {{2,2}}};
211 int global_src_mask[12] = {1, 1, 1, 0,
212 1, 1, 1, 1,
213 0, 1, 1, 1};
214 int with_halo = 0;
215 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
216 coordinates_x[is_tgt][i] *= YAC_RAD;
217 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
218 coordinates_y[is_tgt][i] *= YAC_RAD;
219
220 struct yac_basic_grid_data grid_data =
221 yac_generate_basic_grid_data_reg2d(
222 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
223 local_start[is_tgt][split_comm_rank],
224 local_count[is_tgt][split_comm_rank], with_halo);
225
226 struct yac_basic_grid * grids[2] =
227 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
228 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
229 int * src_mask = NULL;
230 if (!is_tgt) {
231 src_mask = xmalloc(grid_data.num_vertices * sizeof(*src_mask));
232 for (size_t i = 0; i < grid_data.num_vertices; ++i)
233 src_mask[i] = global_src_mask[grid_data.vertex_ids[i]];
234 yac_basic_grid_add_mask_nocpy(grids[0], YAC_LOC_CORNER, src_mask, NULL);
235 }
236
237 struct yac_dist_grid_pair * grid_pair =
238 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
239 struct yac_interp_field src_fields[] =
240 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
241 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
242 struct yac_interp_field tgt_field =
243 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
244
245 struct yac_interp_grid * interp_grid =
246 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
247 num_src_fields, src_fields, tgt_field);
248
249 struct interp_method * method_stack[] =
250 {yac_interp_method_avg_new(YAC_INTERP_AVG_ARITHMETIC, 0),
251 yac_interp_method_fixed_new(-1.0), NULL};
252
253 struct yac_interp_weights * weights =
254 yac_interp_method_do_search(method_stack, interp_grid);
255
256 for (size_t i = 0; i < num_reorder_types; ++i) {
257
258 struct yac_interpolation * interpolation =
259 yac_interp_weights_get_interpolation(
260 weights, reorder_types[i], 1,
261 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
262
263 // check generated interpolation
264 {
265 double * src_field = NULL;
266 double ** src_fields = &src_field;
267 double * tgt_field = NULL;
268 double * ref_tgt_field = NULL;
269 double ref_global_tgt_field[9] =
270 {2.5, 3.5, -1, -1, 7.5, 8.5, -1, -1, -1};
271
272 if (is_tgt) {
273 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
274 ref_tgt_field = xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
275 for (size_t i = 0; i < grid_data.num_vertices; ++i)
276 ref_tgt_field[i] = ref_global_tgt_field[grid_data.vertex_ids[i]];
277 } else {
278 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
279 for (size_t i = 0; i < grid_data.num_vertices; ++i)
280 src_field[i] = (double)(grid_data.vertex_ids[i]);
281 }
282
283 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
284
285 if (is_tgt)
286 for (size_t i = 0; i < grid_data.num_vertices; ++i)
287 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-9)
288 PUT_ERR("wrong interpolation result");
289
290 free(src_field);
291 free(tgt_field);
292 free(ref_tgt_field);
293 }
294
295 yac_interpolation_delete(interpolation);
296 }
297
298 yac_interp_weights_delete(weights);
299 yac_interp_method_delete(method_stack);
300 yac_interp_grid_delete(interp_grid);
301 yac_dist_grid_pair_delete(grid_pair);
302 yac_basic_grid_delete(grids[1]);
303 yac_basic_grid_delete(grids[0]);
304 }
305
306 {// Test 3
307 // linear distance weighted (no partial coverage) point interpolation
308 // with two source processes and a single target process
309
310 // the global grid is a 3x3 grid:
311 // 12--21--13--22--14--23--15
312 // | | | |
313 // 15 06 17 07 19 08 20
314 // | | | |
315 // 08--14--09--16--10--18--11
316 // | | | |
317 // 08 03 10 04 12 05 13
318 // | | | |
319 // 04--07--05--09--06--11--07
320 // | | | |
321 // 01 00 03 01 05 02 06
322 // | | | |
323 // 00--00--01--02--02--04--03
324 //
325 // source mask:
326 // 0-------1-------1-------0
327 // | | | |
328 // | | | |
329 // | | | |
330 // 1-------1-------1-------1
331 // | | | |
332 // | | | |
333 // | | | |
334 // 1-------1-------1-------1
335 // | | | |
336 // | | | |
337 // | | | |
338 // 0-------1-------1-------0
339 //
340 //---------------
341 // setup
342 //---------------
343
344 int is_tgt = split_comm_size == 1;
345 double coordinates_x[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
346 double coordinates_y[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
347 size_t const num_cells[2][2] = {{3,3}, {3,3}};
348 size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
349 size_t local_count[2][2][2] = {{{2,3},{2,3}}, {{3,3}}};
350 int global_src_mask[16] = {0, 1, 1, 0,
351 1, 1, 1, 1,
352 1, 1, 1, 1,
353 0, 1, 1, 0};
354 int with_halo = 0;
355 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
356 coordinates_x[is_tgt][i] *= YAC_RAD;
357 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
358 coordinates_y[is_tgt][i] *= YAC_RAD;
359
360 struct yac_basic_grid_data grid_data =
361 yac_generate_basic_grid_data_reg2d(
362 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
363 local_start[is_tgt][split_comm_rank],
364 local_count[is_tgt][split_comm_rank], with_halo);
365
366 struct yac_basic_grid * grids[2] =
367 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
368 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
369
370 int * src_mask = NULL;
371 if (!is_tgt) {
372 src_mask = xmalloc(grid_data.num_vertices * sizeof(*src_mask));
373 for (size_t i = 0; i < grid_data.num_vertices; ++i)
374 src_mask[i] = global_src_mask[grid_data.vertex_ids[i]];
375 yac_basic_grid_add_mask_nocpy(grids[0], YAC_LOC_CORNER, src_mask, NULL);
376 }
377
378 yac_coordinate_pointer tgt_cell_coordinates = NULL;
379 if (is_tgt) {
380 double cell_coordinates_x[3] = {0.75,1.75,2.75};
381 double cell_coordinates_y[3] = {0.75,1.75,2.75};
382 tgt_cell_coordinates = xmalloc(9 * sizeof(*tgt_cell_coordinates));
383 for (int i = 0, k = 0; i < 3; ++i)
384 for (int j = 0; j < 3; ++j, ++k)
385 LLtoXYZ_deg(
386 cell_coordinates_x[j], cell_coordinates_y[i],
387 tgt_cell_coordinates[k]);
388 yac_basic_grid_add_coordinates_nocpy(grids[0], YAC_LOC_CELL, tgt_cell_coordinates);
389 }
390
391 struct yac_dist_grid_pair * grid_pair =
392 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
393
394 struct yac_interp_field src_fields[] =
395 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
396 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
397 struct yac_interp_field tgt_field =
398 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX};
399
400 struct yac_interp_grid * interp_grid =
401 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
402 num_src_fields, src_fields, tgt_field);
403
404 struct interp_method * method_stack[] =
405 {yac_interp_method_avg_new(YAC_INTERP_AVG_DIST, 0),
406 yac_interp_method_fixed_new(-1.0), NULL};
407
408 struct yac_interp_weights * weights =
409 yac_interp_method_do_search(method_stack, interp_grid);
410
411 for (size_t i = 0; i < num_reorder_types; ++i) {
412
413 struct yac_interpolation * interpolation =
414 yac_interp_weights_get_interpolation(
415 weights, reorder_types[i], 1,
416 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
417
418 // check generated interpolation
419 {
420 double * src_field = NULL;
421 double ** src_fields = &src_field;
422 double * tgt_field = NULL;
423 double * ref_tgt_field = NULL;
424 double ref_global_tgt_field[9] =
425 {-1.0,//0.14963* 0.0 + 0.20075*( 1.0+ 4.0) + 0.44888* 5.0,
426 0.14963* 1.0 + 0.20075*( 2.0+ 5.0) + 0.44888* 6.0,
427 -1.0,//0.14963* 2.0 + 0.20075*( 3.0+ 6.0) + 0.44888* 7.0,
428 0.14963* 4.0 + 0.20075*( 5.0+ 8.0) + 0.44888* 9.0,
429 0.14963* 5.0 + 0.20075*( 6.0+ 9.0) + 0.44888*10.0,
430 0.14963* 6.0 + 0.20075*( 7.0+10.0) + 0.44888*11.0,
431 -1.0,//0.14963* 8.0 + 0.20075*( 9.0+12.0) + 0.44888*13.0,
432 0.14963* 9.0 + 0.20075*(10.0+13.0) + 0.44888*14.0,
433 -1.0};//0.14963*10.0 + 0.20075*(11.0+14.0) + 0.44888*15.0};
434
435 if (is_tgt) {
436 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
437 ref_tgt_field = xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
438 for (size_t i = 0; i < grid_data.num_vertices; ++i)
439 ref_tgt_field[i] = ref_global_tgt_field[grid_data.vertex_ids[i]];
440 } else {
441 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
442 for (size_t i = 0; i < grid_data.num_vertices; ++i)
443 src_field[i] = (double)(grid_data.vertex_ids[i]);
444 }
445
446 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
447
448 if (is_tgt)
449 for (size_t i = 0; i < grid_data.num_cells; ++i)
450 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-3)
451 PUT_ERR("wrong interpolation result");
452
453 free(src_field);
454 free(tgt_field);
455 free(ref_tgt_field);
456 }
457
458 yac_interpolation_delete(interpolation);
459 }
460
461 yac_interp_weights_delete(weights);
462 yac_interp_method_delete(method_stack);
463 yac_interp_grid_delete(interp_grid);
464 yac_dist_grid_pair_delete(grid_pair);
465 yac_basic_grid_delete(grids[1]);
466 yac_basic_grid_delete(grids[0]);
467 }
468
469 {// Test 4
470 // linear distance weighted (with partial coverage) point interpolation
471 // with two source processes and a single target process
472
473 // the global grid is a 3x3 grid:
474 // 12--21--13--22--14--23--15
475 // | | | |
476 // 15 06 17 07 19 08 20
477 // | | | |
478 // 08--14--09--16--10--18--11
479 // | | | |
480 // 08 03 10 04 12 05 13
481 // | | | |
482 // 04--07--05--09--06--11--07
483 // | | | |
484 // 01 00 03 01 05 02 06
485 // | | | |
486 // 00--00--01--02--02--04--03
487 //
488 // source mask:
489 // 0-------1-------1-------0
490 // | | | |
491 // | | | |
492 // | | | |
493 // 1-------0-------1-------1
494 // | | | |
495 // | | | |
496 // | | | |
497 // 1-------1-------0-------0
498 // | | | |
499 // | | | |
500 // | | | |
501 // 0-------1-------0-------0
502 //
503 //---------------
504 // setup
505 //---------------
506
507 int is_tgt = split_comm_size == 1;
508 double coordinates_x[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
509 double coordinates_y[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
510 size_t const num_cells[2][2] = {{3,3}, {3,3}};
511 size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
512 size_t local_count[2][2][2] = {{{2,3},{2,3}}, {{3,3}}};
513 int global_src_mask[16] = {0, 1, 0, 0,
514 1, 1, 0, 0,
515 1, 0, 1, 1,
516 0, 1, 1, 0};
517 int with_halo = 0;
518 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
519 coordinates_x[is_tgt][i] *= YAC_RAD;
520 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
521 coordinates_y[is_tgt][i] *= YAC_RAD;
522
523 struct yac_basic_grid_data grid_data =
524 yac_generate_basic_grid_data_reg2d(
525 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
526 local_start[is_tgt][split_comm_rank],
527 local_count[is_tgt][split_comm_rank], with_halo);
528
529 struct yac_basic_grid * grids[2] =
530 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
531 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
532
533 int * src_mask = NULL;
534 if (!is_tgt) {
535 src_mask = xmalloc(grid_data.num_vertices * sizeof(*src_mask));
536 for (size_t i = 0; i < grid_data.num_vertices; ++i)
537 src_mask[i] = global_src_mask[grid_data.vertex_ids[i]];
538 yac_basic_grid_add_mask_nocpy(grids[0], YAC_LOC_CORNER, src_mask, NULL);
539 }
540
541 yac_coordinate_pointer tgt_cell_coordinates = NULL;
542 if (is_tgt) {
543 double cell_coordinates_x[3] = {0.75,1.75,2.75};
544 double cell_coordinates_y[3] = {0.75,1.75,2.75};
545 tgt_cell_coordinates = xmalloc(9 * sizeof(*tgt_cell_coordinates));
546 for (int i = 0, k = 0; i < 3; ++i)
547 for (int j = 0; j < 3; ++j, ++k)
548 LLtoXYZ_deg(
549 cell_coordinates_x[j], cell_coordinates_y[i],
550 tgt_cell_coordinates[k]);
551 yac_basic_grid_add_coordinates_nocpy(grids[0], YAC_LOC_CELL, tgt_cell_coordinates);
552 }
553
554 struct yac_dist_grid_pair * grid_pair =
555 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
556
557 struct yac_interp_field src_fields[] =
558 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
559 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
560 struct yac_interp_field tgt_field =
561 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX};
562
563 struct yac_interp_grid * interp_grid =
564 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
565 num_src_fields, src_fields, tgt_field);
566
567 struct interp_method * method_stack[] =
568 {yac_interp_method_avg_new(YAC_INTERP_AVG_DIST, 1),
569 yac_interp_method_fixed_new(-1.0), NULL};
570
571 struct yac_interp_weights * weights =
572 yac_interp_method_do_search(method_stack, interp_grid);
573
574 for (size_t i = 0; i < num_reorder_types; ++i) {
575
576 struct yac_interpolation * interpolation =
577 yac_interp_weights_get_interpolation(
578 weights, reorder_types[i], 1,
579 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
580
581 // check generated interpolation
582 {
583 double * src_field = NULL;
584 double ** src_fields = &src_field;
585 double * tgt_field = NULL;
586 double * ref_tgt_field = NULL;
587 double inv_dist_a = sqrt(8.0);
588 double inv_dist_b = sqrt(8.0/5.0);
589 double inv_dist_c = sqrt(8.0/9.0);
590 double ref_global_tgt_field[9] =
591 {(inv_dist_b*( 1.0+ 4.0) + inv_dist_a* 5.0)/(inv_dist_a + 2.0 * inv_dist_b),
592 (inv_dist_c* 1.0 + inv_dist_b* 5.0)/(inv_dist_b + inv_dist_c),
593 -1.0,
594 (inv_dist_c* 4.0 + inv_dist_b*( 5.0+ 8.0))/(2.0 * inv_dist_b + inv_dist_c),
595 (inv_dist_c* 5.0 + inv_dist_a*10.0)/(inv_dist_a + inv_dist_c),
596 (inv_dist_b*10.0 + inv_dist_a*11.0)/(inv_dist_a + inv_dist_b),
597 (inv_dist_c* 8.0 + inv_dist_a*13.0)/(inv_dist_a + inv_dist_c),
598 (inv_dist_b*(10.0+13.0) + inv_dist_a*14.0)/(inv_dist_a + 2.0 * inv_dist_b),
599 (inv_dist_c*10.0 + inv_dist_b*(11.0+14.0))/(2.0 * inv_dist_b + inv_dist_c)};
600
601 if (is_tgt) {
602 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
603 ref_tgt_field = xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
604 for (size_t i = 0; i < grid_data.num_vertices; ++i)
605 ref_tgt_field[i] = ref_global_tgt_field[grid_data.vertex_ids[i]];
606 } else {
607 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
608 for (size_t i = 0; i < grid_data.num_vertices; ++i)
609 src_field[i] = (double)(grid_data.vertex_ids[i]);
610 }
611
612 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
613
614 if (is_tgt)
615 for (size_t i = 0; i < grid_data.num_cells; ++i)
616 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-3)
617 PUT_ERR("wrong interpolation result");
618
619 free(src_field);
620 free(tgt_field);
621 free(ref_tgt_field);
622 }
623
624 yac_interpolation_delete(interpolation);
625 }
626
627 yac_interp_weights_delete(weights);
628 yac_interp_method_delete(method_stack);
629 yac_interp_grid_delete(interp_grid);
630 yac_dist_grid_pair_delete(grid_pair);
631 yac_basic_grid_delete(grids[1]);
632 yac_basic_grid_delete(grids[0]);
633 }
634
635 {// Test 5
636 // linear distance weighted (with partial coverage) point interpolation
637 // with two source processes and a single target process
638 // (target points are on vertices of the source grid
639
640 // the global grid is a 3x3 grid:
641 // 12--21--13--22--14--23--15
642 // | | | |
643 // 15 06 17 07 19 08 20
644 // | | | |
645 // 08--14--09--16--10--18--11
646 // | | | |
647 // 08 03 10 04 12 05 13
648 // | | | |
649 // 04--07--05--09--06--11--07
650 // | | | |
651 // 01 00 03 01 05 02 06
652 // | | | |
653 // 00--00--01--02--02--04--03
654 //
655 //---------------
656 // setup
657 //---------------
658
659 int is_tgt = split_comm_size == 1;
660 double coordinates_x[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
661 double coordinates_y[2][4] = {{0.0,1.0,2.0,3.0}, {0.0,1.0,2.0,3.0}};
662 size_t const num_cells[2][2] = {{3,3}, {3,3}};
663 size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
664 size_t local_count[2][2][2] = {{{2,3},{2,3}}, {{3,3}}};
665 int with_halo = 0;
666 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
667 coordinates_x[is_tgt][i] *= YAC_RAD;
668 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
669 coordinates_y[is_tgt][i] *= YAC_RAD;
670
671 struct yac_basic_grid_data grid_data =
672 yac_generate_basic_grid_data_reg2d(
673 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
674 local_start[is_tgt][split_comm_rank],
675 local_count[is_tgt][split_comm_rank], with_halo);
676
677 struct yac_basic_grid * grids[2] =
678 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
679 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
680
681 yac_coordinate_pointer tgt_cell_coordinates = NULL;
682 if (is_tgt) {
683 double cell_coordinates_x[3] = {1.0,2.0,3.0};
684 double cell_coordinates_y[3] = {1.0,2.0,3.0};
685 tgt_cell_coordinates = xmalloc(9 * sizeof(*tgt_cell_coordinates));
686 for (int i = 0, k = 0; i < 3; ++i)
687 for (int j = 0; j < 3; ++j, ++k)
688 LLtoXYZ_deg(
689 cell_coordinates_x[j], cell_coordinates_y[i],
690 tgt_cell_coordinates[k]);
691 yac_basic_grid_add_coordinates_nocpy(grids[0], YAC_LOC_CELL, tgt_cell_coordinates);
692 }
693
694 struct yac_dist_grid_pair * grid_pair =
695 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
696
697 struct yac_interp_field src_fields[] =
698 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
699 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
700 struct yac_interp_field tgt_field =
701 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX};
702
703 struct yac_interp_grid * interp_grid =
704 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
705 num_src_fields, src_fields, tgt_field);
706
707 struct interp_method * method_stack[] =
708 {yac_interp_method_avg_new(YAC_INTERP_AVG_DIST, 1),
709 yac_interp_method_fixed_new(-1.0), NULL};
710
711 struct yac_interp_weights * weights =
712 yac_interp_method_do_search(method_stack, interp_grid);
713
714 for (size_t i = 0; i < num_reorder_types; ++i) {
715
716 struct yac_interpolation * interpolation =
717 yac_interp_weights_get_interpolation(
718 weights, reorder_types[i], 1,
719 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
720
721 // check generated interpolation
722 {
723 double * src_field = NULL;
724 double ** src_fields = &src_field;
725 double * tgt_field = NULL;
726 double * ref_tgt_field = NULL;
727 double ref_global_tgt_field[9] =
728 {5.0, 6.0, 7.0, 9.0, 10.0, 11.0, 13.0, 14.0, 15.0};
729
730 if (is_tgt) {
731 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
732 ref_tgt_field = xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
733 for (size_t i = 0; i < grid_data.num_vertices; ++i)
734 ref_tgt_field[i] = ref_global_tgt_field[grid_data.vertex_ids[i]];
735 } else {
736 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
737 for (size_t i = 0; i < grid_data.num_vertices; ++i)
738 src_field[i] = (double)(grid_data.vertex_ids[i]);
739 }
740
741 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
742
743 if (is_tgt)
744 for (size_t i = 0; i < grid_data.num_cells; ++i)
745 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-9)
746 PUT_ERR("wrong interpolation result");
747
748 free(src_field);
749 free(tgt_field);
750 free(ref_tgt_field);
751 }
752
753 yac_interpolation_delete(interpolation);
754 }
755
756 yac_interp_weights_delete(weights);
757 yac_interp_method_delete(method_stack);
758 yac_interp_grid_delete(interp_grid);
759 yac_dist_grid_pair_delete(grid_pair);
760 yac_basic_grid_delete(grids[1]);
761 yac_basic_grid_delete(grids[0]);
762 }
763
764 {// Test 6
765 // tests cell based source field
766
767 // the global grid is a 7x7 grid:
768 // 56-----57-----58-----59-----60-----61-----62-----63
769 // | | | | | | | |
770 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
771 // | | | | | | | |
772 // 48-----49-----50-----51-----52-----53-----54-----55
773 // | | | | | | | |
774 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
775 // | | | | | | | |
776 // 40-----41-----42-----43-----44-----45-----46-----47
777 // | | | | | | | |
778 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
779 // | | | | | | | |
780 // 32-----33-----34-----35-----36-----37-----38-----39
781 // | | | | | | | |
782 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
783 // | | | | | | | |
784 // 24-----25-----26-----27-----28-----29-----30-----31
785 // | | | | | | | |
786 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
787 // | | | | | | | |
788 // 16-----17-----18-----19-----20-----21-----22-----23
789 // | | | | | | | |
790 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
791 // | | | | | | | |
792 // 08-----09-----10-----11-----12-----13-----14-----15
793 // | | | | | | | |
794 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
795 // | | | | | | | |
796 // 00-----01-----02-----03-----04-----05-----06-----07
797 //
798 //---------------
799 // setup
800 //---------------
801
802 int is_tgt = split_comm_size == 1;
803 double coordinates_x[8] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
804 double coordinates_y[8] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
805 size_t const num_cells[2] = {7,7};
806 size_t local_start[2][2][2] = {{{0,0},{0,4}}, {{0,0}}};
807 size_t local_count[2][2][2] = {{{7,4},{7,3}}, {{7,7}}};
808 int with_halo = 0;
809 for (size_t i = 0; i <= num_cells[0]; ++i)
810 coordinates_x[i] *= YAC_RAD;
811 for (size_t i = 0; i <= num_cells[1]; ++i)
812 coordinates_y[i] *= YAC_RAD;
813
814 struct yac_basic_grid_data grid_data =
815 yac_generate_basic_grid_data_reg2d(
816 coordinates_x, coordinates_y, num_cells,
817 local_start[is_tgt][split_comm_rank],
818 local_count[is_tgt][split_comm_rank], with_halo);
819
820 struct yac_basic_grid * grids[2] =
821 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
822 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
823
824 yac_coordinate_pointer src_point_coordinates = NULL;
825 if (!is_tgt) {
826 src_point_coordinates =
827 xmalloc(grid_data.num_cells * sizeof(*src_point_coordinates));
828 for (size_t i = 0; i < grid_data.num_cells; ++i) {
829 double * middle_point = src_point_coordinates[i];
830 for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
831 size_t * curr_vertices =
832 grid_data.cell_to_vertex + grid_data.cell_to_vertex_offsets[i];
833 size_t curr_num_vertices = grid_data.num_vertices_per_cell[i];
834 for (size_t j = 0; j < curr_num_vertices; ++j) {
835 double * curr_vertex_coord =
836 grid_data.vertex_coordinates[curr_vertices[j]];
837 for (size_t k = 0; k < 3; ++k)
838 middle_point[k] += curr_vertex_coord[k];
839 }
840 normalise_vector(middle_point);
841 }
842 yac_basic_grid_add_coordinates_nocpy(grids[0], YAC_LOC_CELL, src_point_coordinates);
843 }
844
845 struct yac_dist_grid_pair * grid_pair =
846 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
847 struct yac_interp_field src_fields[] =
848 {{.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX}};
849 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
850 struct yac_interp_field tgt_field =
851 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
852
853 struct yac_interp_grid * interp_grid =
854 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
855 num_src_fields, src_fields, tgt_field);
856
857 struct interp_method * method_stack[] =
858 {yac_interp_method_avg_new(YAC_INTERP_AVG_ARITHMETIC, 0),
859 yac_interp_method_fixed_new(-1.0), NULL};
860
861 struct yac_interp_weights * weights =
862 yac_interp_method_do_search(method_stack, interp_grid);
863
864 for (size_t i = 0; i < num_reorder_types; ++i) {
865
866 struct yac_interpolation * interpolation =
867 yac_interp_weights_get_interpolation(
868 weights, reorder_types[i], 1,
869 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
870
871 // check generated interpolation
872 {
873 double * src_field = NULL;
874 double ** src_fields = &src_field;
875 double * tgt_field = NULL;
876 double * ref_tgt_field = NULL;
877 double ref_global_tgt_field[64] =
878 {-4,-4,-4,-4,-4,-4,-4,-4,
879 -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,
880 -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,
881 -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,
882 -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,
883 -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,
884 -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,
885 -4,-4,-4,-4,-4,-4,-4,-4};
886 for (size_t i = 0;
887 i < sizeof(ref_global_tgt_field) / sizeof(ref_global_tgt_field[0]);
888 ++i)
889 ref_global_tgt_field[i] /= 4.0;
890
891 if (is_tgt) {
892 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
893 ref_tgt_field = xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
894 for (size_t i = 0; i < grid_data.num_vertices; ++i)
895 ref_tgt_field[i] = ref_global_tgt_field[grid_data.vertex_ids[i]];
896 } else {
897 src_field = xmalloc(grid_data.num_cells * sizeof(*src_field));
898 for (size_t i = 0; i < grid_data.num_cells; ++i)
899 src_field[i] = (double)(grid_data.cell_ids[i]);
900 }
901
902 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
903
904 if (is_tgt)
905 for (size_t i = 0; i < grid_data.num_vertices; ++i)
906 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-9)
907 PUT_ERR("wrong interpolation result");
908
909 free(src_field);
910 free(tgt_field);
911 free(ref_tgt_field);
912 }
913
914 yac_interpolation_delete(interpolation);
915 }
916
917 yac_interp_weights_delete(weights);
918 yac_interp_method_delete(method_stack);
919 yac_interp_grid_delete(interp_grid);
920 yac_dist_grid_pair_delete(grid_pair);
921 yac_basic_grid_delete(grids[1]);
922 yac_basic_grid_delete(grids[0]);
923 }
924
925 {// Test 7
926 // barycentric coordinates weighted (with partial coverage) point
927 // interpolation with two source processes and a single target process
928
929 // the global source grid is a 1x1 grid:
930 // 02--03--03
931 // | |
932 // 01 00 02
933 // | |
934 // 00--00--01
935 //
936 //---------------
937 // setup
938 //---------------
939
940 int is_tgt = split_comm_size == 1;
941 double src_coordinates_x[2][2] = {{-0.5,0.5}, {0.5,-0.5}};
942 double src_coordinates_y[2][2][2] =
943 {{{-0.5,0.5}, {0.5,-0.5}}, {{89,90}, {90,89}}};
944 double src_coordinates_x_unstruct[4];
945 double src_coordinates_y_unstruct[4];
946 size_t const src_num_cells[2] = {1,1};
947 size_t src_local_start[2] = {0,0};
948 size_t src_local_count[2] = {1,1};
949 int src_cell_mask[4] = {1,1, 0,1};
950 yac_int src_reorder[2][2][4] =
951 {{{0,1,2,3},{2,3,0,1}},{{1,0,3,2},{3,2,1,0}}};
952 double tgt_coordinates_x[5] = {-0.5,-0.25,0.0,0.25,0.5};
953 double tgt_coordinates_y[2][5] =
954 {{-0.5,-0.25,0.0,0.25,0.5}, {89.00,89.25,89.50,89.75,90.00}};
955 size_t const tgt_num_cells[2] = {4,4};
956 size_t tgt_local_start[2] = {0,0};
957 size_t tgt_local_count[2] = {4,4};
958
959 for (size_t i = 0; i < 2 * 2; ++i)
960 (&(src_coordinates_x[0][0]))[i] *= YAC_RAD;
961 for (size_t i = 0; i < 2 * 2 * 2; ++i)
962 (&(src_coordinates_y[0][0][0]))[i] *= YAC_RAD;
963 for (size_t i = 0; i < 5; ++i) tgt_coordinates_x[i] *= YAC_RAD;
964 for (size_t i = 0; i < 2 * 5; ++i)
965 (&(tgt_coordinates_y[0][0]))[i] *= YAC_RAD;
966
967 for (int partial_coverage = 0; partial_coverage < 2; ++partial_coverage) {
968 for (int is_reg = 0; is_reg < 2; ++is_reg) {
969 for (int at_pole = 0; at_pole <= is_reg; ++at_pole) {
970 for (int with_mask = 0; with_mask < 2; ++with_mask) {
971 for (int order_x = 0; order_x < 2; ++order_x) {
972 for (int order_y = 0; order_y < 2; ++order_y) {
973
974 struct yac_basic_grid_data grid_data;
975
976 if (is_tgt) {
977 grid_data =
978 yac_generate_basic_grid_data_reg2d(
979 tgt_coordinates_x, tgt_coordinates_y[at_pole],
980 tgt_num_cells, tgt_local_start, tgt_local_count, 0);
981 } else {
982 if (is_reg) {
983 grid_data =
984 yac_generate_basic_grid_data_reg2d(
985 src_coordinates_x[order_x],
986 src_coordinates_y[at_pole][order_y], src_num_cells,
987 src_local_start, src_local_count, 0);
988 } else {
989 for (int i = 0; i < 4; ++i) {
990 src_coordinates_x_unstruct[i] =
991 src_coordinates_x[order_x][i&1];
992 src_coordinates_y_unstruct[i] =
993 src_coordinates_y[at_pole][order_y][i>>1];
994 }
995 grid_data =
996 yac_generate_basic_grid_data_unstruct(
997 4, 1, (int[]){4},
998 src_coordinates_x_unstruct,
999 src_coordinates_y_unstruct, (int[]){0,1,3,2});
1000 }
1001 }
1002
1003 if (!is_tgt) {
1004
1005 yac_int * vertex_ids =
1006 xmalloc(grid_data.num_vertices * sizeof(*vertex_ids));
1007 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1008 vertex_ids[i] = (yac_int)(src_reorder[order_x][order_y][i]);
1009 grid_data.vertex_ids = vertex_ids;
1010 }
1011
1012 struct yac_basic_grid * grids[2] =
1013 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
1014 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
1015
1016 if (with_mask && !is_tgt) {
1017 int * mask = xmalloc(grid_data.num_vertices * sizeof(*mask));
1018 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1019 mask[i] =
1020 src_cell_mask[src_reorder[order_x][order_y][i]];
1021 yac_basic_grid_add_mask_nocpy(grids[0], YAC_LOC_CORNER, mask, NULL);
1022 }
1023
1024 struct yac_dist_grid_pair * grid_pair =
1025 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
1026
1027 struct yac_interp_field src_fields[] =
1028 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX,
1029 .masks_idx = (with_mask)?0:SIZE_MAX}};
1030 size_t num_src_fields =
1031 sizeof(src_fields) / sizeof(src_fields[0]);
1032 struct yac_interp_field tgt_field =
1033 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX,
1034 .masks_idx = SIZE_MAX};
1035
1036 struct yac_interp_grid * interp_grid =
1037 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
1038 num_src_fields, src_fields, tgt_field);
1039
1040 struct interp_method * method_stack[] =
1041 {yac_interp_method_avg_new(
1042 YAC_INTERP_AVG_BARY, partial_coverage),
1043 yac_interp_method_fixed_new(-1.0), NULL};
1044
1045 struct yac_interp_weights * weights =
1046 yac_interp_method_do_search(method_stack, interp_grid);
1047
1048 for (size_t i = 0; i < num_reorder_types; ++i) {
1049
1050 struct yac_interpolation * interpolation =
1051 yac_interp_weights_get_interpolation(
1052 weights, reorder_types[i], 1,
1053 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1054
1055 // check generated interpolation
1056 {
1057 double * src_field = NULL;
1058 double ** src_fields = &src_field;
1059 double * tgt_field = NULL;
1060 double * ref_tgt_field = NULL;
1061 double ref_global_tgt_field[25];
1062 double ref_weights[2][2][2][25][4] =
1063 {{{{{1.00, 0.00, 0.00, 0.00},
1064 {0.75, 0.25, 0.00, 0.00},
1065 {0.50, 0.50, 0.00, 0.00},
1066 {0.25, 0.75, 0.00, 0.00},
1067 {0.00, 1.00, 0.00, 0.00},
1068 {0.75, 0.00, 0.25, 0.00},
1069 {0.75, 0.00, 0.00, 0.25},
1070 {0.50, 0.25, 0.00, 0.25},
1071 {0.25, 0.50, 0.00, 0.25},
1072 {0.00, 0.75, 0.00, 0.25},
1073 {0.50, 0.00, 0.50, 0.00},
1074 {0.50, 0.00, 0.25, 0.25},
1075 {0.50, 0.00, 0.00, 0.50},
1076 {0.25, 0.25, 0.00, 0.50},
1077 {0.00, 0.50, 0.00, 0.50},
1078 {0.25, 0.00, 0.75, 0.00},
1079 {0.25, 0.00, 0.50, 0.25},
1080 {0.25, 0.00, 0.25, 0.50},
1081 {0.25, 0.00, 0.00, 0.75},
1082 {0.00, 0.25, 0.00, 0.75},
1083 {0.00, 0.00, 1.00, 0.00},
1084 {0.00, 0.00, 0.75, 0.25},
1085 {0.00, 0.00, 0.50, 0.50},
1086 {0.00, 0.00, 0.25, 0.75},
1087 {0.00, 0.00, 0.00, 1.00}},
1088 {{1.00, 0.00, 0.00, 0.00},
1089 {0.75, 0.25, 0.00, 0.00},
1090 {0.50, 0.50, 0.00, 0.00},
1091 {0.25, 0.75, 0.00, 0.00},
1092 {0.00, 1.00, 0.00, 0.00},
1093 {0.00, 0.00, 0.00, 0.00},//{0.75, 0.00, 0.25, 0.00}
1094 {0.75, 0.00, 0.00, 0.25},
1095 {0.50, 0.25, 0.00, 0.25},
1096 {0.25, 0.50, 0.00, 0.25},
1097 {0.00, 0.75, 0.00, 0.25},
1098 {0.00, 0.00, 0.00, 0.00},//{0.50, 0.00, 0.50, 0.00}
1099 {0.00, 0.00, 0.00, 0.00},//{0.50, 0.00, 0.25, 0.25}
1100 {0.50, 0.00, 0.00, 0.50},
1101 {0.25, 0.25, 0.00, 0.50},
1102 {0.00, 0.50, 0.00, 0.50},
1103 {0.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.75, 0.00}
1104 {0.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.50, 0.25}
1105 {0.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.25, 0.50}
1106 {0.25, 0.00, 0.00, 0.75},
1107 {0.00, 0.25, 0.00, 0.75},
1108 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 1.00, 0.00}
1109 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 0.75, 0.25}
1110 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 0.50, 0.50}
1111 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 0.25, 0.75}
1112 {0.00, 0.00, 0.00, 1.00}}},
1113 {{{1.00, 0.00, 0.00, 0.00},
1114 {0.75, 0.25, 0.00, 0.00},
1115 {0.50, 0.50, 0.00, 0.00},
1116 {0.25, 0.75, 0.00, 0.00},
1117 {0.00, 1.00, 0.00, 0.00},
1118 {0.75, 0.00, 0.25, 0.00},
1119 {0.75, 0.00, 0.00, 0.25},
1120 {0.50, 0.25, 0.00, 0.25},
1121 {0.25, 0.50, 0.00, 0.25},
1122 {0.00, 0.75, 0.00, 0.25},
1123 {0.50, 0.00, 0.50, 0.00},
1124 {0.50, 0.00, 0.25, 0.25},
1125 {0.50, 0.00, 0.00, 0.50},
1126 {0.25, 0.25, 0.00, 0.50},
1127 {0.00, 0.50, 0.00, 0.50},
1128 {0.25, 0.00, 0.75, 0.00},
1129 {0.25, 0.00, 0.50, 0.25},
1130 {0.25, 0.00, 0.25, 0.50},
1131 {0.25, 0.00, 0.00, 0.75},
1132 {0.00, 0.25, 0.00, 0.75},
1133 {0.00, 0.00, 1.00, 0.00},
1134 {0.00, 0.00, 0.75, 0.25},
1135 {0.00, 0.00, 0.50, 0.50},
1136 {0.00, 0.00, 0.25, 0.75},
1137 {0.00, 0.00, 0.00, 1.00}},
1138 {{1.00, 0.00, 0.00, 0.00},
1139 {0.75, 0.25, 0.00, 0.00},
1140 {0.50, 0.50, 0.00, 0.00},
1141 {0.25, 0.75, 0.00, 0.00},
1142 {0.00, 1.00, 0.00, 0.00},
1143 {1.00, 0.00, 0.00, 0.00},//{0.75, 0.00, 0.25, 0.00}
1144 {0.75, 0.00, 0.00, 0.25},
1145 {0.50, 0.25, 0.00, 0.25},
1146 {0.25, 0.50, 0.00, 0.25},
1147 {0.00, 0.75, 0.00, 0.25},
1148 {1.00, 0.00, 0.00, 0.00},//{0.50, 0.00, 0.50, 0.00}
1149 {0.50/0.75, 0.00, 0.00, 0.25/0.75},//{0.50, 0.00, 0.25, 0.25}
1150 {0.50, 0.00, 0.00, 0.50},
1151 {0.25, 0.25, 0.00, 0.50},
1152 {0.00, 0.50, 0.00, 0.50},
1153 {1.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.75, 0.00}
1154 {0.50, 0.00, 0.00, 0.50},//{0.25, 0.00, 0.50, 0.25}
1155 {0.25/0.75, 0.00, 0.00, 0.50/0.75},//{0.25, 0.00, 0.25, 0.50}
1156 {0.25, 0.00, 0.00, 0.75},
1157 {0.00, 0.25, 0.00, 0.75},
1158 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 1.00, 0.00}
1159 {0.00, 0.00, 0.00, 1.00},//{0.00, 0.00, 0.75, 0.25}
1160 {0.00, 0.00, 0.00, 1.00},//{0.00, 0.00, 0.50, 0.50}
1161 {0.00, 0.00, 0.00, 1.00},//{0.00, 0.00, 0.25, 0.75}
1162 {0.00, 0.00, 0.00, 1.00}}}},
1163 {{{{1.00, 0.00, 0.00, 0.00},
1164 {0.75, 0.25, 0.00, 0.00},
1165 {0.50, 0.50, 0.00, 0.00},
1166 {0.25, 0.75, 0.00, 0.00},
1167 {0.00, 1.00, 0.00, 0.00},
1168 {0.75, 0.00, 0.25, 0.00},
1169 {0.75, 0.00, 0.00, 0.25},
1170 {0.50, 0.25, 0.00, 0.25},
1171 {0.25, 0.50, 0.00, 0.25},
1172 {0.00, 0.75, 0.00, 0.25},
1173 {0.50, 0.00, 0.50, 0.00},
1174 {0.50, 0.00, 0.25, 0.25},
1175 {0.50, 0.00, 0.00, 0.50},
1176 {0.25, 0.25, 0.00, 0.50},
1177 {0.00, 0.50, 0.00, 0.50},
1178 {0.25, 0.00, 0.75, 0.00},
1179 {0.25, 0.00, 0.50, 0.25},
1180 {0.25, 0.00, 0.25, 0.50},
1181 {0.25, 0.00, 0.00, 0.75},
1182 {0.00, 0.25, 0.00, 0.75},
1183 {0.00, 0.00, 0.50, 0.50},
1184 {0.00, 0.00, 0.50, 0.50},
1185 {0.00, 0.00, 0.50, 0.50},
1186 {0.00, 0.00, 0.50, 0.50},
1187 {0.00, 0.00, 0.50, 0.50}},
1188 {{1.00, 0.00, 0.00, 0.00},
1189 {0.75, 0.25, 0.00, 0.00},
1190 {0.50, 0.50, 0.00, 0.00},
1191 {0.25, 0.75, 0.00, 0.00},
1192 {0.00, 1.00, 0.00, 0.00},
1193 {0.00, 0.00, 0.00, 0.00},//{0.75, 0.00, 0.25, 0.00}
1194 {0.75, 0.00, 0.00, 0.25},
1195 {0.50, 0.25, 0.00, 0.25},
1196 {0.25, 0.50, 0.00, 0.25},
1197 {0.00, 0.75, 0.00, 0.25},
1198 {0.00, 0.00, 0.00, 0.00},//{0.50, 0.00, 0.50, 0.00}
1199 {0.00, 0.00, 0.00, 0.00},//{0.50, 0.00, 0.25, 0.25}
1200 {0.50, 0.00, 0.00, 0.50},
1201 {0.25, 0.25, 0.00, 0.50},
1202 {0.00, 0.50, 0.00, 0.50},
1203 {0.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.75, 0.00}
1204 {0.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.50, 0.25}
1205 {0.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.25, 0.50}
1206 {0.25, 0.00, 0.00, 0.75},
1207 {0.00, 0.25, 0.00, 0.75},
1208 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 1.00, 0.00}
1209 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 0.75, 0.25}
1210 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 0.50, 0.50}
1211 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 0.25, 0.75}
1212 {0.00, 0.00, 0.00, 0.00}}},//{0.00, 0.00, 0.00, 1.00}
1213 {{{1.00, 0.00, 0.00, 0.00},
1214 {0.75, 0.25, 0.00, 0.00},
1215 {0.50, 0.50, 0.00, 0.00},
1216 {0.25, 0.75, 0.00, 0.00},
1217 {0.00, 1.00, 0.00, 0.00},
1218 {0.75, 0.00, 0.25, 0.00},
1219 {0.75, 0.00, 0.00, 0.25},
1220 {0.50, 0.25, 0.00, 0.25},
1221 {0.25, 0.50, 0.00, 0.25},
1222 {0.00, 0.75, 0.00, 0.25},
1223 {0.50, 0.00, 0.50, 0.00},
1224 {0.50, 0.00, 0.25, 0.25},
1225 {0.50, 0.00, 0.00, 0.50},
1226 {0.25, 0.25, 0.00, 0.50},
1227 {0.00, 0.50, 0.00, 0.50},
1228 {0.25, 0.00, 0.75, 0.00},
1229 {0.25, 0.00, 0.50, 0.25},
1230 {0.25, 0.00, 0.25, 0.50},
1231 {0.25, 0.00, 0.00, 0.75},
1232 {0.00, 0.25, 0.00, 0.75},
1233 {0.00, 0.00, 0.50, 0.50},
1234 {0.00, 0.00, 0.50, 0.50},
1235 {0.00, 0.00, 0.50, 0.50},
1236 {0.00, 0.00, 0.50, 0.50},
1237 {0.00, 0.00, 0.50, 0.50}},
1238 {{1.00, 0.00, 0.00, 0.00},
1239 {0.75, 0.25, 0.00, 0.00},
1240 {0.50, 0.50, 0.00, 0.00},
1241 {0.25, 0.75, 0.00, 0.00},
1242 {0.00, 1.00, 0.00, 0.00},
1243 {1.00, 0.00, 0.00, 0.00},//{0.75, 0.00, 0.25, 0.00}
1244 {0.75, 0.00, 0.00, 0.25},
1245 {0.50, 0.25, 0.00, 0.25},
1246 {0.25, 0.50, 0.00, 0.25},
1247 {0.00, 0.75, 0.00, 0.25},
1248 {1.00, 0.00, 0.00, 0.00},//{0.50, 0.00, 0.50, 0.00}
1249 {0.50/0.75, 0.00, 0.00, 0.25/0.75},//{0.50, 0.00, 0.25, 0.25}
1250 {0.50, 0.00, 0.00, 0.50},
1251 {0.25, 0.25, 0.00, 0.50},
1252 {0.00, 0.50, 0.00, 0.50},
1253 {1.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.75, 0.00}
1254 {0.50, 0.00, 0.00, 0.50},//{0.25, 0.00, 0.50, 0.25}
1255 {0.25/0.75, 0.00, 0.00, 0.50/0.75},//{0.25, 0.00, 0.25, 0.50}
1256 {0.25, 0.00, 0.00, 0.75},
1257 {0.00, 0.25, 0.00, 0.75},
1258 {0.00, 0.00, 0.00, 1.00},//{0.00, 0.00, 1.00, 0.00}
1259 {0.00, 0.00, 0.00, 1.00},//{0.00, 0.00, 0.75, 0.25}
1260 {0.00, 0.00, 0.00, 1.00},//{0.00, 0.00, 0.50, 0.50}
1261 {0.00, 0.00, 0.00, 1.00},//{0.00, 0.00, 0.25, 0.75}
1262 {0.00, 0.00, 0.00, 1.00}}}}};
1263
1264 for (size_t j = 0; j < 25; ++j) {
1265 ref_global_tgt_field[j] =
1266 (double)(src_reorder[order_x][order_y][0] + 1) *
1267 ref_weights[at_pole][partial_coverage][with_mask][j][0] +
1268 (double)(src_reorder[order_x][order_y][1] + 1) *
1269 ref_weights[at_pole][partial_coverage][with_mask][j][1] +
1270 (double)(src_reorder[order_x][order_y][2] + 1) *
1271 ref_weights[at_pole][partial_coverage][with_mask][j][2] +
1272 (double)(src_reorder[order_x][order_y][3] + 1) *
1273 ref_weights[at_pole][partial_coverage][with_mask][j][3];
1274 if (ref_global_tgt_field[j] == 0.0)
1275 ref_global_tgt_field[j] = -1.0;
1276 }
1277
1278 if (is_tgt) {
1279 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
1280 ref_tgt_field = xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
1281 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1282 ref_tgt_field[i] = ref_global_tgt_field[i];
1283 } else {
1284 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
1285 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1286 src_field[i] = (double)(i) + 1.0;
1287 }
1288
1289 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
1290
1291 if (is_tgt)
1292 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1293 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-3)
1294 PUT_ERR("wrong interpolation result");
1295
1296 free(src_field);
1297 free(tgt_field);
1298 free(ref_tgt_field);
1299 }
1300
1301 yac_interpolation_delete(interpolation);
1302 }
1303
1304 yac_interp_weights_delete(weights);
1305 yac_interp_method_delete(method_stack);
1306 yac_interp_grid_delete(interp_grid);
1307 yac_dist_grid_pair_delete(grid_pair);
1308 yac_basic_grid_delete(grids[1]);
1309 yac_basic_grid_delete(grids[0]);
1310 }
1311 }
1312 }
1313 }
1314 }
1315 }
1316 }
1317
1318 {// Test 8
1319 // barycentric coordinates weighted (with partial coverage) point
1320 // interpolation with two source processes and a single target process
1321
1322 // the global source grid:
1323 // 02------03
1324 // | / |
1325 // | / |
1326 // | / |
1327 // 00------01
1328 //
1329 //---------------
1330 // setup
1331 //---------------
1332
1333 int is_tgt = split_comm_size == 1;
1334 size_t src_num_vertices = 4;
1335 size_t src_num_cells = 2;
1336 int src_num_vertices_per_cell[] = {3, 3};
1337 int src_cell_to_vertex[] = {0,1,3, 0,2,3};
1338 double src_coordinates_x[4] = {-0.5,0.5,-0.5,0.5};
1339 double src_coordinates_y[4] = {-0.5,-0.5,0.5,0.5};
1340 size_t tgt_num_vertices[2] = {5,5};
1341 int tgt_cyclic[2] = {0,0};
1342 double tgt_coordinates_x[5] = {-0.5,-0.25,0.0,0.25,0.5};
1343 double tgt_coordinates_y[5] = {-0.5,-0.25,0.0,0.25,0.5};
1344 int src_cell_mask[4] = {1,1, 0,1};
1345
1346 for (int with_mask = 0; with_mask < 2; ++with_mask) {
1347
1348 struct yac_basic_grid_data grid_data;
1349 if (is_tgt) {
1350 grid_data =
1351 yac_generate_basic_grid_data_reg_2d_deg(
1352 tgt_num_vertices, tgt_cyclic,
1353 tgt_coordinates_x, tgt_coordinates_y);
1354 } else {
1355 grid_data =
1356 yac_generate_basic_grid_data_unstruct_deg(
1357 src_num_vertices, src_num_cells, src_num_vertices_per_cell,
1358 src_coordinates_x, src_coordinates_y, src_cell_to_vertex);
1359 }
1360
1361 struct yac_basic_grid * grids[2] =
1362 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
1363 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
1364
1365 if (with_mask && !is_tgt) {
1366 int * mask = xmalloc(grid_data.num_vertices * sizeof(*mask));
1367 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1368 mask[i] = src_cell_mask[i];
1369 yac_basic_grid_add_mask_nocpy(grids[0], YAC_LOC_CORNER, mask, NULL);
1370 }
1371
1372 struct yac_dist_grid_pair * grid_pair =
1373 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
1374
1375 struct yac_interp_field src_fields[] =
1376 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX,
1377 .masks_idx = (with_mask)?0:SIZE_MAX}};
1378 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1379 struct yac_interp_field tgt_field =
1380 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX,
1381 .masks_idx = SIZE_MAX};
1382
1383 struct yac_interp_grid * interp_grid =
1384 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
1385 num_src_fields, src_fields, tgt_field);
1386
1387 for (int partial_coverage = 0; partial_coverage < 2; ++partial_coverage) {
1388
1389 struct interp_method * method_stack[] =
1390 {yac_interp_method_avg_new(YAC_INTERP_AVG_BARY, partial_coverage),
1391 yac_interp_method_fixed_new(-1.0), NULL};
1392
1393 struct yac_interp_weights * weights =
1394 yac_interp_method_do_search(method_stack, interp_grid);
1395
1396 for (size_t i = 0; i < num_reorder_types; ++i) {
1397
1398 struct yac_interpolation * interpolation =
1399 yac_interp_weights_get_interpolation(
1400 weights, reorder_types[i], 1,
1401 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1402
1403 // check generated interpolation
1404 {
1405 double * src_field = NULL;
1406 double ** src_fields = &src_field;
1407 double * tgt_field = NULL;
1408 double * ref_tgt_field = NULL;
1409 double ref_global_tgt_field[25];
1410 double ref_weights[2][2][25][4] =
1411 {{{{1.00, 0.00, 0.00, 0.00},
1412 {0.75, 0.25, 0.00, 0.00},
1413 {0.50, 0.50, 0.00, 0.00},
1414 {0.25, 0.75, 0.00, 0.00},
1415 {0.00, 1.00, 0.00, 0.00},
1416 {0.75, 0.00, 0.25, 0.00},
1417 {0.75, 0.00, 0.00, 0.25},
1418 {0.50, 0.25, 0.00, 0.25},
1419 {0.25, 0.50, 0.00, 0.25},
1420 {0.00, 0.75, 0.00, 0.25},
1421 {0.50, 0.00, 0.50, 0.00},
1422 {0.50, 0.00, 0.25, 0.25},
1423 {0.50, 0.00, 0.00, 0.50},
1424 {0.25, 0.25, 0.00, 0.50},
1425 {0.00, 0.50, 0.00, 0.50},
1426 {0.25, 0.00, 0.75, 0.00},
1427 {0.25, 0.00, 0.50, 0.25},
1428 {0.25, 0.00, 0.25, 0.50},
1429 {0.25, 0.00, 0.00, 0.75},
1430 {0.00, 0.25, 0.00, 0.75},
1431 {0.00, 0.00, 1.00, 0.00},
1432 {0.00, 0.00, 0.75, 0.25},
1433 {0.00, 0.00, 0.50, 0.50},
1434 {0.00, 0.00, 0.25, 0.75},
1435 {0.00, 0.00, 0.00, 1.00}},
1436 {{1.00, 0.00, 0.00, 0.00},
1437 {0.75, 0.25, 0.00, 0.00},
1438 {0.50, 0.50, 0.00, 0.00},
1439 {0.25, 0.75, 0.00, 0.00},
1440 {0.00, 1.00, 0.00, 0.00},
1441 {0.00, 0.00, 0.00, 0.00},//{0.75, 0.00, 0.25, 0.00}
1442 {0.75, 0.00, 0.00, 0.25},
1443 {0.50, 0.25, 0.00, 0.25},
1444 {0.25, 0.50, 0.00, 0.25},
1445 {0.00, 0.75, 0.00, 0.25},
1446 {0.00, 0.00, 0.00, 0.00},//{0.50, 0.00, 0.50, 0.00}
1447 {0.00, 0.00, 0.00, 0.00},//{0.50, 0.00, 0.25, 0.25}
1448 {0.50, 0.00, 0.00, 0.50},
1449 {0.25, 0.25, 0.00, 0.50},
1450 {0.00, 0.50, 0.00, 0.50},
1451 {0.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.75, 0.00}
1452 {0.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.50, 0.25}
1453 {0.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.25, 0.50}
1454 {0.25, 0.00, 0.00, 0.75},
1455 {0.00, 0.25, 0.00, 0.75},
1456 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 1.00, 0.00}
1457 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 0.75, 0.25}
1458 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 0.50, 0.50}
1459 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 0.25, 0.75}
1460 {0.00, 0.00, 0.00, 1.00}}},
1461 {{{1.00, 0.00, 0.00, 0.00},
1462 {0.75, 0.25, 0.00, 0.00},
1463 {0.50, 0.50, 0.00, 0.00},
1464 {0.25, 0.75, 0.00, 0.00},
1465 {0.00, 1.00, 0.00, 0.00},
1466 {0.75, 0.00, 0.25, 0.00},
1467 {0.75, 0.00, 0.00, 0.25},
1468 {0.50, 0.25, 0.00, 0.25},
1469 {0.25, 0.50, 0.00, 0.25},
1470 {0.00, 0.75, 0.00, 0.25},
1471 {0.50, 0.00, 0.50, 0.00},
1472 {0.50, 0.00, 0.25, 0.25},
1473 {0.50, 0.00, 0.00, 0.50},
1474 {0.25, 0.25, 0.00, 0.50},
1475 {0.00, 0.50, 0.00, 0.50},
1476 {0.25, 0.00, 0.75, 0.00},
1477 {0.25, 0.00, 0.50, 0.25},
1478 {0.25, 0.00, 0.25, 0.50},
1479 {0.25, 0.00, 0.00, 0.75},
1480 {0.00, 0.25, 0.00, 0.75},
1481 {0.00, 0.00, 1.00, 0.00},
1482 {0.00, 0.00, 0.75, 0.25},
1483 {0.00, 0.00, 0.50, 0.50},
1484 {0.00, 0.00, 0.25, 0.75},
1485 {0.00, 0.00, 0.00, 1.00}},
1486 {{1.00, 0.00, 0.00, 0.00},
1487 {0.75, 0.25, 0.00, 0.00},
1488 {0.50, 0.50, 0.00, 0.00},
1489 {0.25, 0.75, 0.00, 0.00},
1490 {0.00, 1.00, 0.00, 0.00},
1491 {1.00, 0.00, 0.00, 0.00},//{0.75, 0.00, 0.25, 0.00}
1492 {0.75, 0.00, 0.00, 0.25},
1493 {0.50, 0.25, 0.00, 0.25},
1494 {0.25, 0.50, 0.00, 0.25},
1495 {0.00, 0.75, 0.00, 0.25},
1496 {1.00, 0.00, 0.00, 0.00},//{0.50, 0.00, 0.50, 0.00}
1497 {0.50/0.75, 0.00, 0.00, 0.25/0.75},//{0.50, 0.00, 0.25, 0.25}
1498 {0.50, 0.00, 0.00, 0.50},
1499 {0.25, 0.25, 0.00, 0.50},
1500 {0.00, 0.50, 0.00, 0.50},
1501 {1.00, 0.00, 0.00, 0.00},//{0.25, 0.00, 0.75, 0.00}
1502 {0.50, 0.00, 0.00, 0.50},//{0.25, 0.00, 0.50, 0.25}
1503 {0.25/0.75, 0.00, 0.00, 0.50/0.75},//{0.25, 0.00, 0.25, 0.50}
1504 {0.25, 0.00, 0.00, 0.75},
1505 {0.00, 0.25, 0.00, 0.75},
1506 {0.00, 0.00, 0.00, 0.00},//{0.00, 0.00, 1.00, 0.00}
1507 {0.00, 0.00, 0.00, 1.00},//{0.00, 0.00, 0.75, 0.25}
1508 {0.00, 0.00, 0.00, 1.00},//{0.00, 0.00, 0.50, 0.50}
1509 {0.00, 0.00, 0.00, 1.00},//{0.00, 0.00, 0.25, 0.75}
1510 {0.00, 0.00, 0.00, 1.00}}}};
1511
1512 for (size_t j = 0; j < 25; ++j) {
1513 ref_global_tgt_field[j] =
1514 1.0 * ref_weights[partial_coverage][with_mask][j][0] +
1515 2.0 * ref_weights[partial_coverage][with_mask][j][1] +
1516 3.0 * ref_weights[partial_coverage][with_mask][j][2] +
1517 4.0 * ref_weights[partial_coverage][with_mask][j][3];
1518 if (ref_global_tgt_field[j] == 0.0)
1519 ref_global_tgt_field[j] = -1.0;
1520 }
1521
1522 if (is_tgt) {
1523 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
1524 ref_tgt_field = xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
1525 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1526 ref_tgt_field[i] = ref_global_tgt_field[i];
1527 } else {
1528 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
1529 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1530 src_field[i] = (double)(i) + 1.0;
1531 }
1532
1533 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
1534
1535 if (is_tgt)
1536 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1537 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-3)
1538 PUT_ERR("wrong interpolation result");
1539
1540 free(src_field);
1541 free(tgt_field);
1542 free(ref_tgt_field);
1543 }
1544
1545 yac_interpolation_delete(interpolation);
1546 }
1547
1548 yac_interp_weights_delete(weights);
1549 yac_interp_method_delete(method_stack);
1550 } // partial_coverage
1551 yac_interp_grid_delete(interp_grid);
1552 yac_dist_grid_pair_delete(grid_pair);
1553 yac_basic_grid_delete(grids[1]);
1554 yac_basic_grid_delete(grids[0]);
1555 }
1556 }
1557
1558
1559
1560 {// Test 9
1561 // simple interpolation test with various configuration options
1562
1563 // the global source grid:
1564 // 06--10--07--11--08
1565 // | | |
1566 // 06 02 08 03 09
1567 // | | |
1568 // 03--05--04--07--05
1569 // | | |
1570 // 01 00 03 01 04
1571 // | | |
1572 // 00--00--01--02--02
1573 //
1574 //---------------
1575 // setup
1576 //---------------
1577
1578 int is_tgt = split_comm_size == 1;
1579 double coordinates_x[2][7] = {{0.0,0.1,0.2}, {0.0,0.05,0.1,0.15,0.2}};
1580 double coordinates_y[2][7] = {{0.0,0.1,0.2}, {0.0,0.05,0.1,0.15,0.2}};
1581 size_t const num_cells[2][2] = {{2,2}, {4,4}};
1582 size_t local_start[2][2][2] = {{{0,0},{0,1}}, {{0,0}}};
1583 size_t local_count[2][2][2] = {{{2,1},{2,1}}, {{4,4}}};
1584 int global_mask[2][5*5] = {{1,1,1,
1585 0,1,1,
1586 0,0,1},
1587 {1,1,1,1,1,
1588 1,1,1,1,0,
1589 1,1,1,0,0,
1590 1,1,0,0,0,
1591 1,0,0,0,0}};
1592 int with_halo = 0;
1593 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
1594 coordinates_x[is_tgt][i] *= YAC_RAD;
1595 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
1596 coordinates_y[is_tgt][i] *= YAC_RAD;
1597
1598 struct yac_basic_grid_data grid_data =
1599 yac_generate_basic_grid_data_reg2d(
1600 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
1601 local_start[is_tgt][split_comm_rank],
1602 local_count[is_tgt][split_comm_rank], with_halo);
1603
1604 struct yac_basic_grid * grids[2] =
1605 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
1606 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
1607
1608 int field_mask[5*5];
1609 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1610 field_mask[i] = global_mask[is_tgt][grid_data.vertex_ids[i]];
1611 yac_basic_grid_add_mask(
1612 grids[0], YAC_LOC_CORNER, field_mask, grid_data.num_vertices, NULL);
1613
1614 struct yac_dist_grid_pair * grid_pair =
1615 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
1616
1617 for (int with_src_mask = 0; with_src_mask < 2; ++with_src_mask) {
1618 for (int with_tgt_mask = 0; with_tgt_mask < 2; ++with_tgt_mask) {
1619
1620 struct yac_interp_field src_fields[] =
1621 {{.location = YAC_LOC_CORNER,
1622 .coordinates_idx = SIZE_MAX,
1623 .masks_idx = (with_src_mask)?0:SIZE_MAX}};
1624 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1625 struct yac_interp_field tgt_field =
1626 {.location = YAC_LOC_CORNER,
1627 .coordinates_idx = SIZE_MAX,
1628 .masks_idx = (with_tgt_mask)?0:SIZE_MAX};
1629
1630 struct yac_interp_grid * interp_grid =
1631 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
1632 num_src_fields, src_fields, tgt_field);
1633
1634 enum yac_interp_avg_weight_type weight_types[] =
1635 {YAC_INTERP_AVG_ARITHMETIC,
1636 YAC_INTERP_AVG_DIST,
1637 YAC_INTERP_AVG_BARY};
1638 size_t const num_weight_types =
1639 sizeof(weight_types) / sizeof(weight_types[0]);
1640
1641 for (size_t weight_types_idx = 0; weight_types_idx < num_weight_types;
1642 ++weight_types_idx) {
1643 for (int partial_coverage = 0; partial_coverage < 2;
1644 ++partial_coverage) {
1645
1646 struct interp_method * method_stack[] =
1647 {yac_interp_method_avg_new(
1648 weight_types[weight_types_idx], partial_coverage),
1649 yac_interp_method_fixed_new(-1.0), NULL};
1650
1651 struct yac_interp_weights * weights =
1652 yac_interp_method_do_search(method_stack, interp_grid);
1653
1654 for (size_t i = 0; i < num_reorder_types; ++i) {
1655
1656 struct yac_interpolation * interpolation =
1657 yac_interp_weights_get_interpolation(
1658 weights, reorder_types[i], 1,
1659 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1660
1661 // check generated interpolation
1662 {
1663 double src_field_data[3*3];
1664 double * src_field = src_field_data;
1665 double ** src_fields = &src_field;
1666 double tgt_field_data[5*5] =
1667 {-2.0,-2.0,-2.0,-2.0,-2.0,
1668 -2.0,-2.0,-2.0,-2.0,-2.0,
1669 -2.0,-2.0,-2.0,-2.0,-2.0,
1670 -2.0,-2.0,-2.0,-2.0,-2.0,
1671 -2.0,-2.0,-2.0,-2.0,-2.0};
1672 double * tgt_field = tgt_field_data;
1673
1674 if (!is_tgt)
1675 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1676 src_field[i] = (double)(grid_data.vertex_ids[i]);
1677
1678 yac_interpolation_execute(
1679 interpolation, &src_fields, &tgt_field);
1680
1681 if (is_tgt) {
1682 size_t src_per_tgt[5*5][4] =
1683 {{0,1,3,4},{0,1,3,4},{0,1,3,4},{1,2,4,5},{1,2,4,5},
1684 {0,1,3,4},{0,1,3,4},{0,1,3,4},{1,2,4,5},{1,2,4,5},
1685 {0,1,3,4},{0,1,3,4},{0,1,3,4},{1,2,4,5},{1,2,4,5},
1686 {3,4,6,7},{3,4,6,7},{3,4,6,7},{4,5,7,8},{4,5,7,8},
1687 {3,4,6,7},{3,4,6,7},{3,4,6,7},{4,5,7,8},{4,5,7,8}};
1688 double A = 0.345491502812526;
1689 double B = 0.154508497187473;
1690 double C = 0.25;
1691 double D = 0.5;
1692 double weights[3][5*5][4] =
1693 {{{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},
1694 {C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},
1695 {C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},
1696 {C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},
1697 {C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C},{C,C,C,C}},
1698 {{1,0,0,0},{A,A,B,B},{0,1,0,0},{A,A,B,B},{0,1,0,0},
1699 {A,B,A,B},{C,C,C,C},{B,A,B,A},{C,C,C,C},{B,A,B,A},
1700 {0,0,1,0},{B,B,A,A},{0,0,0,1},{B,B,A,A},{0,0,0,1},
1701 {A,B,A,B},{C,C,C,C},{B,A,B,A},{C,C,C,C},{B,A,B,A},
1702 {0,0,1,0},{B,B,A,A},{0,0,0,1},{B,B,A,A},{0,0,0,1}},
1703 {{1,0,0,0},{D,D,0,0},{0,1,0,0},{D,D,0,0},{0,1,0,0},
1704 {D,0,D,0},{D,0,0,D},{0,D,0,D},{D,0,0,D},{0,D,0,D},
1705 {0,0,1,0},{0,0,D,D},{0,0,0,1},{0,0,D,D},{0,0,0,1},
1706 {D,0,D,0},{D,0,0,D},{0,D,0,D},{D,0,0,D},{0,D,0,D},
1707 {0,0,1,0},{0,0,D,D},{0,0,0,1},{0,0,D,D},{0,0,0,1}}};
1708
1709 for (size_t i = 0; i < grid_data.num_vertices; ++i) {
1710
1711 size_t tgt_idx = grid_data.vertex_ids[i];
1712 double src_values[4];
1713 double w[4];
1714
1715 double ref_tgt_value = 0.0;
1716
1717 int contains_masked_src = 0;
1718 int contains_unmasked_src = 0;
1719
1720 for (int j = 0; j < 4; ++j) {
1721 size_t src_idx = src_per_tgt[tgt_idx][j];
1722 w[j] = weights[weight_types_idx][tgt_idx][j];
1723 if (with_src_mask &&
1724 !global_mask[0][src_idx] &&
1725 (w[j] != 0.0)) {
1726 src_values[j] = -1.0;
1727 w[j] = 0.0;
1728 contains_masked_src = 1;
1729 } else {
1730 src_values[j] = (double)src_idx;
1731 contains_unmasked_src |= w[j] != 0.0;
1732 }
1733 }
1734
1735 if (with_tgt_mask && !global_mask[1][tgt_idx]) {
1736 ref_tgt_value = -2.0;
1737 } else if (contains_masked_src) {
1738 if (partial_coverage && contains_unmasked_src) {
1739 yac_correct_weights(4, w);
1740 for (int j = 0; j < 4; ++j)
1741 ref_tgt_value += w[j] * src_values[j];
1742 } else {
1743 ref_tgt_value = -1.0;
1744 }
1745 } else {
1746 for (int j = 0; j < 4; ++j)
1747 ref_tgt_value += w[j] * src_values[j];
1748 }
1749
1750 // some special cases
1751 if (with_src_mask && partial_coverage &&
1752 (weight_types[weight_types_idx] ==
1753 YAC_INTERP_AVG_DIST) &&
1754 (!with_tgt_mask || global_mask[1][tgt_idx])) {
1755 switch (tgt_idx) {
1756 case(10):
1757 ref_tgt_value = 1.73879612503;
1758 break;
1759 case(20):
1760 case(22):
1761 ref_tgt_value = 4.0;
1762 default:
1763 break;
1764 }
1765 }
1766
1767 if (fabs(ref_tgt_value - tgt_field[i]) > 1e-3)
1768 PUT_ERR("wrong interpolation result");
1769 }
1770 }
1771 }
1772
1773 yac_interpolation_delete(interpolation);
1774 } // reorder type
1775 yac_interp_weights_delete(weights);
1776 yac_interp_method_delete(method_stack);
1777 } // partial coverage
1778 } // weight type
1779 yac_interp_grid_delete(interp_grid);
1780 } // with tgt mask
1781 } // with src mask
1782 yac_dist_grid_pair_delete(grid_pair);
1783 yac_basic_grid_delete(grids[1]);
1784 yac_basic_grid_delete(grids[0]);
1785 }
1786
1787 {// Test 10
1788 // barycentric coordinates weighted (cell touches pole and has duplicated
1789 // longitude coordinate at pole)
1790 // interpolation with two source processes and a single target process
1791
1792 // the global source grid is a 1x1 grid:
1793 // 02--03--03
1794 // | |
1795 // 01 00 02
1796 // | |
1797 // 00--00--01
1798 //
1799 //---------------
1800 // setup
1801 //---------------
1802
1803 int is_tgt = split_comm_size == 1;
1804 double src_coordinates_x[2][2] = {{-0.5,0.5}, {0.5,-0.5}};
1805 double src_coordinates_y[2][2] = {{89,90}, {90,89}};
1806 size_t const src_num_cells[2] = {1,1};
1807 size_t src_local_start[2] = {0,0};
1808 size_t src_local_count[2] = {1,1};
1809 yac_int src_reorder[2][2][4] =
1810 {{{0,1,2,3},{2,3,0,1}},{{1,0,3,2},{3,2,1,0}}};
1811 double tgt_coordinates_x[5] = {-0.5,-0.25,0.0,0.25,0.5};
1812 double tgt_coordinates_y[5] = {89.00,89.25,89.50,89.75,90.00};
1813 size_t const tgt_num_cells[2] = {4,4};
1814 size_t tgt_local_start[2] = {0,0};
1815 size_t tgt_local_count[2] = {4,4};
1816
1817 for (size_t i = 0; i < 2 * 2; ++i)
1818 (&(src_coordinates_x[0][0]))[i] *= YAC_RAD;
1819 for (size_t i = 0; i < 2 * 2; ++i)
1820 (&(src_coordinates_y[0][0]))[i] *= YAC_RAD;
1821 for (size_t i = 0; i < 5; ++i) tgt_coordinates_x[i] *= YAC_RAD;
1822 for (size_t i = 0; i < 5; ++i)
1823 (&(tgt_coordinates_y[0]))[i] *= YAC_RAD;
1824
1825 for (int order_x = 0; order_x < 2; ++order_x) {
1826 for (int order_y = 0; order_y < 2; ++order_y) {
1827
1828 struct yac_basic_grid_data grid_data;
1829
1830 if (is_tgt) {
1831 grid_data =
1832 yac_generate_basic_grid_data_reg2d(
1833 tgt_coordinates_x, tgt_coordinates_y,
1834 tgt_num_cells, tgt_local_start, tgt_local_count, 0);
1835 } else {
1836 grid_data =
1837 yac_generate_basic_grid_data_reg2d(
1838 src_coordinates_x[order_x],
1839 src_coordinates_y[order_y], src_num_cells,
1840 src_local_start, src_local_count, 0);
1841 }
1842
1843 if (!is_tgt) {
1844
1845 yac_int * vertex_ids =
1846 xmalloc(grid_data.num_vertices * sizeof(*vertex_ids));
1847 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1848 vertex_ids[i] = (yac_int)(src_reorder[order_x][order_y][i]);
1849 grid_data.vertex_ids = vertex_ids;
1850 }
1851
1852 struct yac_basic_grid * grids[2] =
1853 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
1854 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
1855
1856 if (!is_tgt) {
1857
1858 yac_coordinate_pointer vertex_coordinates =
1859 xmalloc(grid_data.num_vertices * sizeof(*vertex_coordinates));
1860 memcpy(
1861 vertex_coordinates, grid_data.vertex_coordinates,
1862 grid_data.num_vertices * sizeof(*vertex_coordinates));
1863 // set x and y coordinate of all pole points to 0.0 in order to
1864 // reproduce an error that occured in ICON
1865 for (size_t i = 0; i < grid_data.num_vertices; ++i) {
1866 if (fabs(1.0 - fabs(vertex_coordinates[i][2])) < 1e-9) {
1867 vertex_coordinates[i][0] = 0.0;
1868 vertex_coordinates[i][1] = 0.0;
1869 }
1870 }
1871 yac_basic_grid_add_coordinates_nocpy(
1872 grids[0], YAC_LOC_CORNER, vertex_coordinates);
1873 }
1874
1875 struct yac_dist_grid_pair * grid_pair =
1876 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
1877
1878 struct yac_interp_field src_fields[] =
1879 {{.location = YAC_LOC_CORNER, .coordinates_idx = 0,
1880 .masks_idx = SIZE_MAX}};
1881 size_t num_src_fields =
1882 sizeof(src_fields) / sizeof(src_fields[0]);
1883 struct yac_interp_field tgt_field =
1884 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX,
1885 .masks_idx = SIZE_MAX};
1886
1887 struct yac_interp_grid * interp_grid =
1888 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
1889 num_src_fields, src_fields, tgt_field);
1890
1891 int partial_coverage = 0;
1892 struct interp_method * method_stack[] =
1893 {yac_interp_method_avg_new(
1894 YAC_INTERP_AVG_BARY, partial_coverage),
1895 yac_interp_method_fixed_new(-1.0), NULL};
1896
1897 struct yac_interp_weights * weights =
1898 yac_interp_method_do_search(method_stack, interp_grid);
1899
1900 for (size_t i = 0; i < num_reorder_types; ++i) {
1901
1902 struct yac_interpolation * interpolation =
1903 yac_interp_weights_get_interpolation(
1904 weights, reorder_types[i], 1,
1905 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1906
1907 // check generated interpolation
1908 {
1909 double * src_field = NULL;
1910 double ** src_fields = &src_field;
1911 double * tgt_field = NULL;
1912 double * ref_tgt_field = NULL;
1913 double ref_global_tgt_field[25];
1914 double ref_weights[25][4] =
1915 {{1.00, 0.00, 0.00, 0.00},
1916 {0.75, 0.25, 0.00, 0.00},
1917 {0.50, 0.50, 0.00, 0.00},
1918 {0.25, 0.75, 0.00, 0.00},
1919 {0.00, 1.00, 0.00, 0.00},
1920 {0.75, 0.00, 0.25, 0.00},
1921 {0.75, 0.00, 0.00, 0.25},
1922 {0.50, 0.25, 0.00, 0.25},
1923 {0.25, 0.50, 0.00, 0.25},
1924 {0.00, 0.75, 0.00, 0.25},
1925 {0.50, 0.00, 0.50, 0.00},
1926 {0.50, 0.00, 0.25, 0.25},
1927 {0.50, 0.00, 0.00, 0.50},
1928 {0.25, 0.25, 0.00, 0.50},
1929 {0.00, 0.50, 0.00, 0.50},
1930 {0.25, 0.00, 0.75, 0.00},
1931 {0.25, 0.00, 0.50, 0.25},
1932 {0.25, 0.00, 0.25, 0.50},
1933 {0.25, 0.00, 0.00, 0.75},
1934 {0.00, 0.25, 0.00, 0.75},
1935 {0.00, 0.00, 0.50, 0.50},
1936 {0.00, 0.00, 0.50, 0.50},
1937 {0.00, 0.00, 0.50, 0.50},
1938 {0.00, 0.00, 0.50, 0.50},
1939 {0.00, 0.00, 0.50, 0.50}};
1940
1941 for (size_t j = 0; j < 25; ++j) {
1942 ref_global_tgt_field[j] =
1943 (double)(src_reorder[order_x][order_y][0] + 1) *
1944 ref_weights[j][0] +
1945 (double)(src_reorder[order_x][order_y][1] + 1) *
1946 ref_weights[j][1] +
1947 (double)(src_reorder[order_x][order_y][2] + 1) *
1948 ref_weights[j][2] +
1949 (double)(src_reorder[order_x][order_y][3] + 1) *
1950 ref_weights[j][3];
1951 if (ref_global_tgt_field[j] == 0.0)
1952 ref_global_tgt_field[j] = -1.0;
1953 }
1954
1955 if (is_tgt) {
1956 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
1957 ref_tgt_field = xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
1958 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1959 ref_tgt_field[i] = ref_global_tgt_field[i];
1960 } else {
1961 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
1962 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1963 src_field[i] = (double)(i) + 1.0;
1964 }
1965
1966 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
1967
1968 if (is_tgt)
1969 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1970 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-3)
1971 PUT_ERR("wrong interpolation result");
1972
1973 free(src_field);
1974 free(tgt_field);
1975 free(ref_tgt_field);
1976 }
1977
1978 yac_interpolation_delete(interpolation);
1979 }
1980
1981 yac_interp_weights_delete(weights);
1982 yac_interp_method_delete(method_stack);
1983 yac_interp_grid_delete(interp_grid);
1984 yac_dist_grid_pair_delete(grid_pair);
1985 yac_basic_grid_delete(grids[1]);
1986 yac_basic_grid_delete(grids[0]);
1987 }
1988 }
1989 }
1990
1991 MPI_Comm_free(&split_comm);
1992
1993 xt_finalize();
1994
1995 MPI_Finalize();
1996
1997 return TEST_EXIT_CODE;
1998}
yac_basic_grid_new
struct yac_basic_grid * yac_basic_grid_new(char const *name, struct yac_basic_grid_data grid_data)
Definition basic_grid.c:50
yac_basic_grid_add_mask
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)
Definition basic_grid.c:284
yac_basic_grid_add_coordinates_nocpy
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:208
yac_basic_grid_empty_new
struct yac_basic_grid * yac_basic_grid_empty_new(char const *name)
Definition basic_grid.c:63
yac_basic_grid_delete
void yac_basic_grid_delete(struct yac_basic_grid *grid)
Definition basic_grid.c:70
yac_basic_grid_add_mask_nocpy
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:258
yac_generate_basic_grid_data_unstruct_deg
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)
Definition grid_unstruct.c:275
yac_generate_basic_grid_data_unstruct
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)
Definition grid_unstruct.c:265
yac_generate_basic_grid_data_reg_2d_deg
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)
Definition grid_reg2d.c:74
yac_correct_weights
void yac_correct_weights(size_t nSourceCells, double *weight)
correct interpolation weights
Definition clipping.c:1445
yac_dist_grid_pair_delete
void yac_dist_grid_pair_delete(struct yac_dist_grid_pair *grid_pair)
Definition dist_grid.c:2313
yac_dist_grid_pair_new
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:2061
yac_generate_basic_grid_data_reg2d
struct yac_basic_grid_data yac_generate_basic_grid_data_reg2d(double const *global_coords_x, double const *global_coords_y, size_t const num_global_cells_[2], size_t const local_start[2], size_t const local_count[2], int with_halo)
Definition dist_grid_utils.c:200
dist_grid_utils.h
YAC_RAD
#define YAC_RAD
Definition dynamic_config.c:12
LLtoXYZ_deg
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:269
normalise_vector
static void normalise_vector(double v[])
Definition geometry.h:635
yac_interp_grid_delete
void yac_interp_grid_delete(struct yac_interp_grid *interp_grid)
Definition interp_grid.c:575
yac_interp_grid_new
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
yac_interp_method_delete
void yac_interp_method_delete(struct interp_method **method)
Definition interp_method.c:54
yac_interp_method_do_search
struct yac_interp_weights * yac_interp_method_do_search(struct interp_method **method, struct yac_interp_grid *interp_grid)
Definition interp_method.c:15
interp_method.h
yac_interp_method_avg_new
struct interp_method * yac_interp_method_avg_new(enum yac_interp_avg_weight_type weight_type, int partial_coverage)
Definition interp_method_avg.c:933
interp_method_avg.h
yac_interp_avg_weight_type
yac_interp_avg_weight_type
Definition interp_method_avg.h:12
YAC_INTERP_AVG_DIST
@ YAC_INTERP_AVG_DIST
Definition interp_method_avg.h:14
YAC_INTERP_AVG_ARITHMETIC
@ YAC_INTERP_AVG_ARITHMETIC
Definition interp_method_avg.h:13
YAC_INTERP_AVG_BARY
@ YAC_INTERP_AVG_BARY
Definition interp_method_avg.h:15
yac_interp_method_fixed_new
struct interp_method * yac_interp_method_fixed_new(double value)
Definition interp_method_fixed.c:52
interp_method_fixed.h
yac_interp_weights_get_interpolation
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, int is_source, int is_target)
Definition interp_weights.c:4835
yac_interp_weights_delete
void yac_interp_weights_delete(struct yac_interp_weights *weights)
Definition interp_weights.c:6414
yac_interp_weights_reorder_type
yac_interp_weights_reorder_type
Definition interp_weights.h:42
YAC_MAPPING_ON_TGT
@ YAC_MAPPING_ON_TGT
weights will be applied at target processes
Definition interp_weights.h:44
YAC_MAPPING_ON_SRC
@ YAC_MAPPING_ON_SRC
weights will be applied at source processes
Definition interp_weights.h:43
yac_interpolation_execute
void yac_interpolation_execute(struct yac_interpolation *interp, double ***src_fields, double **tgt_field)
Execute interpolation synchronously and write results to the target field.
Definition interpolation.c:251
yac_interpolation_delete
void yac_interpolation_delete(struct yac_interpolation *interp)
Free an interpolation object and release all resources.
Definition interpolation.c:371
YAC_FRAC_MASK_NO_VALUE
double const YAC_FRAC_MASK_NO_VALUE
Definition interpolation.c:25
YAC_LOC_CORNER
@ YAC_LOC_CORNER
Definition location.h:15
YAC_LOC_CELL
@ YAC_LOC_CELL
Definition location.h:14
test_core.weights
weights
Definition test_core.py:63
test_core.src_field
src_field
Definition test_core.py:55
test_raw_exchange.i
int i
Definition test_raw_exchange.py:89
test_raw_exchange.tgt_field
tgt_field
Definition test_raw_exchange.py:64
yac_generate_weights.src_mask
src_mask
Definition yac_generate_weights.py:108
xmalloc
#define xmalloc(size)
Definition ppm_xfuncs.h:66
grid_data
Definition test_dummy_coupling3_c.c:34
grid_data::num_vertices_per_cell
int * num_vertices_per_cell
Definition test_dummy_coupling3_c.c:40
grid_data::cell_to_vertex
int * cell_to_vertex
Definition test_dummy_coupling3_c.c:41
interp_method
Definition interp_method_internal.h:36
yac_basic_grid_data
Definition basic_grid_data.h:13
yac_basic_grid_data::vertex_ids
yac_int * vertex_ids
Definition basic_grid_data.h:16
yac_basic_grid
Definition basic_grid.c:17
yac_dist_grid_pair
Definition dist_grid.c:119
yac_interp_field
Definition basic_grid.h:15
yac_interp_field::location
enum yac_location location
Definition basic_grid.h:16
yac_interp_grid
Definition interp_grid.c:22
yac_interp_grid::tgt_field
struct yac_interp_field tgt_field
Definition interp_grid.c:26
yac_interp_grid::num_src_fields
size_t num_src_fields
Definition interp_grid.c:27
yac_interp_grid::grid_pair
struct yac_dist_grid_pair * grid_pair
Definition interp_grid.c:25
yac_interp_grid::src_fields
struct yac_interp_field src_fields[]
Definition interp_grid.c:28
yac_interp_weights
Definition interp_weights.c:116
yac_interpolation
Definition interpolation.c:28
test_common.h
grid_names
static char const * grid_names[2]
Definition test_interp_method_avg_parallel.c:30
num_reorder_types
size_t num_reorder_types
Definition test_interp_method_avg_parallel.c:29
reorder_types
enum yac_interp_weights_reorder_type reorder_types[]
Definition test_interp_method_avg_parallel.c:27
split_comm
static MPI_Comm split_comm
Definition test_interp_method_file_parallel.c:56
num_weight_types
size_t num_weight_types
Definition test_interp_method_spmap_parallel.c:32
weight_types
enum yac_interp_spmap_weight_type weight_types[]
Definition test_interp_method_spmap_parallel.c:28
tgt_field_data
static double tgt_field_data[MAX_COLLECTION_SIZE][NUM_TGT_POINTS]
Definition test_interp_operator_fixed.c:30
coordinates_x
double coordinates_x[]
Definition test_interpolation_parallel3.c:104
num_cells
size_t num_cells[2]
Definition test_interpolation_parallel3.c:103
cell_coordinates_y
double cell_coordinates_y[]
Definition test_interpolation_parallel3.c:107
coordinates_y
double coordinates_y[]
Definition test_interpolation_parallel3.c:105
ref_weights
double ref_weights[4 *36]
Definition test_interpolation_parallel3.c:63
with_halo
int with_halo
Definition test_interpolation_parallel3.c:109
cell_coordinates_x
double cell_coordinates_x[]
Definition test_interpolation_parallel3.c:106
mask
static int mask[16]
Definition test_read_icon_common.h:34
TEST_EXIT_CODE
#define TEST_EXIT_CODE
Definition tests.h:14
PUT_ERR
#define PUT_ERR(string)
Definition tests.h:10
grids
struct yac_basic_grid ** grids
Definition yac.c:152
yac_mpi.h
yac_int
YAC_INT yac_int
Definition yac_types.h:15
yac_coordinate_pointer
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19