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