YAC 3.21.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_interp_method_nnn_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 <stdio.h>
7#include <math.h>
8
9#include "tests.h"
13#include "dist_grid_utils.h"
14#include "yac_mpi.h"
15#include "geometry.h"
16
17#include <mpi.h>
18#include <yaxt.h>
19
25char const src_grid_name[] = "src_grid";
26char const tgt_grid_name[] = "tgt_grid";
27
28static void utest_target_main(MPI_Comm target_comm);
29static void utest_source_main(MPI_Comm source_comm);
30
31int main (void) {
32
33 MPI_Init(NULL, NULL);
34
35 xt_initialize(MPI_COMM_WORLD);
36
37 int comm_rank, comm_size;
38 MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
39 MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
40 MPI_Barrier(MPI_COMM_WORLD);
41
42 if (comm_size != 5) {
43 PUT_ERR("ERROR: wrong number of processes");
44 xt_finalize();
45 MPI_Finalize();
46 return TEST_EXIT_CODE;
47 }
48
49 int tgt_flag = comm_rank < 1;
50
51 MPI_Comm split_comm;
52 MPI_Comm_split(MPI_COMM_WORLD, tgt_flag, 0, &split_comm);
53
54 if (tgt_flag) utest_target_main(split_comm);
55 else utest_source_main(split_comm);
56
57 MPI_Comm_free(&split_comm);
58 xt_finalize();
59 MPI_Finalize();
60
61 return TEST_EXIT_CODE;
62}
63
64static void utest_source_main(MPI_Comm source_comm) {
65
66 { // 1 and 5 nearest neighbours per target point with fixed value fallback
67
68 // corner and cell ids for a 7 x 7 grid (x == target point position)
69
70 // 56-----57-----58-----59-----60-----61-----62-----63
71 // | x| x| x| x| x| x| x|
72 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
73 // | | | | | | | |
74 // 48-----49-----50-----51-----52-----53-----54-----55
75 // | x| x| x| x| x| x| x|
76 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
77 // | | | | | | | |
78 // 40-----41-----42-----43-----44-----45-----46-----47
79 // | x| x| x| x| x| x| x|
80 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
81 // | | | | | | | |
82 // 32-----33-----34-----35-----36-----37-----38-----39
83 // | x| x| x| x| x| x| x|
84 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
85 // | | | | | | | |
86 // 24-----25-----26-----27-----28-----29-----30-----31
87 // | x| x| x| x| x| x| x|
88 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
89 // | | | | | | | |
90 // 16-----17-----18-----19-----20-----21-----22-----23
91 // | x| x| x| x| x| x| x|
92 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
93 // | | | | | | | |
94 // 08-----09-----10-----11-----12-----13-----14-----15
95 // | x| x| x| x| x| x| x|
96 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
97 // | | | | | | | |
98 // 00-----01-----02-----03-----04-----05-----06-----07
99 //
100 // the grid is distributed among the processes as follows:
101 // (index == process)
102 //
103 // 3---3---3---3---3---3---3---3
104 // | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
105 // 3---3---3---3---3---3---3---3
106 // | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
107 // 3---3---3---1---2---2---3---3
108 // | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
109 // 1---1---1---2---2---2---2---2
110 // | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
111 // 1---1---1---1---2---2---2---2
112 // | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
113 // 1---1---1---0---0---0---2---2
114 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
115 // 0---0---0---0---0---0---0---0
116 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
117 // 0---0---0---0---0---0---0---0
118 //
119 // the source mask looks as follows (# == masked point)
120 //
121 // +---+---+---+---+---+---+---#
122 // | | | | | | | |
123 // +---+---+---+---+---+---#---+
124 // | | | | | | | |
125 // +---+---+---+---+---#---+---+
126 // | | | | | | | |
127 // +---+---+---+---#---+---+---+
128 // | | | | | | | |
129 // +---+---+---#---+---+---+---+
130 // | | | | | | | |
131 // +---+---#---+---+---+---+---+
132 // | | | | | | | |
133 // +---#---+---+---+---+---+---+
134 // | | | | | | | |
135 // #---+---+---+---+---+---+---+
136
137 int my_source_rank;
138 MPI_Comm_rank(source_comm, &my_source_rank);
139
140 double coordinates_x[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
141 double coordinates_y[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
142 size_t const num_cells[2] = {7,7};
143 size_t local_start[4][2] = {{0,0},{0,2},{3,2},{0,5}};
144 size_t local_count[4][2] = {{7,2},{3,3},{4,3},{7,2}};
145 int with_halo = 0;
146 int global_corner_mask[8][8] = {
147 {0,1,1,1,1,1,1,1},
148 {1,0,1,1,1,1,1,1},
149 {1,1,0,1,1,1,1,1},
150 {1,1,1,0,1,1,1,1},
151 {1,1,1,1,0,1,1,1},
152 {1,1,1,1,1,0,1,1},
153 {1,1,1,1,1,1,0,1},
154 {1,1,1,1,1,1,1,0}};
155 for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
156 for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
157
159 utest_generate_basic_grid_data_reg2d(
161 local_start[my_source_rank], local_count[my_source_rank], with_halo);
162
163 int * src_corner_mask =
164 xmalloc(grid_data.num_vertices * sizeof(*src_corner_mask));
165 for (size_t i = 0; i < grid_data.num_vertices; ++i)
166 src_corner_mask[i] =
167 ((int*)(&(global_corner_mask[0][0])))[grid_data.vertex_ids[i]];
168
169 struct yac_basic_grid * src_grid =
172 src_grid, YAC_LOC_CORNER, src_corner_mask, NULL);
174
175 struct yac_dist_grid_pair * grid_pair =
176 yac_dist_grid_pair_new(src_grid, tgt_grid, MPI_COMM_WORLD);
177
178 struct yac_interp_field src_fields[] =
179 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
180 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
182 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX};
183
184 struct yac_interp_grid * interp_grid =
187
188 struct interp_method * method_stack[][3] = {
191 .type = YAC_INTERP_NNN_AVG, .n = 1,
193 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
197 .type = YAC_INTERP_NNN_AVG, .n = 5,
199 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
203 .type = YAC_INTERP_NNN_DIST, .n = 5,
205 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
209 .type = YAC_INTERP_NNN_GAUSS, .n = 5,
211 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT,
212 .data.gauss_scale = YAC_INTERP_NNN_GAUSS_SCALE_DEFAULT}),
216 .type = YAC_INTERP_NNN_ZERO, .n = 1,
218 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
222 .type = YAC_INTERP_NNN_AVG, .n = 1,
224 .max_search_distance = 0.4 * YAC_RAD}),
228 .type = YAC_INTERP_NNN_AVG, .n = 5,
230 .max_search_distance = 1.3 * YAC_RAD}),
232 // n_min < n: only 2 or 3 of the n source points are within the search
233 // distance; all target points are interpolated using the source points
234 // that were found
237 .type = YAC_INTERP_NNN_AVG, .n = 5, .n_min = 2,
238 .max_search_distance = 0.85 * YAC_RAD}),
240 // same as above, but with n_min == 3; target points that have only 2
241 // source points within the search distance are not interpolated
244 .type = YAC_INTERP_NNN_AVG, .n = 5, .n_min = 3,
245 .max_search_distance = 0.85 * YAC_RAD}),
246 yac_interp_method_fixed_new(-1), NULL}};
247 enum {NUM_STACKS = sizeof(method_stack) / sizeof(method_stack[0])};
248
249 struct yac_interp_weights * weights[NUM_STACKS];
250 for (size_t i = 0; i < NUM_STACKS; ++i) {
251 weights[i] = yac_interp_method_do_search(method_stack[i], interp_grid);
252 yac_interp_method_delete(method_stack[i]);
253 }
254
255 enum yac_interp_weights_reorder_type reorder_type[2] =
257
258 for (size_t i = 0; i < NUM_STACKS; ++i) {
259 for (size_t j = 0; j < 2; ++j) {
260 for (size_t collection_size = 1; collection_size < 4;
261 collection_size += 2) {
262
263 struct yac_interpolation * interpolation =
265 weights[i], reorder_type[j], collection_size,
266 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
267
268 // check generated interpolation
269 {
270 double *** src_data = xmalloc(collection_size * sizeof(*src_data));
271
272 for (size_t collection_idx = 0; collection_idx < collection_size;
273 ++collection_idx) {
274 // only one field
275 src_data[collection_idx] = xmalloc(1 * sizeof(**src_data));
276 src_data[collection_idx][0] =
277 xmalloc(grid_data.num_vertices * sizeof(***src_data));
278 for (size_t i = 0; i < grid_data.num_vertices; ++i)
279 src_data[collection_idx][0][i] =
280 (double)(grid_data.vertex_ids[i]) +
281 (double)(collection_idx * 49);
282 }
283
284 yac_interpolation_execute(interpolation, src_data, NULL);
285
286 for (size_t collection_idx = 0; collection_idx < collection_size;
287 ++collection_idx) {
288 free(src_data[collection_idx][0]);
289 free(src_data[collection_idx]);
290 }
291 free(src_data);
292 }
293
294 yac_interpolation_delete(interpolation);
295 }
296 }
297 }
298
299 for (size_t i = 0; i < NUM_STACKS; ++i)
300 yac_interp_weights_delete(weights[i]);
301 yac_interp_grid_delete(interp_grid);
302 yac_dist_grid_pair_delete(grid_pair);
303 yac_basic_grid_delete(tgt_grid);
304 yac_basic_grid_delete(src_grid);
305 }
306
307 { // 1 and 5 nearest neighbours per target point with fixed value fallback
308
309 // corner and cell ids for a 7 x 7 grid (x == target point position)
310
311 // 56-----57-----58-----59-----60-----61-----62-----63
312 // | x| x| x| x| x| x| x|
313 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
314 // | | | | | | | |
315 // 48-----49-----50-----51-----52-----53-----54-----55
316 // | x| x| x| x| x| x| x|
317 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
318 // | | | | | | | |
319 // 40-----41-----42-----43-----44-----45-----46-----47
320 // | x| x| x| x| x| x| x|
321 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
322 // | | | | | | | |
323 // 32-----33-----34-----35-----36-----37-----38-----39
324 // | x| x| x| x| x| x| x|
325 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
326 // | | | | | | | |
327 // 24-----25-----26-----27-----28-----29-----30-----31
328 // | x| x| x| x| x| x| x|
329 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
330 // | | | | | | | |
331 // 16-----17-----18-----19-----20-----21-----22-----23
332 // | x| x| x| x| x| x| x|
333 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
334 // | | | | | | | |
335 // 08-----09-----10-----11-----12-----13-----14-----15
336 // | x| x| x| x| x| x| x|
337 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
338 // | | | | | | | |
339 // 00-----01-----02-----03-----04-----05-----06-----07
340 //
341 // the grid is distributed among the processes as follows:
342 // (index == process)
343 //
344 // 3---3---3---3---3---3---3---3
345 // | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
346 // 3---3---3---3---3---3---3---3
347 // | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
348 // 3---3---3---1---2---2---3---3
349 // | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
350 // 1---1---1---2---2---2---2---2
351 // | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
352 // 1---1---1---1---2---2---2---2
353 // | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
354 // 1---1---1---0---0---0---2---2
355 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
356 // 0---0---0---0---0---0---0---0
357 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
358 // 0---0---0---0---0---0---0---0
359 //
360 // the source mask looks as follows (# == masked point)
361 //
362 // +---+---#---#---#---#---#---#
363 // | | | | | | | |
364 // +---+---#---#---#---#---#---#
365 // | | | | | | | |
366 // #---#---#---#---#---#---#---#
367 // | | | | | | | |
368 // #---#---#---#---#---#---#---#
369 // | | | | | | | |
370 // #---#---#---#---#---#---#---#
371 // | | | | | | | |
372 // #---#---#---#---#---#---#---#
373 // | | | | | | | |
374 // +---#---#---#---#---#---#---+
375 // | | | | | | | |
376 // +---+---#---#---#---#---+---+
377 //
378 // the target mask looks as follows (# == masked point)
379 //
380 // +---+---+---+---+---+---+---+
381 // | | | | | | | |
382 // +---+---+---+---+---+---+---+
383 // | | | | | | | |
384 // +---+---#---#---#---#---+---+
385 // | | | | | | | |
386 // +---+---#---#---#---#---+---+
387 // | | | | | | | |
388 // +---+---#---#---#---#---+---+
389 // | | | | | | | |
390 // +---+---#---#---#---#---+---+
391 // | | | | | | | |
392 // +---+---+---+---+---+---+---+
393 // | | | | | | | |
394 // +---+---+---+---+---+---+---+
395 //--------------------
396 // set up global search
397 //--------------------
398
399 int my_source_rank;
400 MPI_Comm_rank(source_comm, &my_source_rank);
401
402 double coordinates_x[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
403 double coordinates_y[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
404 size_t const num_cells[2] = {7,7};
405 size_t local_start[4][2] = {{0,0},{0,2},{3,2},{0,5}};
406 size_t local_count[4][2] = {{7,2},{3,3},{4,3},{7,2}};
407 int with_halo = 0;
408 int global_corner_mask[8][8] = {
409 {1,1,0,0,0,0,1,1},
410 {1,0,0,0,0,0,0,1},
411 {0,0,0,0,0,0,0,0},
412 {0,0,0,0,0,0,0,0},
413 {0,0,0,0,0,0,0,0},
414 {0,0,0,0,0,0,0,0},
415 {1,1,0,0,0,0,0,0},
416 {1,1,0,0,0,0,0,0}};
417 for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
418 for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
419
421 utest_generate_basic_grid_data_reg2d(
423 local_start[my_source_rank], local_count[my_source_rank], with_halo);
424
425 int * src_corner_mask =
426 xmalloc(grid_data.num_vertices * sizeof(*src_corner_mask));
427 for (size_t i = 0; i < grid_data.num_vertices; ++i)
428 src_corner_mask[i] =
429 ((int*)(&(global_corner_mask[0][0])))[grid_data.vertex_ids[i]];
430
431 struct yac_basic_grid * src_grid =
434 src_grid, YAC_LOC_CORNER, src_corner_mask, NULL);
436
437 struct yac_dist_grid_pair * grid_pair =
438 yac_dist_grid_pair_new(src_grid, tgt_grid, MPI_COMM_WORLD);
439
440 struct yac_interp_field src_fields[] =
441 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
442 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
444 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
445
446 struct yac_interp_grid * interp_grid =
449
450 struct interp_method * method_stack[2] = {
453 .type = YAC_INTERP_NNN_AVG, .n = 4,
455 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
456 NULL};
457
458 struct yac_interp_weights * weights =
459 yac_interp_method_do_search(method_stack, interp_grid);
460
461 yac_interp_method_delete(method_stack);
462
463 enum yac_interp_weights_reorder_type reorder_type[2] =
465
466 for (size_t i = 0; i < 2; ++i) {
467
468 struct yac_interpolation * interpolation =
470 weights, reorder_type[i], 1,
471 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
472
473 // check generated interpolation
474 {
475 double * src_field = NULL;
476 double ** src_fields = &src_field;
477
478 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
479 for (size_t i = 0; i < grid_data.num_vertices; ++i)
480 src_field[i] = (double)(grid_data.vertex_ids[i]);
481
482 yac_interpolation_execute(interpolation, &src_fields, NULL);
483
484 free(src_field);
485 }
486
487 yac_interpolation_delete(interpolation);
488 }
489
491 yac_interp_grid_delete(interp_grid);
492 yac_dist_grid_pair_delete(grid_pair);
493 yac_basic_grid_delete(tgt_grid);
494 yac_basic_grid_delete(src_grid);
495 }
496
497 { // 4 nearest neighbours per target point with fixed value fallback
498 // with matching vertices between source and target
499
500 double coordinates_x[] = {0.0, 1.0, 2.0, 3.0};
501 double coordinates_y[] = {0.0, 1.0, 2.0, 3.0};
502 size_t const num_cells[2] = {3,3};
503 size_t local_start[2] = {0,0};
504 size_t local_count[2] = {3,3};
505 int with_halo = 0;
506 for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
507 for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
508
510 utest_generate_basic_grid_data_reg2d(
512 local_start, local_count, with_halo);
513
514 struct yac_basic_grid * src_grid =
517
518 struct yac_dist_grid_pair * grid_pair =
519 yac_dist_grid_pair_new(src_grid, tgt_grid, MPI_COMM_WORLD);
520
521 struct yac_interp_field src_fields[] =
522 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
523 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
525 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0};
526
527 struct yac_interp_grid * interp_grid =
530
531 enum {NUM_STACKS = 1};
532 struct interp_method * method_stack[NUM_STACKS][3] = {
535 .type = YAC_INTERP_NNN_DIST, .n = 4,
537 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
538 yac_interp_method_fixed_new(-1), NULL}};
539
540 struct yac_interp_weights * weights[NUM_STACKS];
541 for (size_t i = 0; i < NUM_STACKS; ++i) {
542 weights[i] = yac_interp_method_do_search(method_stack[i], interp_grid);
543 yac_interp_method_delete(method_stack[i]);
544 }
545
546 for (size_t i = 0; i < NUM_STACKS; ++i) {
547 for (size_t collection_size = 1; collection_size < 4;
548 collection_size += 2) {
549
550 struct yac_interpolation * interpolation =
553 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
554
555 // check generated interpolation
556 {
557 double *** src_data = xmalloc(collection_size * sizeof(*src_data));
558
559 for (size_t collection_idx = 0; collection_idx < collection_size;
560 ++collection_idx) {
561 // only one field
562 src_data[collection_idx] = xmalloc(1 * sizeof(**src_data));
563 src_data[collection_idx][0] =
564 xmalloc(grid_data.num_vertices * sizeof(***src_data));
565 for (size_t i = 0; i < grid_data.num_vertices; ++i)
566 src_data[collection_idx][0][i] =
567 (double)(grid_data.vertex_ids[i]) +
568 (double)(collection_idx * 16);
569 }
570
571 yac_interpolation_execute(interpolation, src_data, NULL);
572
573 for (size_t collection_idx = 0; collection_idx < collection_size;
574 ++collection_idx) {
575 free(src_data[collection_idx][0]);
576 free(src_data[collection_idx]);
577 }
578 free(src_data);
579 }
580
581 yac_interpolation_delete(interpolation);
582 }
583 }
584
585 for (size_t i = 0; i < NUM_STACKS; ++i)
586 yac_interp_weights_delete(weights[i]);
587 yac_interp_grid_delete(interp_grid);
588 yac_dist_grid_pair_delete(grid_pair);
589 yac_basic_grid_delete(tgt_grid);
590 yac_basic_grid_delete(src_grid);
591 }
592
593 { // 4 nearest neighbours per target point
594 // where the 4 source points are far away from the target point
595
596 double coordinates_x[] = {0.0, 1.0, 2.0, 3.0};
597 double coordinates_y[] = {0.0, 1.0, 2.0, 3.0};
598 size_t const num_cells[2] = {3,3};
599 size_t local_start[2] = {0,0};
600 size_t local_count[2] = {3,3};
601 int with_halo = 0;
602 for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
603 for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
604
606 utest_generate_basic_grid_data_reg2d(
608 local_start, local_count, with_halo);
609
610 struct yac_basic_grid * src_grid =
613
614 struct yac_dist_grid_pair * grid_pair =
615 yac_dist_grid_pair_new(src_grid, tgt_grid, MPI_COMM_WORLD);
616
617 struct yac_interp_field src_fields[] =
618 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
619 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
621 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0};
622
623 struct yac_interp_grid * interp_grid =
626
627 enum {NUM_STACKS = 2};
628 struct interp_method * method_stack[NUM_STACKS][3] = {
631 .type = YAC_INTERP_NNN_DIST, .n = 4,
633 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
637 .type = YAC_INTERP_NNN_GAUSS, .n = 4,
639 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT,
640 .data.gauss_scale = YAC_INTERP_NNN_GAUSS_SCALE_DEFAULT}),
641 yac_interp_method_fixed_new(-1), NULL}};
642
643 struct yac_interp_weights * weights[NUM_STACKS];
644 for (size_t i = 0; i < NUM_STACKS; ++i) {
645 weights[i] = yac_interp_method_do_search(method_stack[i], interp_grid);
646 yac_interp_method_delete(method_stack[i]);
647 }
648
649 for (size_t i = 0; i < NUM_STACKS; ++i) {
650 for (size_t collection_size = 1; collection_size < 4;
651 collection_size += 2) {
652
653 struct yac_interpolation * interpolation =
656 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
657
658 // check generated interpolation
659 {
660 double *** src_data = xmalloc(collection_size * sizeof(*src_data));
661
662 for (size_t collection_idx = 0; collection_idx < collection_size;
663 ++collection_idx) {
664 // only one field
665 src_data[collection_idx] = xmalloc(1 * sizeof(**src_data));
666 src_data[collection_idx][0] =
667 xmalloc(grid_data.num_vertices * sizeof(***src_data));
668 for (size_t i = 0; i < grid_data.num_vertices; ++i)
669 src_data[collection_idx][0][i] =
670 (double)(grid_data.vertex_ids[i]) +
671 (double)(collection_idx * 16);
672 }
673
674 yac_interpolation_execute(interpolation, src_data, NULL);
675
676 for (size_t collection_idx = 0; collection_idx < collection_size;
677 ++collection_idx) {
678 free(src_data[collection_idx][0]);
679 free(src_data[collection_idx]);
680 }
681 free(src_data);
682 }
683
684 yac_interpolation_delete(interpolation);
685 }
686 }
687
688 for (size_t i = 0; i < NUM_STACKS; ++i)
689 yac_interp_weights_delete(weights[i]);
690 yac_interp_grid_delete(interp_grid);
691 yac_dist_grid_pair_delete(grid_pair);
692 yac_basic_grid_delete(tgt_grid);
693 yac_basic_grid_delete(src_grid);
694 }
695
696 { // 4 nearest neighbours per target point
697 // (from edge to vertices)
698
699 double coordinates_x[] = {0.0, 1.0, 2.0, 3.0};
700 double coordinates_y[] = {0.0, 1.0, 2.0, 3.0};
701 size_t const num_cells[2] = {3,3};
702 size_t local_start[2] = {0,0};
703 size_t local_count[2] = {3,3};
704 int with_halo = 0;
705 double edge_coordinates_x[] = {0.5,0.0,1.5,1.0,2.5,2.0,3.0,
706 0.5,0.0,1.5,1.0,2.5,2.0,3.0,
707 0.5,0.0,1.5,1.0,2.5,2.0,3.0,
708 0.5,1.5,2.5};
709 double edge_coordinates_y[] = {0.0,0.5,0.0,0.5,0.0,0.5,0.5,
710 1.0,1.5,1.0,1.5,1.0,1.5,1.5,
711 2.0,2.5,2.0,2.5,2.0,2.5,2.5,
712 3.0,3.0,3.0};
713 enum {num_edges = sizeof(edge_coordinates_x) / sizeof(edge_coordinates_x[0])};
714 yac_coordinate_pointer edge_coordinates =
715 xmalloc(num_edges * sizeof(*edge_coordinates));
716 for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
717 for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
718 for (size_t i = 0; i < num_edges; ++i)
720 edge_coordinates_x[i], edge_coordinates_y[i], edge_coordinates[i]);
721
723 utest_generate_basic_grid_data_reg2d(
725 local_start, local_count, with_halo);
726
727 struct yac_basic_grid * src_grid =
730 src_grid, YAC_LOC_EDGE, edge_coordinates);
732
733 struct yac_dist_grid_pair * grid_pair =
734 yac_dist_grid_pair_new(src_grid, tgt_grid, MPI_COMM_WORLD);
735
736 struct yac_interp_field src_fields[] =
737 {{.location = YAC_LOC_EDGE, .coordinates_idx = 0, .masks_idx = SIZE_MAX}};
738 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
740 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
741
742 struct yac_interp_grid * interp_grid =
745
746 enum {NUM_STACKS = 1};
747 struct interp_method * method_stack[NUM_STACKS][3] = {
750 .type = YAC_INTERP_NNN_AVG, .n = 4,
752 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
753 yac_interp_method_fixed_new(-1), NULL}};
754
755 struct yac_interp_weights * weights[NUM_STACKS];
756 for (size_t i = 0; i < NUM_STACKS; ++i) {
757 weights[i] = yac_interp_method_do_search(method_stack[i], interp_grid);
758 yac_interp_method_delete(method_stack[i]);
759 }
760
761 for (size_t i = 0; i < NUM_STACKS; ++i) {
762 for (size_t collection_size = 1; collection_size <= 1;
763 ++collection_size) {
764
765 struct yac_interpolation * interpolation =
768 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
769
770 // check generated interpolation
771 {
772 double *** src_data = xmalloc(collection_size * sizeof(*src_data));
773
774 for (size_t collection_idx = 0; collection_idx < collection_size;
775 ++collection_idx) {
776 // only one field
777 src_data[collection_idx] = xmalloc(1 * sizeof(**src_data));
778 src_data[collection_idx][0] =
779 xmalloc(grid_data.num_edges * sizeof(***src_data));
780 for (size_t i = 0; i < grid_data.num_edges; ++i)
781 src_data[collection_idx][0][i] =
782 (double)(grid_data.edge_ids[i]) +
783 (double)(collection_idx * 16);
784 }
785
786 yac_interpolation_execute(interpolation, src_data, NULL);
787
788 for (size_t collection_idx = 0; collection_idx < collection_size;
789 ++collection_idx) {
790 free(src_data[collection_idx][0]);
791 free(src_data[collection_idx]);
792 }
793 free(src_data);
794 }
795
796 yac_interpolation_delete(interpolation);
797 }
798 }
799
800 for (size_t i = 0; i < NUM_STACKS; ++i)
801 yac_interp_weights_delete(weights[i]);
802 yac_interp_grid_delete(interp_grid);
803 yac_dist_grid_pair_delete(grid_pair);
804 yac_basic_grid_delete(tgt_grid);
805 yac_basic_grid_delete(src_grid);
806 }
807}
808
809static double utest_compute_dist_result_(
810 size_t n, size_t * src_indices, size_t tgt_index,
811 double * src_coordinates_x, double * src_coordinates_y, size_t src_size_x,
812 double * tgt_coordinates_x, double * tgt_coordinates_y, size_t tgt_size_x) {
813
814 double tgt_coord[3];
816 tgt_coordinates_x[tgt_index%tgt_size_x],
817 tgt_coordinates_y[tgt_index/tgt_size_x], tgt_coord);
818
819 double result = 0.0;
820 double weights_sum = 0.0;
821
822 for (size_t i = 0; i < n; ++i) {
823
824 double src_coord[3];
825 LLtoXYZ_deg(src_coordinates_x[src_indices[i]%src_size_x],
826 src_coordinates_y[src_indices[i]/src_size_x], src_coord);
827
828 double weight = 1.0 / get_vector_angle(tgt_coord, src_coord);
829 result += weight * (double)src_indices[i];
830 weights_sum += weight;
831 }
832 return result / weights_sum;
833}
834
835static double utest_compute_dist_result(
836 size_t n, size_t * src_indices, size_t tgt_index) {
837
838 return
839 utest_compute_dist_result_(
840 n, src_indices, tgt_index,
841 (double[]){0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},
842 (double[]){0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0}, 8,
843 (double[]){0.75,1.75,2.75,3.75,4.75,5.75,6.75},
844 (double[]){0.75,1.75,2.75,3.75,4.75,5.75,6.75}, 7);
845}
846
847static double utest_compute_dist_result_2(
848 size_t n, size_t * src_indices, size_t tgt_index) {
849
850 return
851 utest_compute_dist_result_(
852 n, src_indices, tgt_index,
853 (double[]){0.0,1.0,2.0,3.0},
854 (double[]){0.0,1.0,2.0,3.0}, 4,
855 (double[]){1.0,1.5,2.0},
856 (double[]){1.0,1.5,2.0}, 3);
857}
858
859static double utest_compute_dist_result_3(
860 size_t n, size_t * src_indices, size_t tgt_index) {
861
862 return
863 utest_compute_dist_result_(
864 n, src_indices, tgt_index,
865 (double[]){0.0,1.0,2.0,3.0},
866 (double[]){0.0,1.0,2.0,3.0}, 4,
867 (double[]){1.0+90.0,1.5+90.0,2.0+90.0},
868 (double[]){1.0,1.5,2.0}, 3);
869}
870
871static double utest_compute_gauss_result_(
872 size_t n, size_t * src_indices, size_t tgt_index,
873 double * src_coordinates_x, double * src_coordinates_y, size_t src_size_x,
874 double * tgt_coordinates_x, double * tgt_coordinates_y, size_t tgt_size_x) {
875
876 double src_coords[n][3];
877
878 for (size_t i = 0; i < n; ++i)
879 LLtoXYZ_deg(src_coordinates_x[src_indices[i]%src_size_x],
880 src_coordinates_y[src_indices[i]/src_size_x], src_coords[i]);
881
882 // "- n" because we do not count the diagonal
883 double src_distances_sum = 0.0;
884 for (size_t i = 0; i < n; ++i)
885 for (size_t j = 0; j < n; ++j)
886 src_distances_sum += get_vector_angle(src_coords[i], src_coords[j]);
887
888 size_t src_distance_count = n * n - n;
889
890 // compute mean distance
891 double src_distance_mean = src_distances_sum / (double)src_distance_count;
892 double src_distance_mean_sq = src_distance_mean * src_distance_mean;
893
894 // compute weights
895 double tgt_coord[3];
896 LLtoXYZ_deg(tgt_coordinates_x[tgt_index%tgt_size_x],
897 tgt_coordinates_y[tgt_index/tgt_size_x], tgt_coord);
898 double weights[n];
899 for (size_t i = 0; i < n; ++i) {
900 double tgt_distance = get_vector_angle(src_coords[i], tgt_coord);
901 weights[i] =
902 exp(-1.0 * (tgt_distance * tgt_distance) /
903 (YAC_INTERP_NNN_GAUSS_SCALE_DEFAULT * src_distance_mean_sq));
904 }
905
906 // compute sum of weights
907 double weights_sum = 0.0;
908 for (size_t i = 0; i < n; ++i) weights_sum += weights[i];
909 for (size_t i = 0; i < n; ++i) weights[i] /= weights_sum;
910
911 // compute interpolation results
912 double result = 0.0;
913 for (size_t i = 0; i < n; ++i)
914 result += weights[i] * (double)src_indices[i];
915
916 return result;
917}
918
919static double utest_compute_gauss_result(
920 size_t n, size_t * src_indices, size_t tgt_index) {
921
922 return
923 utest_compute_gauss_result_(
924 n, src_indices, tgt_index,
925 (double[]){0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},
926 (double[]){0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0}, 8,
927 (double[]){0.75,1.75,2.75,3.75,4.75,5.75,6.75},
928 (double[]){0.75,1.75,2.75,3.75,4.75,5.75,6.75}, 7);
929}
930
931static void utest_target_main(MPI_Comm target_comm) {
932
933 UNUSED(target_comm);
934
935 {
936 double coordinates_x[] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
937 double coordinates_y[] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
938 double cell_coordinates_x[] = {0.75,1.75,2.75,3.75,4.75,5.75,6.75};
939 double cell_coordinates_y[] = {0.75,1.75,2.75,3.75,4.75,5.75,6.75};
940 yac_coordinate_pointer cell_coordinates = xmalloc(49 * sizeof(*cell_coordinates));
941 size_t const num_cells[2] = {7,7};
942 size_t local_start[2] = {0,0};
943 size_t local_count[2] = {7,7};
944 int with_halo = 0;
945 for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
946 for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
947 for (size_t i = 0, k = 0; i < num_cells[1]; ++i)
948 for (size_t j = 0; j < num_cells[0]; ++j, ++k)
950 cell_coordinates_x[j], cell_coordinates_y[i], cell_coordinates[k]);
951
953 utest_generate_basic_grid_data_reg2d(
955 local_start, local_count, with_halo);
956
957 struct yac_basic_grid * tgt_grid =
960 tgt_grid, YAC_LOC_CELL, cell_coordinates);
962
963 struct yac_dist_grid_pair * grid_pair =
964 yac_dist_grid_pair_new(tgt_grid, src_grid, MPI_COMM_WORLD);
965
966 struct yac_interp_field src_fields[] =
967 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
968 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
970 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX};
971
972 struct yac_interp_grid * interp_grid =
975
976 struct interp_method * method_stack[][3] = {
979 .type = YAC_INTERP_NNN_AVG, .n = 1,
981 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
985 .type = YAC_INTERP_NNN_AVG, .n = 5,
987 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
991 .type = YAC_INTERP_NNN_DIST, .n = 5,
993 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
997 .type = YAC_INTERP_NNN_GAUSS, .n = 5,
999 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT,
1000 .data.gauss_scale = YAC_INTERP_NNN_GAUSS_SCALE_DEFAULT}),
1001 yac_interp_method_fixed_new(-1), NULL},
1004 .type = YAC_INTERP_NNN_ZERO, .n = 1,
1006 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
1007 yac_interp_method_fixed_new(-1), NULL},
1010 .type = YAC_INTERP_NNN_AVG, .n = 1,
1012 .max_search_distance = 0.4 * YAC_RAD}),
1013 yac_interp_method_fixed_new(-1), NULL},
1016 .type = YAC_INTERP_NNN_AVG, .n = 5,
1018 .max_search_distance = 1.3 * YAC_RAD}),
1019 yac_interp_method_fixed_new(-1), NULL},
1020 // n_min < n: only 2 or 3 of the n source points are within the search
1021 // distance; all target points are interpolated using the source points
1022 // that were found
1025 .type = YAC_INTERP_NNN_AVG, .n = 5, .n_min = 2,
1026 .max_search_distance = 0.85 * YAC_RAD}),
1027 yac_interp_method_fixed_new(-1), NULL},
1028 // same as above, but with n_min == 3; target points that have only 2
1029 // source points within the search distance are not interpolated
1032 .type = YAC_INTERP_NNN_AVG, .n = 5, .n_min = 3,
1033 .max_search_distance = 0.85 * YAC_RAD}),
1034 yac_interp_method_fixed_new(-1), NULL}};
1035 enum {NUM_STACKS = sizeof(method_stack) / sizeof(method_stack[0])};
1036
1037
1038 struct yac_interp_weights * weights[NUM_STACKS];
1039 for (size_t i = 0; i < NUM_STACKS; ++i) {
1040 weights[i] = yac_interp_method_do_search(method_stack[i], interp_grid);
1041 yac_interp_method_delete(method_stack[i]);
1042 }
1043
1044 enum yac_interp_weights_reorder_type reorder_type[2] =
1046
1047 for (size_t i = 0; i < NUM_STACKS; ++i) {
1048 for (size_t j = 0; j < 2; ++j) {
1049 for (size_t collection_size = 1; collection_size < 4;
1050 collection_size += 2) {
1051
1052 struct yac_interpolation * interpolation =
1054 weights[i], reorder_type[j], collection_size,
1055 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1056
1057 // check generated interpolation
1058 {
1059 double ref_tgt_field[][49] =
1060 {{8,10,11,12,13,14,15,
1061 17,17,19,20,21,22,23,
1062 25,26,26,28,29,30,31,
1063 33,34,35,35,37,38,39,
1064 41,42,43,44,44,46,47,
1065 49,50,51,52,53,53,55,
1066 57,58,59,60,61,62,62},
1067 {(1+8+10+16+17)/5.0,(1+2+10+11+17)/5.0,
1068 (2+3+10+11+12)/5.0,(3+4+11+12+13)/5.0,
1069 (4+5+12+13+14)/5.0,(5+6+13+14+15)/5.0,
1070 (6+7+14+15+23)/5.0,(8+16+17+24+25)/5.0,
1071 (10+17+19+25+26)/5.0,(10+11+19+20+26)/5.0,
1072 (11+12+19+20+21)/5.0,(12+13+20+21+22)/5.0,
1073 (13+14+21+22+23)/5.0,(14+15+22+23+31)/5.0,
1074 (16+17+24+25+26)/5.0,(17+19+25+26+34)/5.0,
1075 (19+20+26+28+35)/5.0,(19+20+21+28+29)/5.0,
1076 (20+21+28+29+30)/5.0,(21+22+29+30+31)/5.0,
1077 (22+23+30+31+39)/5.0,(24+25+32+33+34)/5.0,
1078 (25+26+33+34+35)/5.0,(26+28+34+35+43)/5.0,
1079 (28+29+35+37+44)/5.0,(28+29+30+37+38)/5.0,
1080 (29+30+37+38+39)/5.0,(30+31+38+39+47)/5.0,
1081 (32+33+40+41+42)/5.0,(33+34+41+42+43)/5.0,
1082 (34+35+42+43+44)/5.0,(35+37+43+44+52)/5.0,
1083 (37+38+44+46+53)/5.0,(37+38+39+46+47)/5.0,
1084 (38+39+46+47+55)/5.0,(40+41+48+49+50)/5.0,
1085 (41+42+49+50+51)/5.0,(42+43+50+51+52)/5.0,
1086 (43+44+51+52+53)/5.0,(44+46+52+53+61)/5.0,
1087 (46+47+53+55+62)/5.0,(46+47+53+55+62)/5.0,
1088 (48+49+56+57+58)/5.0,(49+50+57+58+59)/5.0,
1089 (50+51+58+59+60)/5.0,(51+52+59+60+61)/5.0,
1090 (52+53+60+61+62)/5.0,(53+55+60+61+62)/5.0,
1091 (47+53+55+61+62)/5.0},
1092 {utest_compute_dist_result(5, (size_t[]){1,8,10,16,17}, 0),
1093 utest_compute_dist_result(5, (size_t[]){1,2,10,11,17}, 1),
1094 utest_compute_dist_result(5, (size_t[]){2,3,10,11,12}, 2),
1095 utest_compute_dist_result(5, (size_t[]){3,4,11,12,13}, 3),
1096 utest_compute_dist_result(5, (size_t[]){4,5,12,13,14}, 4),
1097 utest_compute_dist_result(5, (size_t[]){5,6,13,14,15}, 5),
1098 utest_compute_dist_result(5, (size_t[]){6,7,14,15,23}, 6),
1099 utest_compute_dist_result(5, (size_t[]){8,16,17,24,25}, 7),
1100 utest_compute_dist_result(5, (size_t[]){10,17,19,25,26}, 8),
1101 utest_compute_dist_result(5, (size_t[]){10,11,19,20,26}, 9),
1102 utest_compute_dist_result(5, (size_t[]){11,12,19,20,21}, 10),
1103 utest_compute_dist_result(5, (size_t[]){12,13,20,21,22}, 11),
1104 utest_compute_dist_result(5, (size_t[]){13,14,21,22,23}, 12),
1105 utest_compute_dist_result(5, (size_t[]){14,15,22,23,31}, 13),
1106 utest_compute_dist_result(5, (size_t[]){16,17,24,25,26}, 14),
1107 utest_compute_dist_result(5, (size_t[]){17,19,25,26,34}, 15),
1108 utest_compute_dist_result(5, (size_t[]){19,20,26,28,35}, 16),
1109 utest_compute_dist_result(5, (size_t[]){19,20,21,28,29}, 17),
1110 utest_compute_dist_result(5, (size_t[]){20,21,28,29,30}, 18),
1111 utest_compute_dist_result(5, (size_t[]){21,22,29,30,31}, 19),
1112 utest_compute_dist_result(5, (size_t[]){22,23,30,31,39}, 20),
1113 utest_compute_dist_result(5, (size_t[]){24,25,32,33,34}, 21),
1114 utest_compute_dist_result(5, (size_t[]){25,26,33,34,35}, 22),
1115 utest_compute_dist_result(5, (size_t[]){26,28,34,35,43}, 23),
1116 utest_compute_dist_result(5, (size_t[]){28,29,35,37,44}, 24),
1117 utest_compute_dist_result(5, (size_t[]){28,29,30,37,38}, 25),
1118 utest_compute_dist_result(5, (size_t[]){29,30,37,38,39}, 26),
1119 utest_compute_dist_result(5, (size_t[]){30,31,38,39,47}, 27),
1120 utest_compute_dist_result(5, (size_t[]){32,33,40,41,42}, 28),
1121 utest_compute_dist_result(5, (size_t[]){33,34,41,42,43}, 29),
1122 utest_compute_dist_result(5, (size_t[]){34,35,42,43,44}, 30),
1123 utest_compute_dist_result(5, (size_t[]){35,37,43,44,52}, 31),
1124 utest_compute_dist_result(5, (size_t[]){37,38,44,46,53}, 32),
1125 utest_compute_dist_result(5, (size_t[]){37,38,39,46,47}, 33),
1126 utest_compute_dist_result(5, (size_t[]){38,39,46,47,55}, 34),
1127 utest_compute_dist_result(5, (size_t[]){40,41,48,49,50}, 35),
1128 utest_compute_dist_result(5, (size_t[]){41,42,49,50,51}, 36),
1129 utest_compute_dist_result(5, (size_t[]){42,43,50,51,52}, 37),
1130 utest_compute_dist_result(5, (size_t[]){43,44,51,52,53}, 38),
1131 utest_compute_dist_result(5, (size_t[]){44,46,52,53,61}, 39),
1132 utest_compute_dist_result(5, (size_t[]){46,47,53,55,62}, 40),
1133 utest_compute_dist_result(5, (size_t[]){46,47,53,55,62}, 41),
1134 utest_compute_dist_result(5, (size_t[]){48,49,56,57,58}, 42),
1135 utest_compute_dist_result(5, (size_t[]){49,50,57,58,59}, 43),
1136 utest_compute_dist_result(5, (size_t[]){50,51,58,59,60}, 44),
1137 utest_compute_dist_result(5, (size_t[]){51,52,59,60,61}, 45),
1138 utest_compute_dist_result(5, (size_t[]){52,53,60,61,62}, 46),
1139 utest_compute_dist_result(5, (size_t[]){53,55,60,61,62}, 47),
1140 utest_compute_dist_result(5, (size_t[]){47,53,55,61,62}, 48)},
1141 {utest_compute_gauss_result(5, (size_t[]){1,8,10,16,17}, 0),
1142 utest_compute_gauss_result(5, (size_t[]){1,2,10,11,17}, 1),
1143 utest_compute_gauss_result(5, (size_t[]){2,3,10,11,12}, 2),
1144 utest_compute_gauss_result(5, (size_t[]){3,4,11,12,13}, 3),
1145 utest_compute_gauss_result(5, (size_t[]){4,5,12,13,14}, 4),
1146 utest_compute_gauss_result(5, (size_t[]){5,6,13,14,15}, 5),
1147 utest_compute_gauss_result(5, (size_t[]){6,7,14,15,23}, 6),
1148 utest_compute_gauss_result(5, (size_t[]){8,16,17,24,25}, 7),
1149 utest_compute_gauss_result(5, (size_t[]){10,17,19,25,26}, 8),
1150 utest_compute_gauss_result(5, (size_t[]){10,11,19,20,26}, 9),
1151 utest_compute_gauss_result(5, (size_t[]){11,12,19,20,21}, 10),
1152 utest_compute_gauss_result(5, (size_t[]){12,13,20,21,22}, 11),
1153 utest_compute_gauss_result(5, (size_t[]){13,14,21,22,23}, 12),
1154 utest_compute_gauss_result(5, (size_t[]){14,15,22,23,31}, 13),
1155 utest_compute_gauss_result(5, (size_t[]){16,17,24,25,26}, 14),
1156 utest_compute_gauss_result(5, (size_t[]){17,19,25,26,34}, 15),
1157 utest_compute_gauss_result(5, (size_t[]){19,20,26,28,35}, 16),
1158 utest_compute_gauss_result(5, (size_t[]){19,20,21,28,29}, 17),
1159 utest_compute_gauss_result(5, (size_t[]){20,21,28,29,30}, 18),
1160 utest_compute_gauss_result(5, (size_t[]){21,22,29,30,31}, 19),
1161 utest_compute_gauss_result(5, (size_t[]){22,23,30,31,39}, 20),
1162 utest_compute_gauss_result(5, (size_t[]){24,25,32,33,34}, 21),
1163 utest_compute_gauss_result(5, (size_t[]){25,26,33,34,35}, 22),
1164 utest_compute_gauss_result(5, (size_t[]){26,28,34,35,43}, 23),
1165 utest_compute_gauss_result(5, (size_t[]){28,29,35,37,44}, 24),
1166 utest_compute_gauss_result(5, (size_t[]){28,29,30,37,38}, 25),
1167 utest_compute_gauss_result(5, (size_t[]){29,30,37,38,39}, 26),
1168 utest_compute_gauss_result(5, (size_t[]){30,31,38,39,47}, 27),
1169 utest_compute_gauss_result(5, (size_t[]){32,33,40,41,42}, 28),
1170 utest_compute_gauss_result(5, (size_t[]){33,34,41,42,43}, 29),
1171 utest_compute_gauss_result(5, (size_t[]){34,35,42,43,44}, 30),
1172 utest_compute_gauss_result(5, (size_t[]){35,37,43,44,52}, 31),
1173 utest_compute_gauss_result(5, (size_t[]){37,38,44,46,53}, 32),
1174 utest_compute_gauss_result(5, (size_t[]){37,38,39,46,47}, 33),
1175 utest_compute_gauss_result(5, (size_t[]){38,39,46,47,55}, 34),
1176 utest_compute_gauss_result(5, (size_t[]){40,41,48,49,50}, 35),
1177 utest_compute_gauss_result(5, (size_t[]){41,42,49,50,51}, 36),
1178 utest_compute_gauss_result(5, (size_t[]){42,43,50,51,52}, 37),
1179 utest_compute_gauss_result(5, (size_t[]){43,44,51,52,53}, 38),
1180 utest_compute_gauss_result(5, (size_t[]){44,46,52,53,61}, 39),
1181 utest_compute_gauss_result(5, (size_t[]){46,47,53,55,62}, 40),
1182 utest_compute_gauss_result(5, (size_t[]){46,47,53,55,62}, 41),
1183 utest_compute_gauss_result(5, (size_t[]){48,49,56,57,58}, 42),
1184 utest_compute_gauss_result(5, (size_t[]){49,50,57,58,59}, 43),
1185 utest_compute_gauss_result(5, (size_t[]){50,51,58,59,60}, 44),
1186 utest_compute_gauss_result(5, (size_t[]){51,52,59,60,61}, 45),
1187 utest_compute_gauss_result(5, (size_t[]){52,53,60,61,62}, 46),
1188 utest_compute_gauss_result(5, (size_t[]){53,55,60,61,62}, 47),
1189 utest_compute_gauss_result(5, (size_t[]){47,53,55,61,62}, 48)},
1190 {0,0,0,0,0,0,0,
1191 0,0,0,0,0,0,0,
1192 0,0,0,0,0,0,0,
1193 0,0,0,0,0,0,0,
1194 0,0,0,0,0,0,0,
1195 0,0,0,0,0,0,0,
1196 0,0,0,0,0,0,0},
1197 {-1,10,11,12,13,14,15,
1198 17,-1,19,20,21,22,23,
1199 25,26,-1,28,29,30,31,
1200 33,34,35,-1,37,38,39,
1201 41,42,43,44,-1,46,47,
1202 49,50,51,52,53,-1,55,
1203 57,58,59,60,61,62,-1},
1204 {-1.0,-1.0,
1205 (2+3+10+11+12)/5.0,(3+4+11+12+13)/5.0,
1206 (4+5+12+13+14)/5.0,(5+6+13+14+15)/5.0,
1207 (6+7+14+15+23)/5.0,-1.0,
1208 -1.0,-1.0,
1209 (11+12+19+20+21)/5.0,(12+13+20+21+22)/5.0,
1210 (13+14+21+22+23)/5.0,(14+15+22+23+31)/5.0,
1211 (16+17+24+25+26)/5.0,-1.0,
1212 -1.0,-1.0,
1213 (20+21+28+29+30)/5.0,(21+22+29+30+31)/5.0,
1214 (22+23+30+31+39)/5.0,(24+25+32+33+34)/5.0,
1215 (25+26+33+34+35)/5.0,-1.0,
1216 -1.0,-1.0,
1217 (29+30+37+38+39)/5.0,(30+31+38+39+47)/5.0,
1218 (32+33+40+41+42)/5.0,(33+34+41+42+43)/5.0,
1219 (34+35+42+43+44)/5.0,-1.0,
1220 -1.0,-1.0,
1221 (38+39+46+47+55)/5.0,(40+41+48+49+50)/5.0,
1222 (41+42+49+50+51)/5.0,(42+43+50+51+52)/5.0,
1223 (43+44+51+52+53)/5.0,-1.0,
1224 -1.0,-1.0,
1225 (48+49+56+57+58)/5.0,(49+50+57+58+59)/5.0,
1226 (50+51+58+59+60)/5.0,(51+52+59+60+61)/5.0,
1227 (52+53+60+61+62)/5.0,-1.0,
1228 -1.0},
1229 // n = 5, n_min = 2, max_search_distance = 0.85 degree
1230 // (only the three source corners at the upper right, upper left,
1231 // and lower right of a target point are within the search
1232 // distance; one of them is masked for target points close to
1233 // the diagonal, which leaves 2 source points)
1234 {(1+8)/2.0,(2+10)/2.0,(3+10+11)/3.0,
1235 (4+11+12)/3.0,(5+12+13)/3.0,(6+13+14)/3.0,
1236 (7+14+15)/3.0,
1237 (16+17)/2.0,(10+17)/2.0,(11+19)/2.0,
1238 (12+19+20)/3.0,(13+20+21)/3.0,(14+21+22)/3.0,
1239 (15+22+23)/3.0,
1240 (17+24+25)/3.0,(25+26)/2.0,(19+26)/2.0,
1241 (20+28)/2.0,(21+28+29)/3.0,(22+29+30)/3.0,
1242 (23+30+31)/3.0,
1243 (25+32+33)/3.0,(26+33+34)/3.0,(34+35)/2.0,
1244 (28+35)/2.0,(29+37)/2.0,(30+37+38)/3.0,
1245 (31+38+39)/3.0,
1246 (33+40+41)/3.0,(34+41+42)/3.0,(35+42+43)/3.0,
1247 (43+44)/2.0,(37+44)/2.0,(38+46)/2.0,
1248 (39+46+47)/3.0,
1249 (41+48+49)/3.0,(42+49+50)/3.0,(43+50+51)/3.0,
1250 (44+51+52)/3.0,(52+53)/2.0,(46+53)/2.0,
1251 (47+55)/2.0,
1252 (49+56+57)/3.0,(50+57+58)/3.0,(51+58+59)/3.0,
1253 (52+59+60)/3.0,(53+60+61)/3.0,(61+62)/2.0,
1254 (55+62)/2.0},
1255 // same as above, but with n_min = 3; all target points that
1256 // have only 2 source points within the search distance are
1257 // handled by the fixed value fallback
1258 {-1.0,-1.0,(3+10+11)/3.0,
1259 (4+11+12)/3.0,(5+12+13)/3.0,(6+13+14)/3.0,
1260 (7+14+15)/3.0,
1261 -1.0,-1.0,-1.0,
1262 (12+19+20)/3.0,(13+20+21)/3.0,(14+21+22)/3.0,
1263 (15+22+23)/3.0,
1264 (17+24+25)/3.0,-1.0,-1.0,
1265 -1.0,(21+28+29)/3.0,(22+29+30)/3.0,
1266 (23+30+31)/3.0,
1267 (25+32+33)/3.0,(26+33+34)/3.0,-1.0,
1268 -1.0,-1.0,(30+37+38)/3.0,
1269 (31+38+39)/3.0,
1270 (33+40+41)/3.0,(34+41+42)/3.0,(35+42+43)/3.0,
1271 -1.0,-1.0,-1.0,
1272 (39+46+47)/3.0,
1273 (41+48+49)/3.0,(42+49+50)/3.0,(43+50+51)/3.0,
1274 (44+51+52)/3.0,-1.0,-1.0,
1275 -1.0,
1276 (49+56+57)/3.0,(50+57+58)/3.0,(51+58+59)/3.0,
1277 (52+59+60)/3.0,(53+60+61)/3.0,-1.0,
1278 -1.0}};
1279
1280 double ** tgt_data = xmalloc(collection_size * sizeof(*tgt_data));
1281 for (size_t collection_idx = 0; collection_idx < collection_size;
1282 ++collection_idx)
1283 tgt_data[collection_idx] =
1284 xmalloc(grid_data.num_cells * sizeof(**tgt_data));
1285
1286 yac_interpolation_execute(interpolation, NULL, tgt_data);
1287
1288 for (size_t collection_idx = 0; collection_idx < collection_size;
1289 ++collection_idx)
1290 for (size_t j = 0, offset = collection_idx * 49;
1291 j < grid_data.num_cells; ++j) {
1292 if (ref_tgt_field[i][j] == -1.0) {
1293 if (tgt_data[collection_idx][j] != -1.0)
1294 PUT_ERR("wrong interpolation result (fallback)");
1295 } else {
1296 if ((tgt_data[collection_idx][j] != 0.0) &&
1297 (fabs((ref_tgt_field[i][j] + (double)offset) -
1298 tgt_data[collection_idx][j]) > 1e-9))
1299 PUT_ERR("wrong interpolation result");
1300 }
1301 }
1302
1303 for (size_t collection_idx = 0; collection_idx < collection_size;
1304 ++collection_idx)
1305 free(tgt_data[collection_idx]);
1306 free(tgt_data);
1307 }
1308
1309 yac_interpolation_delete(interpolation);
1310 }
1311 }
1312 }
1313
1314 for (size_t i = 0; i < NUM_STACKS; ++i)
1315 yac_interp_weights_delete(weights[i]);
1316 yac_interp_grid_delete(interp_grid);
1317 yac_dist_grid_pair_delete(grid_pair);
1318 yac_basic_grid_delete(src_grid);
1319 yac_basic_grid_delete(tgt_grid);
1320 }
1321
1322 {
1323 double coordinates_x[] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
1324 double coordinates_y[] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
1325 double cell_coordinates_x[] = {0.75,1.75,2.75,3.75,4.75,5.75,6.75};
1326 double cell_coordinates_y[] = {0.75,1.75,2.75,3.75,4.75,5.75,6.75};
1327 yac_coordinate_pointer cell_coordinates = xmalloc(49 * sizeof(*cell_coordinates));
1328 size_t const num_cells[2] = {7,7};
1329 size_t local_start[2] = {0,0};
1330 size_t local_count[2] = {7,7};
1331 int with_halo = 0;
1332 int global_cell_mask[7][7] = {
1333 {1,1,1,1,1,1,1},
1334 {1,1,1,1,1,1,1},
1335 {1,1,0,0,0,1,1},
1336 {1,1,0,0,0,1,1},
1337 {1,1,0,0,0,1,1},
1338 {1,1,1,1,1,1,1},
1339 {1,1,1,1,1,1,1}};
1340 for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
1341 for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
1342 for (size_t i = 0, k = 0; i < num_cells[1]; ++i)
1343 for (size_t j = 0; j < num_cells[0]; ++j, ++k)
1345 cell_coordinates_x[j], cell_coordinates_y[i], cell_coordinates[k]);
1346
1348 utest_generate_basic_grid_data_reg2d(
1350 local_start, local_count, with_halo);
1351
1352 int * tgt_cell_mask =
1353 xmalloc(grid_data.num_cells * sizeof(*tgt_cell_mask));
1354 for (size_t i = 0; i < grid_data.num_cells; ++i)
1355 tgt_cell_mask[i] =
1356 ((int*)(&(global_cell_mask[0][0])))[grid_data.cell_ids[i]];
1357
1358 struct yac_basic_grid * tgt_grid =
1361 tgt_grid, YAC_LOC_CELL, tgt_cell_mask, NULL);
1363 tgt_grid, YAC_LOC_CELL, cell_coordinates);
1365
1366 struct yac_dist_grid_pair * grid_pair =
1367 yac_dist_grid_pair_new(tgt_grid, src_grid, MPI_COMM_WORLD);
1368
1369 struct yac_interp_field src_fields[] =
1370 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
1371 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1373 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
1374
1375 struct yac_interp_grid * interp_grid =
1378
1379 struct interp_method * method_stack[2] = {
1382 .type = YAC_INTERP_NNN_AVG, .n = 4,
1384 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}), NULL};
1385
1386 struct yac_interp_weights * weights =
1387 yac_interp_method_do_search(method_stack, interp_grid);
1388
1389 yac_interp_method_delete(method_stack);
1390
1391 enum yac_interp_weights_reorder_type reorder_type[2] =
1393
1394 for (size_t i = 0; i < 2; ++i) {
1395
1396 struct yac_interpolation * interpolation =
1398 weights, reorder_type[i], 1,
1399 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1400
1401 // check generated interpolation
1402 {
1403 double * tgt_field = NULL;
1404 double ref_tgt_field[49] =
1405 {(0+1+8+49)/4.0,(0+1+8+6)/4.0,(0+1+8+6)/4.0,(1+6+7+15)/4.0,
1406 (1+6+7+15)/4.0,(1+6+7+15)/4.0,(1+6+7+15)/4.0,
1407 (0+1+8+49)/4.0,(0+1+8+49)/4.0,(0+1+8+6)/4.0,(1+6+7+15)/4.0,
1408 (1+6+7+15)/4.0,(1+6+7+15)/4.0,(1+6+7+15)/4.0,
1409 (8+1+0+49)/4.0,(8+1+0+49)/4.0,-1,-1,
1410 -1,(1+6+7+15)/4.0,(1+6+7+15)/4.0,
1411 (48+49+8+57)/4.0,(48+49+8+57)/4.0,-1,-1,
1412 -1,(6+7+15+49)/4.0,(6+7+15+49)/4.0,
1413 (48+49+56+57)/4.0,(48+49+56+57)/4.0,-1,-1,
1414 -1,(6+7+15+49)/4.0,(6+7+15+49)/4.0,
1415 (48+49+56+57)/4.0,(48+49+56+57)/4.0,(48+49+56+57)/4.0,(48+49+56+57)/4.0,
1416 (48+49+56+57)/4.0,(48+49+15+57)/4.0,(6+7+15+49)/4.0,
1417 (48+49+56+57)/4.0,(48+49+56+57)/4.0,(48+49+56+57)/4.0,(48+49+56+57)/4.0,
1418 (48+49+56+57)/4.0,(48+49+56+57)/4.0,(49+57+56+15)/4.0};
1419
1420 tgt_field = xmalloc(grid_data.num_cells * sizeof(*tgt_field));
1421
1422 for (size_t j = 0; j < grid_data.num_cells; ++j) tgt_field[j] = -1.0;
1423
1424 yac_interpolation_execute(interpolation, NULL, &tgt_field);
1425
1426 for (size_t j = 0; j < grid_data.num_cells; ++j)
1427 if (fabs(ref_tgt_field[j] - tgt_field[j]) > 1e-9)
1428 PUT_ERR("wrong interpolation result");
1429
1430 free(tgt_field);
1431 }
1432
1433 yac_interpolation_delete(interpolation);
1434 }
1435
1437 yac_interp_grid_delete(interp_grid);
1438 yac_dist_grid_pair_delete(grid_pair);
1439 yac_basic_grid_delete(src_grid);
1440 yac_basic_grid_delete(tgt_grid);
1441 }
1442
1443 {
1444 double coordinates_x[] = {1.0,1.5,2.0};
1445 double coordinates_y[] = {1.0,1.5,2.0};
1446 size_t const num_cells[2] = {2,2};
1447 size_t local_start[2] = {0,0};
1448 size_t local_count[2] = {2,2};
1449 int with_halo = 0;
1450 int global_corner_mask[3*3] = {1,0,1, 0,1,0, 1,0,1};
1451 for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
1452 for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
1453
1455 utest_generate_basic_grid_data_reg2d(
1457 local_start, local_count, with_halo);
1458
1459 int * tgt_corner_mask =
1460 xmalloc(grid_data.num_vertices * sizeof(*tgt_corner_mask));
1461 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1462 tgt_corner_mask[i] = global_corner_mask[grid_data.vertex_ids[i]];
1463
1464 struct yac_basic_grid * tgt_grid =
1467 tgt_grid, YAC_LOC_CORNER, tgt_corner_mask, NULL);
1469
1470 struct yac_dist_grid_pair * grid_pair =
1471 yac_dist_grid_pair_new(tgt_grid, src_grid, MPI_COMM_WORLD);
1472
1473 struct yac_interp_field src_fields[] =
1474 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
1475 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1477 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0};
1478
1479 struct yac_interp_grid * interp_grid =
1482
1483 enum {NUM_STACKS = 1};
1484 struct interp_method * method_stack[NUM_STACKS][3] = {
1487 .type = YAC_INTERP_NNN_DIST, .n = 4,
1489 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
1490 yac_interp_method_fixed_new(-1), NULL}};
1491
1492
1493 struct yac_interp_weights * weights[NUM_STACKS];
1494 for (size_t i = 0; i < NUM_STACKS; ++i) {
1495 weights[i] = yac_interp_method_do_search(method_stack[i], interp_grid);
1496 yac_interp_method_delete(method_stack[i]);
1497 }
1498
1499 for (size_t i = 0; i < NUM_STACKS; ++i) {
1500 for (size_t collection_size = 1; collection_size < 4;
1501 collection_size += 2) {
1502
1503 struct yac_interpolation * interpolation =
1506 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1507
1508 // check generated interpolation
1509 {
1510 double ref_tgt_field[][3*3] =
1511 {{5.0,-2.0,6.0,
1512 -2.0, utest_compute_dist_result_2(4, (size_t[]){5,6,9,10}, 4),-2.0,
1513 9.0,-2.0,10.0}};
1514
1515 double ** tgt_data = xmalloc(collection_size * sizeof(*tgt_data));
1516 for (size_t collection_idx = 0; collection_idx < collection_size;
1517 ++collection_idx) {
1518 tgt_data[collection_idx] =
1519 xmalloc(grid_data.num_vertices * sizeof(**tgt_data));
1520 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1521 tgt_data[collection_idx][i] = -2.0;
1522 }
1523
1524 yac_interpolation_execute(interpolation, NULL, tgt_data);
1525
1526 for (size_t collection_idx = 0; collection_idx < collection_size;
1527 ++collection_idx)
1528 for (size_t j = 0, offset = collection_idx * 16;
1529 j < grid_data.num_vertices; ++j) {
1530 if (ref_tgt_field[i][j] != -2.0) {
1531 if ((fabs(
1532 (ref_tgt_field[i][j] + (double)offset) -
1533 tgt_data[collection_idx][j]) > 1e-9))
1534 PUT_ERR("wrong interpolation result");
1535 } else {
1536 if (tgt_data[collection_idx][j] != -2.0)
1537 PUT_ERR("wrong interpolation result");
1538 }
1539 }
1540
1541 for (size_t collection_idx = 0; collection_idx < collection_size;
1542 ++collection_idx)
1543 free(tgt_data[collection_idx]);
1544 free(tgt_data);
1545 }
1546
1547 yac_interpolation_delete(interpolation);
1548 }
1549 }
1550
1551 for (size_t i = 0; i < NUM_STACKS; ++i)
1552 yac_interp_weights_delete(weights[i]);
1553 yac_interp_grid_delete(interp_grid);
1554 yac_dist_grid_pair_delete(grid_pair);
1555 yac_basic_grid_delete(src_grid);
1556 yac_basic_grid_delete(tgt_grid);
1557 }
1558
1559 {
1560 double coordinates_x[] = {1.0+90.0,1.5+90.0,2.0+90.0};
1561 double coordinates_y[] = {1.0,1.5,2.0};
1562 size_t const num_cells[2] = {2,2};
1563 size_t local_start[2] = {0,0};
1564 size_t local_count[2] = {2,2};
1565 int with_halo = 0;
1566 int global_corner_mask[3*3] = {1,0,1, 0,1,0, 1,0,1};
1567 for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
1568 for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
1569
1571 utest_generate_basic_grid_data_reg2d(
1573 local_start, local_count, with_halo);
1574
1575 int * tgt_corner_mask =
1576 xmalloc(grid_data.num_vertices * sizeof(*tgt_corner_mask));
1577 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1578 tgt_corner_mask[i] = global_corner_mask[grid_data.vertex_ids[i]];
1579
1580 struct yac_basic_grid * tgt_grid =
1583 tgt_grid, YAC_LOC_CORNER, tgt_corner_mask, NULL);
1585
1586 struct yac_dist_grid_pair * grid_pair =
1587 yac_dist_grid_pair_new(tgt_grid, src_grid, MPI_COMM_WORLD);
1588
1589 struct yac_interp_field src_fields[] =
1590 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
1591 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1593 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0};
1594
1595 struct yac_interp_grid * interp_grid =
1598
1599 enum {NUM_STACKS = 2};
1600 struct interp_method * method_stack[NUM_STACKS][3] = {
1603 .type = YAC_INTERP_NNN_DIST, .n = 4,
1605 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
1606 yac_interp_method_fixed_new(-1), NULL},
1609 .type = YAC_INTERP_NNN_GAUSS, .n = 4,
1611 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT,
1612 .data.gauss_scale = YAC_INTERP_NNN_GAUSS_SCALE_DEFAULT}),
1613 yac_interp_method_fixed_new(-1), NULL}};
1614
1615
1616 struct yac_interp_weights * weights[NUM_STACKS];
1617 for (size_t i = 0; i < NUM_STACKS; ++i) {
1618 weights[i] = yac_interp_method_do_search(method_stack[i], interp_grid);
1619 yac_interp_method_delete(method_stack[i]);
1620 }
1621
1622 for (size_t i = 0; i < NUM_STACKS; ++i) {
1623 for (size_t collection_size = 1; collection_size < 4;
1624 collection_size += 2) {
1625
1626 struct yac_interpolation * interpolation =
1629 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1630
1631 // check generated interpolation
1632 {
1633 double ref_tgt_field[][3*3] =
1634 {{utest_compute_dist_result_3(4, (size_t[]){3,7,11,15}, 0),-2.0,
1635 utest_compute_dist_result_3(4, (size_t[]){3,7,11,15}, 2),
1636 -2.0, utest_compute_dist_result_3(4, (size_t[]){3,7,11,15}, 4),-2.0,
1637 utest_compute_dist_result_3(4, (size_t[]){3,7,11,15}, 6),-2.0,
1638 utest_compute_dist_result_3(4, (size_t[]){3,7,11,15}, 8)},
1639 {utest_compute_dist_result_3(4, (size_t[]){3,7,11,15}, 0),-2.0,
1640 utest_compute_dist_result_3(4, (size_t[]){3,7,11,15}, 2),
1641 -2.0, utest_compute_dist_result_3(4, (size_t[]){3,7,11,15}, 4),-2.0,
1642 utest_compute_dist_result_3(4, (size_t[]){3,7,11,15}, 6),-2.0,
1643 utest_compute_dist_result_3(4, (size_t[]){3,7,11,15}, 8)}};
1644
1645 double ** tgt_data = xmalloc(collection_size * sizeof(*tgt_data));
1646 for (size_t collection_idx = 0; collection_idx < collection_size;
1647 ++collection_idx) {
1648 tgt_data[collection_idx] =
1649 xmalloc(grid_data.num_vertices * sizeof(**tgt_data));
1650 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1651 tgt_data[collection_idx][i] = -2.0;
1652 }
1653
1654 yac_interpolation_execute(interpolation, NULL, tgt_data);
1655
1656 for (size_t collection_idx = 0; collection_idx < collection_size;
1657 ++collection_idx)
1658 for (size_t j = 0, offset = collection_idx * 16;
1659 j < grid_data.num_vertices; ++j) {
1660 if (ref_tgt_field[i][j] != -2.0) {
1661 if ((fabs(
1662 (ref_tgt_field[i][j] + (double)offset) -
1663 tgt_data[collection_idx][j]) > 1e-9) ||
1664 (tgt_data[collection_idx][j] !=
1665 tgt_data[collection_idx][j]))
1666 PUT_ERR("wrong interpolation result");
1667 } else {
1668 if (tgt_data[collection_idx][j] != -2.0)
1669 PUT_ERR("wrong interpolation result");
1670 }
1671 }
1672
1673 for (size_t collection_idx = 0; collection_idx < collection_size;
1674 ++collection_idx)
1675 free(tgt_data[collection_idx]);
1676 free(tgt_data);
1677 }
1678
1679 yac_interpolation_delete(interpolation);
1680 }
1681 }
1682
1683 for (size_t i = 0; i < NUM_STACKS; ++i)
1684 yac_interp_weights_delete(weights[i]);
1685 yac_interp_grid_delete(interp_grid);
1686 yac_dist_grid_pair_delete(grid_pair);
1687 yac_basic_grid_delete(src_grid);
1688 yac_basic_grid_delete(tgt_grid);
1689 }
1690
1691 {
1692 double coordinates_x[] = {0.5,1.5,2.5};
1693 double coordinates_y[] = {0.5,1.5,2.5};
1694 size_t const num_cells[2] = {2,2};
1695 size_t local_start[2] = {0,0};
1696 size_t local_count[2] = {2,2};
1697 int with_halo = 0;
1698 for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
1699 for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
1700
1702 utest_generate_basic_grid_data_reg2d(
1704 local_start, local_count, with_halo);
1705
1706 struct yac_basic_grid * tgt_grid =
1709
1710 struct yac_dist_grid_pair * grid_pair =
1711 yac_dist_grid_pair_new(tgt_grid, src_grid, MPI_COMM_WORLD);
1712
1713 struct yac_interp_field src_fields[] =
1714 {{.location = YAC_LOC_EDGE, .coordinates_idx = 0, .masks_idx = SIZE_MAX}};
1715 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1717 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
1718
1719 struct yac_interp_grid * interp_grid =
1722
1723 enum {NUM_STACKS = 1};
1724 struct interp_method * method_stack[NUM_STACKS][3] = {
1727 .type = YAC_INTERP_NNN_AVG, .n = 4,
1729 .max_search_distance = YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT}),
1730 yac_interp_method_fixed_new(-1), NULL}};
1731
1732
1733 struct yac_interp_weights * weights[NUM_STACKS];
1734 for (size_t i = 0; i < NUM_STACKS; ++i) {
1735 weights[i] = yac_interp_method_do_search(method_stack[i], interp_grid);
1736 yac_interp_method_delete(method_stack[i]);
1737 }
1738
1739 for (size_t i = 0; i < NUM_STACKS; ++i) {
1740 for (size_t collection_size = 1; collection_size <= 1;
1741 collection_size += 2) {
1742
1743 struct yac_interpolation * interpolation =
1746 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1747
1748 // check generated interpolation
1749 {
1750 double ref_tgt_field[][3*3] =
1751 {{(0+3+7+1)*0.25,(2+5+9+3)*0.25,(4+6+11+5)*0.25,
1752 (7+10+14+8)*0.25,(9+12+16+10)*0.25,(11+13+18+12)*0.25,
1753 (14+17+21+15)*0.25,(16+19+22+17)*0.25,(18+20+23+19)*0.25}};
1754
1755 double ** tgt_data = xmalloc(collection_size * sizeof(*tgt_data));
1756 for (size_t collection_idx = 0; collection_idx < collection_size;
1757 ++collection_idx) {
1758 tgt_data[collection_idx] =
1759 xmalloc(grid_data.num_vertices * sizeof(**tgt_data));
1760 for (size_t i = 0; i < grid_data.num_vertices; ++i)
1761 tgt_data[collection_idx][i] = -2.0;
1762 }
1763
1764 yac_interpolation_execute(interpolation, NULL, tgt_data);
1765
1766 for (size_t collection_idx = 0; collection_idx < collection_size;
1767 ++collection_idx)
1768 for (size_t j = 0, offset = collection_idx * 16;
1769 j < grid_data.num_vertices; ++j) {
1770 if (ref_tgt_field[i][j] != -2.0) {
1771 if ((fabs(
1772 (ref_tgt_field[i][j] + (double)offset) -
1773 tgt_data[collection_idx][j]) > 1e-9) ||
1774 (tgt_data[collection_idx][j] !=
1775 tgt_data[collection_idx][j]))
1776 PUT_ERR("wrong interpolation result");
1777 } else {
1778 if (tgt_data[collection_idx][j] != -2.0)
1779 PUT_ERR("wrong interpolation result");
1780 }
1781 }
1782
1783 for (size_t collection_idx = 0; collection_idx < collection_size;
1784 ++collection_idx)
1785 free(tgt_data[collection_idx]);
1786 free(tgt_data);
1787 }
1788
1789 yac_interpolation_delete(interpolation);
1790 }
1791 }
1792
1793 for (size_t i = 0; i < NUM_STACKS; ++i)
1794 yac_interp_weights_delete(weights[i]);
1795 yac_interp_grid_delete(interp_grid);
1796 yac_dist_grid_pair_delete(grid_pair);
1797 yac_basic_grid_delete(src_grid);
1798 yac_basic_grid_delete(tgt_grid);
1799 }
1800}
struct yac_basic_grid * yac_basic_grid_new(char const *name, struct yac_basic_grid_data grid_data)
Definition basic_grid.c:57
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:202
struct yac_basic_grid * yac_basic_grid_empty_new(char const *name)
Definition basic_grid.c:70
void yac_basic_grid_delete(struct yac_basic_grid *grid)
Definition basic_grid.c:77
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:244
#define UNUSED(x)
Definition core.h:72
void yac_dist_grid_pair_delete(struct yac_dist_grid_pair *grid_pair)
Definition dist_grid.c:2377
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:2089
#define YAC_RAD
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:278
static double get_vector_angle(double const a[3], double const b[3])
Definition geometry.h:472
void yac_interp_grid_delete(struct yac_interp_grid *interp_grid)
struct yac_interp_grid * yac_interp_grid_new(struct yac_dist_grid_pair *grid_pair, char const *src_grid_name, char const *tgt_grid_name, size_t num_src_fields, struct yac_interp_field const *src_fields, struct yac_interp_field const tgt_field)
Definition interp_grid.c:31
void yac_interp_method_delete(struct interp_method **method)
Delete an interpolation stack and free its resources (but not the pointer array).
struct yac_interp_weights * yac_interp_method_do_search(struct interp_method **method, struct yac_interp_grid *interp_grid)
Perform weight computation using given interpolation stack and grid.
Defines the interface of the interpolation method "base class" in YAC.
struct interp_method * yac_interp_method_fixed_new(double value)
struct interp_method * yac_interp_method_nnn_new(struct yac_interp_method_nnn_config config)
@ YAC_INTERP_NNN_GAUSS
distance with Gauss weights of n source points
@ YAC_INTERP_NNN_AVG
average of n source points
@ YAC_INTERP_NNN_DIST
distance weighted average of n source points
@ YAC_INTERP_NNN_ZERO
all weights are set to zero
#define YAC_INTERP_NNN_MAX_SEARCH_DISTANCE_DEFAULT
#define YAC_INTERP_NNN_GAUSS_SCALE_DEFAULT
#define YAC_INTERP_NNN_N_MIN_DEFAULT
a value of 0 is interpreted as "use the value of n"
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)
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 applied at source processes
void yac_interpolation_execute(struct yac_interpolation *interp, double ***src_fields, double **tgt_field)
Execute interpolation synchronously and write results to the target field.
void yac_interpolation_delete(struct yac_interpolation *interp)
Free an interpolation object and release all resources.
double const YAC_FRAC_MASK_NO_VALUE
@ YAC_LOC_CORNER
Definition location.h:15
@ YAC_LOC_EDGE
Definition location.h:16
@ YAC_LOC_CELL
Definition location.h:14
#define xmalloc(size)
Definition ppm_xfuncs.h:66
enum yac_location location
Definition basic_grid.h:16
struct yac_interp_field tgt_field
Definition interp_grid.c:26
size_t num_src_fields
Definition interp_grid.c:27
struct yac_dist_grid_pair * grid_pair
Definition interp_grid.c:25
struct yac_interp_field src_fields[]
Definition interp_grid.c:28
int collection_size
static MPI_Comm split_comm
char const src_grid_name[]
char const tgt_grid_name[]
double coordinates_x[]
size_t num_cells[2]
double cell_coordinates_y[]
double coordinates_y[]
double cell_coordinates_x[]
#define TEST_EXIT_CODE
Definition tests.h:15
#define PUT_ERR(string)
Definition tests.h:10
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:21