YAC 3.12.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_interp_method_hcsbb_parallel.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#include <stdlib.h>
6#include <string.h>
7
8#include "tests.h"
9#include "test_common.h"
10#include "geometry.h"
14#include "dist_grid_utils.h"
15#include "yac_mpi.h"
16
17#include <mpi.h>
18#include <yaxt.h>
19#include <netcdf.h>
20
28enum {NUM_REORDER_TYPES = sizeof(reorder_types) / sizeof(reorder_types[0])};
29
30static void utest_compute_field_data_XYZ(
31 double (*vertices)[3], size_t num_vertices, double * f_v);
32static char const * grid_names[2] = {"src_grid", "tgt_grid"};
33
34int main(void) {
35
36 MPI_Init(NULL, NULL);
37
38 xt_initialize(MPI_COMM_WORLD);
39
40 int comm_rank, comm_size;
41 MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
42 MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
43 MPI_Barrier(MPI_COMM_WORLD);
44
45 if (comm_size != 6) {
46 PUT_ERR("ERROR: wrong number of processes");
47 xt_finalize();
48 MPI_Finalize();
49 return TEST_EXIT_CODE;
50 }
51
52 MPI_Comm split_comm;
53 MPI_Comm_split(MPI_COMM_WORLD, comm_rank < 1, 0, &split_comm);
54
55 int split_comm_rank, split_comm_size;
56 MPI_Comm_rank(split_comm, &split_comm_rank);
57 MPI_Comm_size(split_comm, &split_comm_size);
58
59 { // parallel interpolation process
60 // corner and cell ids for a 7 x 7 grid
61 // 56-----57-----58-----59-----60-----61-----62-----63
62 // | | | | | | | |
63 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
64 // | | | | | | | |
65 // 48-----49-----50-----51-----52-----53-----54-----55
66 // | | | | | | | |
67 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
68 // | | | | | | | |
69 // 40-----41-----42-----43-----44-----45-----46-----47
70 // | | | | | | | |
71 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
72 // | | | | | | | |
73 // 32-----33-----34-----35-----36-----37-----38-----39
74 // | | | | | | | |
75 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
76 // | | | | | | | |
77 // 24-----25-----26-----27-----28-----29-----30-----31
78 // | | | | | | | |
79 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
80 // | | | | | | | |
81 // 16-----17-----18-----19-----20-----21-----22-----23
82 // | | | | | | | |
83 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
84 // | | | | | | | |
85 // 08-----09-----10-----11-----12-----13-----14-----15
86 // | | | | | | | |
87 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
88 // | | | | | | | |
89 // 00-----01-----02-----03-----04-----05-----06-----07
90 //
91 // the grid is distributed among the processes as follows:
92 // (index == process)
93 //
94 // case 1)
95 //
96 // 4---4---4---4---4---4---4---4
97 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
98 // 4---4---4---4---4---4---4---4
99 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
100 // 3---3---3---3---4---4---4---4
101 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
102 // 2---2---2---2---3---3---3---3
103 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
104 // 1---1---1---1---2---2---2---2
105 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
106 // 0---0---0---0---1---1---1---1
107 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
108 // 0---0---0---0---0---0---0---0
109 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
110 // 0---0---0---0---0---0---0---0
111 //
112 // case 2)
113 //
114 // 0---0---1---1---3---3---4---4
115 // | 0 | 1 | 1 | 3 | 3 | 4 | 4 |
116 // 0---0---1---1---3---3---4---4
117 // | 0 | 1 | 1 | 3 | 3 | 4 | 4 |
118 // 0---0---1---1---3---3---4---4
119 // | 0 | 1 | 1 | 3 | 3 | 4 | 4 |
120 // 0---0---1---1---3---3---4---4
121 // | 0 | 0 | 1 | 3 | 3 | 4 | 4 |
122 // 0---0---1---1---2---2---4---4
123 // | 0 | 0 | 1 | 2 | 2 | 4 | 4 |
124 // 0---0---0---1---2---2---4---4
125 // | 0 | 0 | 1 | 2 | 2 | 4 | 4 |
126 // 0---0---0---1---2---2---4---4
127 // | 0 | 0 | 1 | 2 | 2 | 4 | 4 |
128 // 0---0---0---1---2---2---4---4
129 //
130 //---------------
131 // setup
132 //---------------
133
134 enum {NUM_SRC_DECOMP = 2};
135 size_t num_points[NUM_SRC_DECOMP][NUM_REORDER_TYPES];
136 double * tgt_fields[NUM_SRC_DECOMP][NUM_REORDER_TYPES];
137 int is_tgt = split_comm_size == 1;
138 for (int case_idx = 0; case_idx < NUM_SRC_DECOMP; ++case_idx) {
139
140 double coordinates_x[NUM_SRC_DECOMP][2][10] =
141 {{{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}},
142 {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}}};
143 double coordinates_y[NUM_SRC_DECOMP][2][10] =
144 {{{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}},
145 {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}}};
146 size_t const num_cells[2][2] = {{7,7}, {9,9}};
147 size_t local_start[NUM_SRC_DECOMP][2][5][2] =
148 {{{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}},
149 {{{0,0},{1,0},{3,0},{3,3},{5,0}}, {{0,0}}}};
150 size_t local_count[NUM_SRC_DECOMP][2][5][2] =
151 {{{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{9,9}}},
152 {{{2,7},{2,7},{2,3},{2,4},{2,7}}, {{9,9}}}};
153 int with_halo = 0;
154 for (size_t i = 0; i < 10; ++i) {
155 coordinates_x[case_idx][1][i] = 2.5 + 2.0 * ((double)i / 9.0);
156 coordinates_y[case_idx][1][i] = 2.5 + 2.0 * ((double)i / 9.0);
157 }
158 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
159 coordinates_x[case_idx][is_tgt][i] *= YAC_RAD;
160 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
161 coordinates_y[case_idx][is_tgt][i] *= YAC_RAD;
162
165 coordinates_x[case_idx][is_tgt], coordinates_y[case_idx][is_tgt],
166 num_cells[is_tgt], local_start[case_idx][is_tgt][split_comm_rank],
167 local_count[case_idx][is_tgt][split_comm_rank], with_halo);
168
169 struct yac_basic_grid * grids[2] =
172
173 struct yac_dist_grid_pair * grid_pair =
174 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
175
176 struct yac_interp_field src_fields[] =
177 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
178 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
180 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
181
182 struct yac_interp_grid * interp_grid =
185
186 struct interp_method * method_stack[] =
188 yac_interp_method_fixed_new(-1.0), NULL};
189
190 struct yac_interp_weights * weights =
191 yac_interp_method_do_search(method_stack, interp_grid);
192
193 for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
194
195 struct yac_interpolation * interpolation =
197 weights, reorder_types[i], 1,
198 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
199
200 // check generated interpolation
201 {
202 double * src_field = NULL;
203 double ** src_fields = &src_field;
204 double * tgt_field = NULL;
205 double * ref_tgt_field = NULL;
206
207 if (is_tgt) {
208 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
209 ref_tgt_field =
210 xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
211 utest_compute_field_data_XYZ(
212 grid_data.vertex_coordinates, 100, ref_tgt_field);
213 for (size_t i = 0; i < 100; ++i) tgt_field[i] = -1;
214 num_points[case_idx][i] = grid_data.num_vertices;
215 tgt_fields[case_idx][i] = tgt_field;
216 } else {
217 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
218 utest_compute_field_data_XYZ(
219 grid_data.vertex_coordinates, grid_data.num_vertices, src_field);
220 }
221
222 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
223
224 if (is_tgt)
225 for (size_t i = 0; i < 100; ++i)
226 if (fabs(1.0 - ref_tgt_field[i] / tgt_field[i]) > 1e-2)
227 PUT_ERR("wrong interpolation result");
228
229 free(ref_tgt_field);
230 free(src_field);
231 }
232
233 yac_interpolation_delete(interpolation);
234 }
235
237 yac_interp_method_delete(method_stack);
238 yac_interp_grid_delete(interp_grid);
239 yac_dist_grid_pair_delete(grid_pair);
242 }
243
244 if (is_tgt) {
245 for (int case_idx = 0; case_idx < NUM_SRC_DECOMP; ++case_idx) {
246 for (int reorder_type = 0; reorder_type < NUM_REORDER_TYPES;
247 ++reorder_type) {
248
249 if (num_points[0][0] != num_points[case_idx][reorder_type])
250 PUT_ERR("error in test setup");
251 for (size_t i = 0; i < num_points[0][0]; ++i)
252 if (tgt_fields[0][0][i] != tgt_fields[case_idx][reorder_type][i])
253 PUT_ERR("interpolation results are not "
254 "decomposition-independant");
255 }
256 }
257 for (int case_idx = 0; case_idx < NUM_SRC_DECOMP; ++case_idx)
258 for (int reorder_type = 0; reorder_type < NUM_REORDER_TYPES;
259 ++reorder_type) free(tgt_fields[case_idx][reorder_type]);
260 }
261 }
262
263 { // parallel interpolation process
264 // corner and cell ids for a 7 x 7 grid
265 // 56-----57-----58-----59-----60-----61-----62-----63
266 // | | | | | | | |
267 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
268 // | | | | | | | |
269 // 48-----49-----50-----51-----52-----53-----54-----55
270 // | | | | | | | |
271 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
272 // | | | | | | | |
273 // 40-----41-----42-----43-----44-----45-----46-----47
274 // | | | | | | | |
275 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
276 // | | | | | | | |
277 // 32-----33-----34-----35-----36-----37-----38-----39
278 // | | | | | | | |
279 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
280 // | | | | | | | |
281 // 24-----25-----26-----27-----28-----29-----30-----31
282 // | | | | | | | |
283 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
284 // | | | | | | | |
285 // 16-----17-----18-----19-----20-----21-----22-----23
286 // | | | | | | | |
287 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
288 // | | | | | | | |
289 // 08-----09-----10-----11-----12-----13-----14-----15
290 // | | | | | | | |
291 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
292 // | | | | | | | |
293 // 00-----01-----02-----03-----04-----05-----06-----07
294 //
295 // the grid is distributed among the processes as follows:
296 // (index == process)
297 //
298 // 4---4---4---4---4---4---4---4
299 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
300 // 4---4---4---4---4---4---4---4
301 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
302 // 3---3---3---3---4---4---4---4
303 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
304 // 2---2---2---2---3---3---3---3
305 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
306 // 1---1---1---1---2---2---2---2
307 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
308 // 0---0---0---0---1---1---1---1
309 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
310 // 0---0---0---0---0---0---0---0
311 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
312 // 0---0---0---0---0---0---0---0
313 //
314 //---------------
315 // setup
316 //---------------
317
318 int is_tgt = split_comm_size == 1;
319 double coordinates_x[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
320 double coordinates_y[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
321 size_t const num_cells[2][2] = {{7,7}, {9,9}};
322 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
323 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{9,9}}};
324 int with_halo = 0;
325 for (size_t i = 0; i < 10; ++i) {
326 coordinates_x[1][i] = 2.0 + 3.0 * ((double)i / 9.0);
327 coordinates_y[1][i] = 2.0 + 3.0 * ((double)i / 9.0);
328 }
329 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
330 coordinates_x[is_tgt][i] *= YAC_RAD;
331 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
332 coordinates_y[is_tgt][i] *= YAC_RAD;
333
336 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
337 local_start[is_tgt][split_comm_rank],
338 local_count[is_tgt][split_comm_rank], with_halo);
339
340 struct yac_basic_grid * grids[2] =
343
344 struct yac_dist_grid_pair * grid_pair =
345 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
346
347 struct yac_interp_field src_fields[] =
348 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
349 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
351 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
352
353 struct yac_interp_grid * interp_grid =
356
357 struct interp_method * method_stack[] =
359 yac_interp_method_fixed_new(-1.0), NULL};
360
361 struct yac_interp_weights * weights =
362 yac_interp_method_do_search(method_stack, interp_grid);
363
364 for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
365
366 struct yac_interpolation * interpolation =
368 weights, reorder_types[i], 1,
369 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
370
371 // check generated interpolation
372 {
373 double * src_field = NULL;
374 double ** src_fields = &src_field;
375 double * tgt_field = NULL;
376 double * ref_tgt_field = NULL;
377
378 if (is_tgt) {
379 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
380 ref_tgt_field =
381 xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
382 utest_compute_field_data_XYZ(
383 grid_data.vertex_coordinates, 100, ref_tgt_field);
384 for (size_t i = 0; i < 100; ++i) tgt_field[i] = -1;
385 } else {
386 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
387 utest_compute_field_data_XYZ(
388 grid_data.vertex_coordinates, grid_data.num_vertices, src_field);
389 }
390
391 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
392
393 if (is_tgt)
394 for (size_t i = 0; i < 100; ++i)
395 if (fabs(1.0 - ref_tgt_field[i] / tgt_field[i]) > 1e-2)
396 PUT_ERR("wrong interpolation result");
397
398 free(ref_tgt_field);
399 free(tgt_field);
400 free(src_field);
401 }
402
403 yac_interpolation_delete(interpolation);
404 }
405
407 yac_interp_method_delete(method_stack);
408 yac_interp_grid_delete(interp_grid);
409 yac_dist_grid_pair_delete(grid_pair);
412 }
413
414 { // parallel interpolation process
415 // (source point location is at the cell centers)
416 // corner and cell ids for a 7 x 7 grid
417 // 56-----57-----58-----59-----60-----61-----62-----63
418 // | | | | | | | |
419 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
420 // | | | | | | | |
421 // 48-----49-----50-----51-----52-----53-----54-----55
422 // | | | | | | | |
423 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
424 // | | | | | | | |
425 // 40-----41-----42-----43-----44-----45-----46-----47
426 // | | | | | | | |
427 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
428 // | | | | | | | |
429 // 32-----33-----34-----35-----36-----37-----38-----39
430 // | | | | | | | |
431 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
432 // | | | | | | | |
433 // 24-----25-----26-----27-----28-----29-----30-----31
434 // | | | | | | | |
435 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
436 // | | | | | | | |
437 // 16-----17-----18-----19-----20-----21-----22-----23
438 // | | | | | | | |
439 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
440 // | | | | | | | |
441 // 08-----09-----10-----11-----12-----13-----14-----15
442 // | | | | | | | |
443 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
444 // | | | | | | | |
445 // 00-----01-----02-----03-----04-----05-----06-----07
446 //
447 // the grid is distributed among the processes as follows:
448 // (index == process)
449 //
450 // 4---4---4---4---4---4---4---4
451 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
452 // 4---4---4---4---4---4---4---4
453 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
454 // 3---3---3---3---4---4---4---4
455 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
456 // 2---2---2---2---3---3---3---3
457 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
458 // 1---1---1---1---2---2---2---2
459 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
460 // 0---0---0---0---1---1---1---1
461 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
462 // 0---0---0---0---0---0---0---0
463 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
464 // 0---0---0---0---0---0---0---0
465 //
466 //---------------
467 // setup
468 //---------------
469
470 int is_tgt = split_comm_size == 1;
471 double coordinates_x[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
472 double coordinates_y[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
473 size_t const num_cells[2][2] = {{7,7}, {9,9}};
474 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
475 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{9,9}}};
476 int with_halo = 0;
477 for (size_t i = 0; i < 10; ++i) {
478 coordinates_x[1][i] = 2.0 + 3.0 * ((double)i / 9.0);
479 coordinates_y[1][i] = 2.0 + 3.0 * ((double)i / 9.0);
480 }
481 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
482 coordinates_x[is_tgt][i] *= YAC_RAD;
483 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
484 coordinates_y[is_tgt][i] *= YAC_RAD;
485
488 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
489 local_start[is_tgt][split_comm_rank],
490 local_count[is_tgt][split_comm_rank], with_halo);
491
492 yac_coordinate_pointer src_point_coordinates = NULL;
493 if (!is_tgt) {
494 src_point_coordinates =
495 xmalloc(grid_data.num_cells * sizeof(*src_point_coordinates));
496 for (size_t i = 0; i < grid_data.num_cells; ++i) {
497 double * middle_point = src_point_coordinates[i];
498 for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
499 size_t * curr_vertices =
500 grid_data.cell_to_vertex + grid_data.cell_to_vertex_offsets[i];
501 size_t curr_num_vertices = grid_data.num_vertices_per_cell[i];
502 for (size_t j = 0; j < curr_num_vertices; ++j) {
503 double * curr_vertex_coord =
504 grid_data.vertex_coordinates[curr_vertices[j]];
505 for (size_t k = 0; k < 3; ++k)
506 middle_point[k] += curr_vertex_coord[k];
507 }
508 normalise_vector(middle_point);
509 }
510 }
511
512 struct yac_basic_grid * grids[2] =
515 if (!is_tgt)
517 grids[0], YAC_LOC_CELL, src_point_coordinates);
518
519 struct yac_dist_grid_pair * grid_pair =
520 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
521 struct yac_interp_field src_fields[] =
522 {{.location = YAC_LOC_CELL, .coordinates_idx = 0, .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 = SIZE_MAX};
526
527 struct yac_interp_grid * interp_grid =
530
531 struct interp_method * method_stack[] =
533 yac_interp_method_fixed_new(-1.0), NULL};
534
535 struct yac_interp_weights * weights =
536 yac_interp_method_do_search(method_stack, interp_grid);
537
538 for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
539
540 struct yac_interpolation * interpolation =
542 weights, reorder_types[i], 1,
543 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
544
545 // check generated interpolation
546 {
547 double * src_field = NULL;
548 double ** src_fields = &src_field;
549 double * tgt_field = NULL;
550 double * ref_tgt_field = NULL;
551
552 if (is_tgt) {
553 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
554 ref_tgt_field =
555 xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
556 utest_compute_field_data_XYZ(
557 grid_data.vertex_coordinates, 100, ref_tgt_field);
558 for (size_t i = 0; i < 100; ++i) tgt_field[i] = -1;
559 } else {
560 src_field = xmalloc(grid_data.num_cells * sizeof(*src_field));
561 utest_compute_field_data_XYZ(
562 src_point_coordinates, grid_data.num_cells, src_field);
563 }
564
565 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
566
567 if (is_tgt)
568 for (size_t i = 0; i < 100; ++i)
569 if (fabs(1.0 - ref_tgt_field[i] / tgt_field[i]) > 1e-2)
570 PUT_ERR("wrong interpolation result");
571
572 free(ref_tgt_field);
573 free(tgt_field);
574 free(src_field);
575 }
576
577 yac_interpolation_delete(interpolation);
578 }
579
581 yac_interp_method_delete(method_stack);
582 yac_interp_grid_delete(interp_grid);
583 yac_dist_grid_pair_delete(grid_pair);
586 }
587
588 { // parallel interpolation process
589 // (source point location is at the cell centers)
590 // (target grid contains cells that are not covered by any source cell)
591 // corner and cell ids for a 7 x 7 grid
592 // 56-----57-----58-----59-----60-----61-----62-----63
593 // | | | | | | | |
594 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
595 // | | | | | | | |
596 // 48-----49-----50-----51-----52-----53-----54-----55
597 // | | | | | | | |
598 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
599 // | | | | | | | |
600 // 40-----41-----42-----43-----44-----45-----46-----47
601 // | | | | | | | |
602 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
603 // | | | | | | | |
604 // 32-----33-----34-----35-----36-----37-----38-----39
605 // | | | | | | | |
606 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
607 // | | | | | | | |
608 // 24-----25-----26-----27-----28-----29-----30-----31
609 // | | | | | | | |
610 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
611 // | | | | | | | |
612 // 16-----17-----18-----19-----20-----21-----22-----23
613 // | | | | | | | |
614 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
615 // | | | | | | | |
616 // 08-----09-----10-----11-----12-----13-----14-----15
617 // | | | | | | | |
618 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
619 // | | | | | | | |
620 // 00-----01-----02-----03-----04-----05-----06-----07
621 //
622 // the grid is distributed among the processes as follows:
623 // (index == process)
624 //
625 // 4---4---4---4---4---4---4---4
626 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
627 // 4---4---4---4---4---4---4---4
628 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
629 // 3---3---3---3---4---4---4---4
630 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
631 // 2---2---2---2---3---3---3---3
632 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
633 // 1---1---1---1---2---2---2---2
634 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
635 // 0---0---0---0---1---1---1---1
636 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
637 // 0---0---0---0---0---0---0---0
638 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
639 // 0---0---0---0---0---0---0---0
640 //
641 //---------------
642 // setup
643 //---------------
644
645 int is_tgt = split_comm_size == 1;
646 double coordinates_x[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
647 double coordinates_y[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
648 size_t const num_cells[2][2] = {{7,7}, {9,9}};
649 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
650 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{9,9}}};
651 int with_halo = 0;
652 for (size_t i = 0; i < 10; ++i) {
653 coordinates_x[1][i] = 4.0 + 4.0 * ((double)i / 9.0);
654 coordinates_y[1][i] = 4.0 + 4.0 * ((double)i / 9.0);
655 }
656 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
657 coordinates_x[is_tgt][i] *= YAC_RAD;
658 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
659 coordinates_y[is_tgt][i] *= YAC_RAD;
660
663 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
664 local_start[is_tgt][split_comm_rank],
665 local_count[is_tgt][split_comm_rank], with_halo);
666
667 yac_coordinate_pointer src_point_coordinates = NULL;
668 if (!is_tgt) {
669 src_point_coordinates =
670 xmalloc(grid_data.num_cells * sizeof(*src_point_coordinates));
671 for (size_t i = 0; i < grid_data.num_cells; ++i) {
672 double * middle_point = src_point_coordinates[i];
673 for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
674 size_t * curr_vertices =
675 grid_data.cell_to_vertex + grid_data.cell_to_vertex_offsets[i];
676 size_t curr_num_vertices = grid_data.num_vertices_per_cell[i];
677 for (size_t j = 0; j < curr_num_vertices; ++j) {
678 double * curr_vertex_coord =
679 grid_data.vertex_coordinates[curr_vertices[j]];
680 for (size_t k = 0; k < 3; ++k)
681 middle_point[k] += curr_vertex_coord[k];
682 }
683 normalise_vector(middle_point);
684 }
685 }
686
687 struct yac_basic_grid * grids[2] =
690 if (!is_tgt)
692 grids[0], YAC_LOC_CELL, src_point_coordinates);
693
694 struct yac_dist_grid_pair * grid_pair =
695 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
696 struct yac_interp_field src_fields[] =
697 {{.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX}};
698 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
700 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
701
702 struct yac_interp_grid * interp_grid =
705
706 struct interp_method * method_stack[] =
708 yac_interp_method_fixed_new(-1.0), NULL};
709
710 struct yac_interp_weights * weights =
711 yac_interp_method_do_search(method_stack, interp_grid);
712
713 for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
714
715 struct yac_interpolation * interpolation =
717 weights, reorder_types[i], 1,
718 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
719
720 // check generated interpolation
721 {
722 double * src_field = NULL;
723 double ** src_fields = &src_field;
724 double * tgt_field = NULL;
725 double * ref_tgt_field = NULL;
726
727 if (is_tgt) {
728 int is_outside_src_grid[] = {
729 0,0,0,0,0,0,1,1,1,1,
730 0,0,0,0,0,0,1,1,1,1,
731 0,0,0,0,0,0,1,1,1,1,
732 0,0,0,0,0,0,1,1,1,1,
733 0,0,0,0,0,0,1,1,1,1,
734 0,0,0,0,0,0,1,1,1,1,
735 1,1,1,1,1,1,1,1,1,1,
736 1,1,1,1,1,1,1,1,1,1,
737 1,1,1,1,1,1,1,1,1,1,
738 1,1,1,1,1,1,1,1,1,1};
739 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
740 ref_tgt_field =
741 xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
742 utest_compute_field_data_XYZ(
743 grid_data.vertex_coordinates, 100, ref_tgt_field);
744 for (size_t i = 0; i < grid_data.num_vertices; ++i)
745 if (is_outside_src_grid[grid_data.vertex_ids[i]])
746 ref_tgt_field[i] = -1.0;
747 for (size_t i = 0; i < 100; ++i) tgt_field[i] = -2.0;
748 } else {
749 src_field = xmalloc(grid_data.num_cells * sizeof(*src_field));
750 utest_compute_field_data_XYZ(
751 src_point_coordinates, grid_data.num_cells, src_field);
752 }
753
754 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
755
756 if (is_tgt)
757 for (size_t i = 0; i < 100; ++i)
758 if (fabs(1.0 - ref_tgt_field[i] / tgt_field[i]) > 1e-2)
759 PUT_ERR("wrong interpolation result");
760
761 free(ref_tgt_field);
762 free(tgt_field);
763 free(src_field);
764 }
765
766 yac_interpolation_delete(interpolation);
767 }
768
770 yac_interp_method_delete(method_stack);
771 yac_interp_grid_delete(interp_grid);
772 yac_dist_grid_pair_delete(grid_pair);
775 }
776
777 { // test in which the source grid goes up to the pole
778 // parallel interpolation process
779 // corner and cell ids for a 7 x 7 grid
780 // 56-----57-----58-----59-----60-----61-----62-----63
781 // | | | | | | | |
782 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
783 // | | | | | | | |
784 // 48-----49-----50-----51-----52-----53-----54-----55
785 // | | | | | | | |
786 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
787 // | | | | | | | |
788 // 40-----41-----42-----43-----44-----45-----46-----47
789 // | | | | | | | |
790 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
791 // | | | | | | | |
792 // 32-----33-----34-----35-----36-----37-----38-----39
793 // | | | | | | | |
794 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
795 // | | | | | | | |
796 // 24-----25-----26-----27-----28-----29-----30-----31
797 // | | | | | | | |
798 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
799 // | | | | | | | |
800 // 16-----17-----18-----19-----20-----21-----22-----23
801 // | | | | | | | |
802 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
803 // | | | | | | | |
804 // 08-----09-----10-----11-----12-----13-----14-----15
805 // | | | | | | | |
806 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
807 // | | | | | | | |
808 // 00-----01-----02-----03-----04-----05-----06-----07
809 //
810 // the grid is distributed among the processes as follows:
811 // (index == process)
812 //
813 // 4---4---4---4---4---4---4---4
814 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
815 // 4---4---4---4---4---4---4---4
816 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
817 // 3---3---3---3---4---4---4---4
818 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
819 // 2---2---2---2---3---3---3---3
820 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
821 // 1---1---1---1---2---2---2---2
822 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
823 // 0---0---0---0---1---1---1---1
824 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
825 // 0---0---0---0---0---0---0---0
826 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
827 // 0---0---0---0---0---0---0---0
828 //
829 //---------------
830 // setup
831 //---------------
832
833 int is_tgt = split_comm_size == 1;
834 double coordinates_x[2][10] =
835 {{72.00,72.25,72.50,72.75,73.00,73.25,73.5,73.75},{-1.0}};
836 double coordinates_y[2][10] =
837 {{90.0,89.75,89.50,89.25,89.00,88.75,88.50,88.25},{-1.0}};
838 size_t const num_cells[2][2] = {{7,7}, {9,9}};
839 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
840 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{9,9}}};
841 int with_halo = 0;
842 for (size_t i = 0; i < 10; ++i) {
843 coordinates_x[1][i] = 72.50 + 0.90 * ((double)i / 9.0);
844 coordinates_y[1][i] = 89.06 + 0.90 * ((double)i / 9.0);
845 }
846 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
847 coordinates_x[is_tgt][i] *= YAC_RAD;
848 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
849 coordinates_y[is_tgt][i] *= YAC_RAD;
850
853 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
854 local_start[is_tgt][split_comm_rank],
855 local_count[is_tgt][split_comm_rank], with_halo);
856
857 struct yac_basic_grid * grids[2] =
860
861 if (!is_tgt) {
862
863 yac_coordinate_pointer vertex_coordinates =
864 xmalloc(grid_data.num_vertices * sizeof(*vertex_coordinates));
865 memcpy(
866 vertex_coordinates, grid_data.vertex_coordinates,
867 grid_data.num_vertices * sizeof(*vertex_coordinates));
868 // set x and y coordinate of all pole points to 0.0 in order to
869 // reproduce an error that occured in ICON
870 for (size_t i = 0; i < grid_data.num_vertices; ++i) {
871 if (fabs(1.0 - fabs(vertex_coordinates[i][2])) < 1e-9) {
872 vertex_coordinates[i][0] = 0.0;
873 vertex_coordinates[i][1] = 0.0;
874 }
875 }
877 grids[0], YAC_LOC_CORNER, vertex_coordinates);
878 }
879
880 struct yac_dist_grid_pair * grid_pair =
881 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
882
883 struct yac_interp_field src_fields[] =
884 {{.location = YAC_LOC_CORNER, .coordinates_idx = 0, .masks_idx = SIZE_MAX}};
885 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
887 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
888
889 struct yac_interp_grid * interp_grid =
892
893 struct interp_method * method_stack[] =
895 yac_interp_method_fixed_new(-1.0), NULL};
896
897 struct yac_interp_weights * weights =
898 yac_interp_method_do_search(method_stack, interp_grid);
899
900 for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
901
902 struct yac_interpolation * interpolation =
904 weights, reorder_types[i], 1,
905 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
906
907 // check generated interpolation
908 {
909 double * src_field = NULL;
910 double ** src_fields = &src_field;
911 double * tgt_field = NULL;
912 double * ref_tgt_field = NULL;
913
914 if (is_tgt) {
915 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
916 ref_tgt_field =
917 xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
918 utest_compute_field_data_XYZ(
919 grid_data.vertex_coordinates, 100, ref_tgt_field);
920 for (size_t i = 0; i < 100; ++i) tgt_field[i] = -1;
921 } else {
922 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
923 utest_compute_field_data_XYZ(
924 grid_data.vertex_coordinates, grid_data.num_vertices, src_field);
925 }
926
927 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
928
929 if (is_tgt)
930 for (size_t i = 0; i < 100; ++i)
931 if (fabs(1.0 - ref_tgt_field[i] / tgt_field[i]) > 1e-2)
932 PUT_ERR("wrong interpolation result");
933
934 free(ref_tgt_field);
935 free(tgt_field);
936 free(src_field);
937 }
938
939 yac_interpolation_delete(interpolation);
940 }
941
943 yac_interp_method_delete(method_stack);
944 yac_interp_grid_delete(interp_grid);
945 yac_dist_grid_pair_delete(grid_pair);
948 }
949
950 MPI_Comm_free(&split_comm);
951
952 xt_finalize();
953
954 MPI_Finalize();
955
956 return TEST_EXIT_CODE;
957}
958
959#if defined __INTEL_COMPILER
960#pragma intel optimization_level 0
961#elif defined _CRAYC
962#pragma _CRI noopt
963#endif
964static void utest_compute_field_data_XYZ(
965 double (*vertices)[3], size_t num_vertices, double * f_v) {
966
967 for (size_t i = 0; i < num_vertices; ++i) {
968
969 double x = vertices[i][0];
970 double y = vertices[i][1];
971 double z = vertices[i][2];
972
973 // f = 1 + x^8 + e^(2y^3) + e^(2x^2) + 10xyz
974 f_v[i] = 1.0 + pow(x, 8.0);
975 f_v[i] += exp(2.0*y*y*y);
976 f_v[i] += exp(2.0*x*x) + 10.0*x*y*z;
977
978 // f = x^3 + 3.0 * x*y^2 + 5 * x*y*z
979 // f_v[i] = x*x*x + 3.0*x*y*y + 5.0*x*y*z + 1.0;
980
981 // f_v[i] = 1.0;
982 }
983}
struct yac_basic_grid * yac_basic_grid_new(char const *name, struct yac_basic_grid_data grid_data)
Definition basic_grid.c:50
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
struct yac_basic_grid * yac_basic_grid_empty_new(char const *name)
Definition basic_grid.c:63
void yac_basic_grid_delete(struct yac_basic_grid *grid)
Definition basic_grid.c:70
void yac_dist_grid_pair_delete(struct yac_dist_grid_pair *grid_pair)
Definition dist_grid.c:2313
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
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)
#define YAC_RAD
static void normalise_vector(double v[])
Definition geometry.h:635
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)
struct yac_interp_weights * yac_interp_method_do_search(struct interp_method **method, struct yac_interp_grid *interp_grid)
struct interp_method * yac_interp_method_fixed_new(double value)
struct interp_method * yac_interp_method_hcsbb_new()
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_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
static MPI_Comm split_comm
static char const * grid_names[2]
enum yac_interp_weights_reorder_type reorder_types[]
double coordinates_x[]
size_t num_cells[2]
double coordinates_y[]
#define TEST_EXIT_CODE
Definition tests.h:14
#define PUT_ERR(string)
Definition tests.h:10
static size_t num_points
Definition yac.c:159
struct yac_basic_grid ** grids
Definition yac.c:152
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19