YAC 3.21.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"
11#include "interpolation/methods/interp_method.h"
12#include "interpolation/methods/interp_method_hcsbb.h"
13#include "interpolation/methods/interp_method_fixed.h"
14#include "dist_grid_utils.h"
15#include "yac_mpi.h"
16
17#include <mpi.h>
18#include <yaxt.h>
19
27enum {NUM_REORDER_TYPES = sizeof(reorder_types) / sizeof(reorder_types[0])};
28
29static void utest_compute_field_data_XYZ(
30 double (*vertices)[3], size_t num_vertices, double * f_v);
31static char const * grid_names[2] = {"src_grid", "tgt_grid"};
32
33int main(void) {
34
35 MPI_Init(NULL, NULL);
36
37 xt_initialize(MPI_COMM_WORLD);
38
39 int comm_rank, comm_size;
40 MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
41 MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
42 MPI_Barrier(MPI_COMM_WORLD);
43
44 if (comm_size != 6) {
45 PUT_ERR("ERROR: wrong number of processes");
46 xt_finalize();
47 MPI_Finalize();
48 return TEST_EXIT_CODE;
49 }
50
51 MPI_Comm split_comm;
52 MPI_Comm_split(MPI_COMM_WORLD, comm_rank < 1, 0, &split_comm);
53
54 int split_comm_rank, split_comm_size;
55 MPI_Comm_rank(split_comm, &split_comm_rank);
56 MPI_Comm_size(split_comm, &split_comm_size);
57
58 { // parallel interpolation process
59 // corner and cell ids for a 7 x 7 grid
60 // 56-----57-----58-----59-----60-----61-----62-----63
61 // | | | | | | | |
62 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
63 // | | | | | | | |
64 // 48-----49-----50-----51-----52-----53-----54-----55
65 // | | | | | | | |
66 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
67 // | | | | | | | |
68 // 40-----41-----42-----43-----44-----45-----46-----47
69 // | | | | | | | |
70 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
71 // | | | | | | | |
72 // 32-----33-----34-----35-----36-----37-----38-----39
73 // | | | | | | | |
74 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
75 // | | | | | | | |
76 // 24-----25-----26-----27-----28-----29-----30-----31
77 // | | | | | | | |
78 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
79 // | | | | | | | |
80 // 16-----17-----18-----19-----20-----21-----22-----23
81 // | | | | | | | |
82 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
83 // | | | | | | | |
84 // 08-----09-----10-----11-----12-----13-----14-----15
85 // | | | | | | | |
86 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
87 // | | | | | | | |
88 // 00-----01-----02-----03-----04-----05-----06-----07
89 //
90 // the grid is distributed among the processes as follows:
91 // (index == process)
92 //
93 // case 1)
94 //
95 // 4---4---4---4---4---4---4---4
96 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
97 // 4---4---4---4---4---4---4---4
98 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
99 // 3---3---3---3---4---4---4---4
100 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
101 // 2---2---2---2---3---3---3---3
102 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
103 // 1---1---1---1---2---2---2---2
104 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
105 // 0---0---0---0---1---1---1---1
106 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
107 // 0---0---0---0---0---0---0---0
108 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
109 // 0---0---0---0---0---0---0---0
110 //
111 // case 2)
112 //
113 // 0---0---1---1---3---3---4---4
114 // | 0 | 1 | 1 | 3 | 3 | 4 | 4 |
115 // 0---0---1---1---3---3---4---4
116 // | 0 | 1 | 1 | 3 | 3 | 4 | 4 |
117 // 0---0---1---1---3---3---4---4
118 // | 0 | 1 | 1 | 3 | 3 | 4 | 4 |
119 // 0---0---1---1---3---3---4---4
120 // | 0 | 0 | 1 | 3 | 3 | 4 | 4 |
121 // 0---0---1---1---2---2---4---4
122 // | 0 | 0 | 1 | 2 | 2 | 4 | 4 |
123 // 0---0---0---1---2---2---4---4
124 // | 0 | 0 | 1 | 2 | 2 | 4 | 4 |
125 // 0---0---0---1---2---2---4---4
126 // | 0 | 0 | 1 | 2 | 2 | 4 | 4 |
127 // 0---0---0---1---2---2---4---4
128 //
129 //---------------
130 // setup
131 //---------------
132
133 enum {NUM_SRC_DECOMP = 2};
134 size_t num_points[NUM_SRC_DECOMP][NUM_REORDER_TYPES];
135 double * tgt_fields[NUM_SRC_DECOMP][NUM_REORDER_TYPES];
136 int is_tgt = split_comm_size == 1;
137 for (int case_idx = 0; case_idx < NUM_SRC_DECOMP; ++case_idx) {
138
139 double coordinates_x[NUM_SRC_DECOMP][2][10] =
140 {{{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}},
141 {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}}};
142 double coordinates_y[NUM_SRC_DECOMP][2][10] =
143 {{{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}},
144 {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}}};
145 size_t const num_cells[2][2] = {{7,7}, {9,9}};
146 size_t local_start[NUM_SRC_DECOMP][2][5][2] =
147 {{{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}},
148 {{{0,0},{1,0},{3,0},{3,3},{5,0}}, {{0,0}}}};
149 size_t local_count[NUM_SRC_DECOMP][2][5][2] =
150 {{{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{9,9}}},
151 {{{2,7},{2,7},{2,3},{2,4},{2,7}}, {{9,9}}}};
152 int with_halo = 0;
153 for (size_t i = 0; i < 10; ++i) {
154 coordinates_x[case_idx][1][i] = 2.5 + 2.0 * ((double)i / 9.0);
155 coordinates_y[case_idx][1][i] = 2.5 + 2.0 * ((double)i / 9.0);
156 }
157 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
158 coordinates_x[case_idx][is_tgt][i] *= YAC_RAD;
159 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
160 coordinates_y[case_idx][is_tgt][i] *= YAC_RAD;
161
162 struct yac_basic_grid_data grid_data =
163 utest_generate_basic_grid_data_reg2d(
164 coordinates_x[case_idx][is_tgt], coordinates_y[case_idx][is_tgt],
165 num_cells[is_tgt], local_start[case_idx][is_tgt][split_comm_rank],
166 local_count[case_idx][is_tgt][split_comm_rank], with_halo);
167
168 struct yac_basic_grid * grids[2] =
169 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
170 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
171
172 struct yac_dist_grid_pair * grid_pair =
173 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
174
175 struct yac_interp_field src_fields[] =
176 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
177 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
178 struct yac_interp_field tgt_field =
179 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
180
181 struct yac_interp_grid * interp_grid =
182 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
183 num_src_fields, src_fields, tgt_field);
184
185 struct interp_method * method_stack[] =
186 {yac_interp_method_hcsbb_new(),
187 yac_interp_method_fixed_new(-1.0), NULL};
188
189 struct yac_interp_weights * weights =
190 yac_interp_method_do_search(method_stack, interp_grid);
191
192 for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
193
194 struct yac_interpolation * interpolation =
195 yac_interp_weights_get_interpolation(
196 weights, reorder_types[i], 1,
197 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
198
199 // check generated interpolation
200 {
201 double * src_field = NULL;
202 double ** src_fields = &src_field;
203 double * tgt_field = NULL;
204 double * ref_tgt_field = NULL;
205
206 if (is_tgt) {
207 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
208 ref_tgt_field =
209 xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
210 utest_compute_field_data_XYZ(
211 grid_data.vertex_coordinates, 100, ref_tgt_field);
212 for (size_t i = 0; i < 100; ++i) tgt_field[i] = -1;
213 num_points[case_idx][i] = grid_data.num_vertices;
214 tgt_fields[case_idx][i] = tgt_field;
215 } else {
216 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
217 utest_compute_field_data_XYZ(
218 grid_data.vertex_coordinates, grid_data.num_vertices, src_field);
219 }
220
221 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
222
223 if (is_tgt)
224 for (size_t i = 0; i < 100; ++i)
225 if (fabs(1.0 - ref_tgt_field[i] / tgt_field[i]) > 1e-2)
226 PUT_ERR("wrong interpolation result");
227
228 free(ref_tgt_field);
229 free(src_field);
230 }
231
232 yac_interpolation_delete(interpolation);
233 }
234
235 yac_interp_weights_delete(weights);
236 yac_interp_method_delete(method_stack);
237 yac_interp_grid_delete(interp_grid);
238 yac_dist_grid_pair_delete(grid_pair);
239 yac_basic_grid_delete(grids[1]);
240 yac_basic_grid_delete(grids[0]);
241 }
242
243 if (is_tgt) {
244 for (int case_idx = 0; case_idx < NUM_SRC_DECOMP; ++case_idx) {
245 for (int reorder_type = 0; reorder_type < NUM_REORDER_TYPES;
246 ++reorder_type) {
247
248 if (num_points[0][0] != num_points[case_idx][reorder_type])
249 PUT_ERR("error in test setup");
250 for (size_t i = 0; i < num_points[0][0]; ++i)
251 if (tgt_fields[0][0][i] != tgt_fields[case_idx][reorder_type][i])
252 PUT_ERR("interpolation results are not "
253 "decomposition-independant");
254 }
255 }
256 for (int case_idx = 0; case_idx < NUM_SRC_DECOMP; ++case_idx)
257 for (int reorder_type = 0; reorder_type < NUM_REORDER_TYPES;
258 ++reorder_type) free(tgt_fields[case_idx][reorder_type]);
259 }
260 }
261
262 { // parallel interpolation process
263 // corner and cell ids for a 7 x 7 grid
264 // 56-----57-----58-----59-----60-----61-----62-----63
265 // | | | | | | | |
266 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
267 // | | | | | | | |
268 // 48-----49-----50-----51-----52-----53-----54-----55
269 // | | | | | | | |
270 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
271 // | | | | | | | |
272 // 40-----41-----42-----43-----44-----45-----46-----47
273 // | | | | | | | |
274 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
275 // | | | | | | | |
276 // 32-----33-----34-----35-----36-----37-----38-----39
277 // | | | | | | | |
278 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
279 // | | | | | | | |
280 // 24-----25-----26-----27-----28-----29-----30-----31
281 // | | | | | | | |
282 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
283 // | | | | | | | |
284 // 16-----17-----18-----19-----20-----21-----22-----23
285 // | | | | | | | |
286 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
287 // | | | | | | | |
288 // 08-----09-----10-----11-----12-----13-----14-----15
289 // | | | | | | | |
290 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
291 // | | | | | | | |
292 // 00-----01-----02-----03-----04-----05-----06-----07
293 //
294 // the grid is distributed among the processes as follows:
295 // (index == process)
296 //
297 // 4---4---4---4---4---4---4---4
298 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
299 // 4---4---4---4---4---4---4---4
300 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
301 // 3---3---3---3---4---4---4---4
302 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
303 // 2---2---2---2---3---3---3---3
304 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
305 // 1---1---1---1---2---2---2---2
306 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
307 // 0---0---0---0---1---1---1---1
308 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
309 // 0---0---0---0---0---0---0---0
310 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
311 // 0---0---0---0---0---0---0---0
312 //
313 //---------------
314 // setup
315 //---------------
316
317 int is_tgt = split_comm_size == 1;
318 double coordinates_x[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
319 double coordinates_y[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
320 size_t const num_cells[2][2] = {{7,7}, {9,9}};
321 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
322 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{9,9}}};
323 int with_halo = 0;
324 for (size_t i = 0; i < 10; ++i) {
325 coordinates_x[1][i] = 2.0 + 3.0 * ((double)i / 9.0);
326 coordinates_y[1][i] = 2.0 + 3.0 * ((double)i / 9.0);
327 }
328 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
329 coordinates_x[is_tgt][i] *= YAC_RAD;
330 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
331 coordinates_y[is_tgt][i] *= YAC_RAD;
332
333 struct yac_basic_grid_data grid_data =
334 utest_generate_basic_grid_data_reg2d(
335 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
336 local_start[is_tgt][split_comm_rank],
337 local_count[is_tgt][split_comm_rank], with_halo);
338
339 struct yac_basic_grid * grids[2] =
340 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
341 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
342
343 struct yac_dist_grid_pair * grid_pair =
344 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
345
346 struct yac_interp_field src_fields[] =
347 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
348 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
349 struct yac_interp_field tgt_field =
350 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
351
352 struct yac_interp_grid * interp_grid =
353 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
354 num_src_fields, src_fields, tgt_field);
355
356 struct interp_method * method_stack[] =
357 {yac_interp_method_hcsbb_new(),
358 yac_interp_method_fixed_new(-1.0), NULL};
359
360 struct yac_interp_weights * weights =
361 yac_interp_method_do_search(method_stack, interp_grid);
362
363 for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
364
365 struct yac_interpolation * interpolation =
366 yac_interp_weights_get_interpolation(
367 weights, reorder_types[i], 1,
368 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
369
370 // check generated interpolation
371 {
372 double * src_field = NULL;
373 double ** src_fields = &src_field;
374 double * tgt_field = NULL;
375 double * ref_tgt_field = NULL;
376
377 if (is_tgt) {
378 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
379 ref_tgt_field =
380 xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
381 utest_compute_field_data_XYZ(
382 grid_data.vertex_coordinates, 100, ref_tgt_field);
383 for (size_t i = 0; i < 100; ++i) tgt_field[i] = -1;
384 } else {
385 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
386 utest_compute_field_data_XYZ(
387 grid_data.vertex_coordinates, grid_data.num_vertices, src_field);
388 }
389
390 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
391
392 if (is_tgt)
393 for (size_t i = 0; i < 100; ++i)
394 if (fabs(1.0 - ref_tgt_field[i] / tgt_field[i]) > 1e-2)
395 PUT_ERR("wrong interpolation result");
396
397 free(ref_tgt_field);
398 free(tgt_field);
399 free(src_field);
400 }
401
402 yac_interpolation_delete(interpolation);
403 }
404
405 yac_interp_weights_delete(weights);
406 yac_interp_method_delete(method_stack);
407 yac_interp_grid_delete(interp_grid);
408 yac_dist_grid_pair_delete(grid_pair);
409 yac_basic_grid_delete(grids[1]);
410 yac_basic_grid_delete(grids[0]);
411 }
412
413 { // parallel interpolation process
414 // (source point location is at the cell centers)
415 // corner and cell ids for a 7 x 7 grid
416 // 56-----57-----58-----59-----60-----61-----62-----63
417 // | | | | | | | |
418 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
419 // | | | | | | | |
420 // 48-----49-----50-----51-----52-----53-----54-----55
421 // | | | | | | | |
422 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
423 // | | | | | | | |
424 // 40-----41-----42-----43-----44-----45-----46-----47
425 // | | | | | | | |
426 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
427 // | | | | | | | |
428 // 32-----33-----34-----35-----36-----37-----38-----39
429 // | | | | | | | |
430 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
431 // | | | | | | | |
432 // 24-----25-----26-----27-----28-----29-----30-----31
433 // | | | | | | | |
434 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
435 // | | | | | | | |
436 // 16-----17-----18-----19-----20-----21-----22-----23
437 // | | | | | | | |
438 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
439 // | | | | | | | |
440 // 08-----09-----10-----11-----12-----13-----14-----15
441 // | | | | | | | |
442 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
443 // | | | | | | | |
444 // 00-----01-----02-----03-----04-----05-----06-----07
445 //
446 // the grid is distributed among the processes as follows:
447 // (index == process)
448 //
449 // 4---4---4---4---4---4---4---4
450 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
451 // 4---4---4---4---4---4---4---4
452 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
453 // 3---3---3---3---4---4---4---4
454 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
455 // 2---2---2---2---3---3---3---3
456 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
457 // 1---1---1---1---2---2---2---2
458 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
459 // 0---0---0---0---1---1---1---1
460 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
461 // 0---0---0---0---0---0---0---0
462 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
463 // 0---0---0---0---0---0---0---0
464 //
465 //---------------
466 // setup
467 //---------------
468
469 int is_tgt = split_comm_size == 1;
470 double coordinates_x[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
471 double coordinates_y[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
472 size_t const num_cells[2][2] = {{7,7}, {9,9}};
473 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
474 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{9,9}}};
475 int with_halo = 0;
476 for (size_t i = 0; i < 10; ++i) {
477 coordinates_x[1][i] = 2.0 + 3.0 * ((double)i / 9.0);
478 coordinates_y[1][i] = 2.0 + 3.0 * ((double)i / 9.0);
479 }
480 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
481 coordinates_x[is_tgt][i] *= YAC_RAD;
482 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
483 coordinates_y[is_tgt][i] *= YAC_RAD;
484
485 struct yac_basic_grid_data grid_data =
486 utest_generate_basic_grid_data_reg2d(
487 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
488 local_start[is_tgt][split_comm_rank],
489 local_count[is_tgt][split_comm_rank], with_halo);
490
491 yac_coordinate_pointer src_point_coordinates = NULL;
492 if (!is_tgt) {
493 src_point_coordinates =
494 xmalloc(grid_data.num_cells * sizeof(*src_point_coordinates));
495 for (size_t i = 0; i < grid_data.num_cells; ++i) {
496 double * middle_point = src_point_coordinates[i];
497 for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
498 size_t * curr_vertices =
499 grid_data.cell_to_vertex + grid_data.cell_to_vertex_offsets[i];
500 size_t curr_num_vertices = grid_data.num_vertices_per_cell[i];
501 for (size_t j = 0; j < curr_num_vertices; ++j) {
502 double * curr_vertex_coord =
503 grid_data.vertex_coordinates[curr_vertices[j]];
504 for (size_t k = 0; k < 3; ++k)
505 middle_point[k] += curr_vertex_coord[k];
506 }
507 normalise_vector(middle_point);
508 }
509 }
510
511 struct yac_basic_grid * grids[2] =
512 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
513 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
514 if (!is_tgt)
515 yac_basic_grid_add_coordinates_nocpy(
516 grids[0], YAC_LOC_CELL, src_point_coordinates);
517
518 struct yac_dist_grid_pair * grid_pair =
519 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
520 struct yac_interp_field src_fields[] =
521 {{.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX}};
522 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
523 struct yac_interp_field tgt_field =
524 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
525
526 struct yac_interp_grid * interp_grid =
527 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
528 num_src_fields, src_fields, tgt_field);
529
530 struct interp_method * method_stack[] =
531 {yac_interp_method_hcsbb_new(),
532 yac_interp_method_fixed_new(-1.0), NULL};
533
534 struct yac_interp_weights * weights =
535 yac_interp_method_do_search(method_stack, interp_grid);
536
537 for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
538
539 struct yac_interpolation * interpolation =
540 yac_interp_weights_get_interpolation(
541 weights, reorder_types[i], 1,
542 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
543
544 // check generated interpolation
545 {
546 double * src_field = NULL;
547 double ** src_fields = &src_field;
548 double * tgt_field = NULL;
549 double * ref_tgt_field = NULL;
550
551 if (is_tgt) {
552 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
553 ref_tgt_field =
554 xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
555 utest_compute_field_data_XYZ(
556 grid_data.vertex_coordinates, 100, ref_tgt_field);
557 for (size_t i = 0; i < 100; ++i) tgt_field[i] = -1;
558 } else {
559 src_field = xmalloc(grid_data.num_cells * sizeof(*src_field));
560 utest_compute_field_data_XYZ(
561 src_point_coordinates, grid_data.num_cells, src_field);
562 }
563
564 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
565
566 if (is_tgt)
567 for (size_t i = 0; i < 100; ++i)
568 if (fabs(1.0 - ref_tgt_field[i] / tgt_field[i]) > 1e-2)
569 PUT_ERR("wrong interpolation result");
570
571 free(ref_tgt_field);
572 free(tgt_field);
573 free(src_field);
574 }
575
576 yac_interpolation_delete(interpolation);
577 }
578
579 yac_interp_weights_delete(weights);
580 yac_interp_method_delete(method_stack);
581 yac_interp_grid_delete(interp_grid);
582 yac_dist_grid_pair_delete(grid_pair);
583 yac_basic_grid_delete(grids[1]);
584 yac_basic_grid_delete(grids[0]);
585 }
586
587 { // parallel interpolation process
588 // (source point location is at the cell centers)
589 // (target grid contains cells that are not covered by any source cell)
590 // corner and cell ids for a 7 x 7 grid
591 // 56-----57-----58-----59-----60-----61-----62-----63
592 // | | | | | | | |
593 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
594 // | | | | | | | |
595 // 48-----49-----50-----51-----52-----53-----54-----55
596 // | | | | | | | |
597 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
598 // | | | | | | | |
599 // 40-----41-----42-----43-----44-----45-----46-----47
600 // | | | | | | | |
601 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
602 // | | | | | | | |
603 // 32-----33-----34-----35-----36-----37-----38-----39
604 // | | | | | | | |
605 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
606 // | | | | | | | |
607 // 24-----25-----26-----27-----28-----29-----30-----31
608 // | | | | | | | |
609 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
610 // | | | | | | | |
611 // 16-----17-----18-----19-----20-----21-----22-----23
612 // | | | | | | | |
613 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
614 // | | | | | | | |
615 // 08-----09-----10-----11-----12-----13-----14-----15
616 // | | | | | | | |
617 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
618 // | | | | | | | |
619 // 00-----01-----02-----03-----04-----05-----06-----07
620 //
621 // the grid is distributed among the processes as follows:
622 // (index == process)
623 //
624 // 4---4---4---4---4---4---4---4
625 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
626 // 4---4---4---4---4---4---4---4
627 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
628 // 3---3---3---3---4---4---4---4
629 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
630 // 2---2---2---2---3---3---3---3
631 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
632 // 1---1---1---1---2---2---2---2
633 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
634 // 0---0---0---0---1---1---1---1
635 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
636 // 0---0---0---0---0---0---0---0
637 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
638 // 0---0---0---0---0---0---0---0
639 //
640 //---------------
641 // setup
642 //---------------
643
644 int is_tgt = split_comm_size == 1;
645 double coordinates_x[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
646 double coordinates_y[2][10] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},{-1.0}};
647 size_t const num_cells[2][2] = {{7,7}, {9,9}};
648 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
649 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{9,9}}};
650 int with_halo = 0;
651 for (size_t i = 0; i < 10; ++i) {
652 coordinates_x[1][i] = 4.0 + 4.0 * ((double)i / 9.0);
653 coordinates_y[1][i] = 4.0 + 4.0 * ((double)i / 9.0);
654 }
655 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
656 coordinates_x[is_tgt][i] *= YAC_RAD;
657 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
658 coordinates_y[is_tgt][i] *= YAC_RAD;
659
660 struct yac_basic_grid_data grid_data =
661 utest_generate_basic_grid_data_reg2d(
662 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
663 local_start[is_tgt][split_comm_rank],
664 local_count[is_tgt][split_comm_rank], with_halo);
665
666 yac_coordinate_pointer src_point_coordinates = NULL;
667 if (!is_tgt) {
668 src_point_coordinates =
669 xmalloc(grid_data.num_cells * sizeof(*src_point_coordinates));
670 for (size_t i = 0; i < grid_data.num_cells; ++i) {
671 double * middle_point = src_point_coordinates[i];
672 for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
673 size_t * curr_vertices =
674 grid_data.cell_to_vertex + grid_data.cell_to_vertex_offsets[i];
675 size_t curr_num_vertices = grid_data.num_vertices_per_cell[i];
676 for (size_t j = 0; j < curr_num_vertices; ++j) {
677 double * curr_vertex_coord =
678 grid_data.vertex_coordinates[curr_vertices[j]];
679 for (size_t k = 0; k < 3; ++k)
680 middle_point[k] += curr_vertex_coord[k];
681 }
682 normalise_vector(middle_point);
683 }
684 }
685
686 struct yac_basic_grid * grids[2] =
687 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
688 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
689 if (!is_tgt)
690 yac_basic_grid_add_coordinates_nocpy(
691 grids[0], YAC_LOC_CELL, src_point_coordinates);
692
693 struct yac_dist_grid_pair * grid_pair =
694 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
695 struct yac_interp_field src_fields[] =
696 {{.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = SIZE_MAX}};
697 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
698 struct yac_interp_field tgt_field =
699 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
700
701 struct yac_interp_grid * interp_grid =
702 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
703 num_src_fields, src_fields, tgt_field);
704
705 struct interp_method * method_stack[] =
706 {yac_interp_method_hcsbb_new(),
707 yac_interp_method_fixed_new(-1.0), NULL};
708
709 struct yac_interp_weights * weights =
710 yac_interp_method_do_search(method_stack, interp_grid);
711
712 for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
713
714 struct yac_interpolation * interpolation =
715 yac_interp_weights_get_interpolation(
716 weights, reorder_types[i], 1,
717 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
718
719 // check generated interpolation
720 {
721 double * src_field = NULL;
722 double ** src_fields = &src_field;
723 double * tgt_field = NULL;
724 double * ref_tgt_field = NULL;
725
726 if (is_tgt) {
727 int is_outside_src_grid[] = {
728 0,0,0,0,0,0,1,1,1,1,
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 1,1,1,1,1,1,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 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
739 ref_tgt_field =
740 xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
741 utest_compute_field_data_XYZ(
742 grid_data.vertex_coordinates, 100, ref_tgt_field);
743 for (size_t i = 0; i < grid_data.num_vertices; ++i)
744 if (is_outside_src_grid[grid_data.vertex_ids[i]])
745 ref_tgt_field[i] = -1.0;
746 for (size_t i = 0; i < 100; ++i) tgt_field[i] = -2.0;
747 } else {
748 src_field = xmalloc(grid_data.num_cells * sizeof(*src_field));
749 utest_compute_field_data_XYZ(
750 src_point_coordinates, grid_data.num_cells, src_field);
751 }
752
753 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
754
755 if (is_tgt)
756 for (size_t i = 0; i < 100; ++i)
757 if (fabs(1.0 - ref_tgt_field[i] / tgt_field[i]) > 1e-2)
758 PUT_ERR("wrong interpolation result");
759
760 free(ref_tgt_field);
761 free(tgt_field);
762 free(src_field);
763 }
764
765 yac_interpolation_delete(interpolation);
766 }
767
768 yac_interp_weights_delete(weights);
769 yac_interp_method_delete(method_stack);
770 yac_interp_grid_delete(interp_grid);
771 yac_dist_grid_pair_delete(grid_pair);
772 yac_basic_grid_delete(grids[1]);
773 yac_basic_grid_delete(grids[0]);
774 }
775
776 { // test in which the source grid goes up to the pole
777 // parallel interpolation process
778 // corner and cell ids for a 7 x 7 grid
779 // 56-----57-----58-----59-----60-----61-----62-----63
780 // | | | | | | | |
781 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
782 // | | | | | | | |
783 // 48-----49-----50-----51-----52-----53-----54-----55
784 // | | | | | | | |
785 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
786 // | | | | | | | |
787 // 40-----41-----42-----43-----44-----45-----46-----47
788 // | | | | | | | |
789 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
790 // | | | | | | | |
791 // 32-----33-----34-----35-----36-----37-----38-----39
792 // | | | | | | | |
793 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
794 // | | | | | | | |
795 // 24-----25-----26-----27-----28-----29-----30-----31
796 // | | | | | | | |
797 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
798 // | | | | | | | |
799 // 16-----17-----18-----19-----20-----21-----22-----23
800 // | | | | | | | |
801 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
802 // | | | | | | | |
803 // 08-----09-----10-----11-----12-----13-----14-----15
804 // | | | | | | | |
805 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
806 // | | | | | | | |
807 // 00-----01-----02-----03-----04-----05-----06-----07
808 //
809 // the grid is distributed among the processes as follows:
810 // (index == process)
811 //
812 // 4---4---4---4---4---4---4---4
813 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
814 // 4---4---4---4---4---4---4---4
815 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
816 // 3---3---3---3---4---4---4---4
817 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
818 // 2---2---2---2---3---3---3---3
819 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
820 // 1---1---1---1---2---2---2---2
821 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
822 // 0---0---0---0---1---1---1---1
823 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
824 // 0---0---0---0---0---0---0---0
825 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
826 // 0---0---0---0---0---0---0---0
827 //
828 //---------------
829 // setup
830 //---------------
831
832 int is_tgt = split_comm_size == 1;
833 double coordinates_x[2][10] =
834 {{72.00,72.25,72.50,72.75,73.00,73.25,73.5,73.75},{-1.0}};
835 double coordinates_y[2][10] =
836 {{90.0,89.75,89.50,89.25,89.00,88.75,88.50,88.25},{-1.0}};
837 size_t const num_cells[2][2] = {{7,7}, {9,9}};
838 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
839 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{9,9}}};
840 int with_halo = 0;
841 for (size_t i = 0; i < 10; ++i) {
842 coordinates_x[1][i] = 72.50 + 0.90 * ((double)i / 9.0);
843 coordinates_y[1][i] = 89.06 + 0.90 * ((double)i / 9.0);
844 }
845 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
846 coordinates_x[is_tgt][i] *= YAC_RAD;
847 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
848 coordinates_y[is_tgt][i] *= YAC_RAD;
849
850 struct yac_basic_grid_data grid_data =
851 utest_generate_basic_grid_data_reg2d(
852 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
853 local_start[is_tgt][split_comm_rank],
854 local_count[is_tgt][split_comm_rank], with_halo);
855
856 struct yac_basic_grid * grids[2] =
857 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
858 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
859
860 if (!is_tgt) {
861
862 yac_coordinate_pointer vertex_coordinates =
863 xmalloc(grid_data.num_vertices * sizeof(*vertex_coordinates));
864 memcpy(
865 vertex_coordinates, grid_data.vertex_coordinates,
866 grid_data.num_vertices * sizeof(*vertex_coordinates));
867 // set x and y coordinate of all pole points to 0.0 in order to
868 // reproduce an error that occured in ICON
869 for (size_t i = 0; i < grid_data.num_vertices; ++i) {
870 if (fabs(1.0 - fabs(vertex_coordinates[i][2])) < 1e-9) {
871 vertex_coordinates[i][0] = 0.0;
872 vertex_coordinates[i][1] = 0.0;
873 }
874 }
875 yac_basic_grid_add_coordinates_nocpy(
876 grids[0], YAC_LOC_CORNER, vertex_coordinates);
877 }
878
879 struct yac_dist_grid_pair * grid_pair =
880 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
881
882 struct yac_interp_field src_fields[] =
883 {{.location = YAC_LOC_CORNER, .coordinates_idx = 0, .masks_idx = SIZE_MAX}};
884 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
885 struct yac_interp_field tgt_field =
886 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
887
888 struct yac_interp_grid * interp_grid =
889 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
890 num_src_fields, src_fields, tgt_field);
891
892 struct interp_method * method_stack[] =
893 {yac_interp_method_hcsbb_new(),
894 yac_interp_method_fixed_new(-1.0), NULL};
895
896 struct yac_interp_weights * weights =
897 yac_interp_method_do_search(method_stack, interp_grid);
898
899 for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
900
901 struct yac_interpolation * interpolation =
902 yac_interp_weights_get_interpolation(
903 weights, reorder_types[i], 1,
904 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
905
906 // check generated interpolation
907 {
908 double * src_field = NULL;
909 double ** src_fields = &src_field;
910 double * tgt_field = NULL;
911 double * ref_tgt_field = NULL;
912
913 if (is_tgt) {
914 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
915 ref_tgt_field =
916 xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
917 utest_compute_field_data_XYZ(
918 grid_data.vertex_coordinates, 100, ref_tgt_field);
919 for (size_t i = 0; i < 100; ++i) tgt_field[i] = -1;
920 } else {
921 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
922 utest_compute_field_data_XYZ(
923 grid_data.vertex_coordinates, grid_data.num_vertices, src_field);
924 }
925
926 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
927
928 if (is_tgt)
929 for (size_t i = 0; i < 100; ++i)
930 if (fabs(1.0 - ref_tgt_field[i] / tgt_field[i]) > 1e-2)
931 PUT_ERR("wrong interpolation result");
932
933 free(ref_tgt_field);
934 free(tgt_field);
935 free(src_field);
936 }
937
938 yac_interpolation_delete(interpolation);
939 }
940
941 yac_interp_weights_delete(weights);
942 yac_interp_method_delete(method_stack);
943 yac_interp_grid_delete(interp_grid);
944 yac_dist_grid_pair_delete(grid_pair);
945 yac_basic_grid_delete(grids[1]);
946 yac_basic_grid_delete(grids[0]);
947 }
948
949 { // parallel interpolation process
950 //
951 // A source point that coincides with one of the gauss points of the
952 // bounding triangle around the HCSBB_D_NNN nearest source points gets the
953 // full weight of that gauss point, because all other source points have
954 // a non-zero distance to it.
955 //
956 // The source grid consists of exactly HCSBB_D_NNN points, so all target
957 // points share a single bounding triangle. The last source point sits
958 // exactly on one of the gauss points of that triangle (it is the fixed
959 // point of "move the point onto the gauss point closest to it").
960
961 enum {NUM_SRC_POINTS = 7};
962 double src_coordinates_x[NUM_SRC_POINTS] =
963 {13.0, 8.5, 8.4999999999999982, 10.6, 9.3000000000000007, 10.9,
964 11.211064122277126};
965 double src_coordinates_y[NUM_SRC_POINTS] =
966 {20.0, 22.598076211353316, 17.401923788646684, 19.399999999999999,
967 20.800000000000001, 20.899999999999999, 20.160148507481018};
968
969 // a Delaunay triangulation of the source points
970 enum {NUM_SRC_CELLS = 9};
971 int src_num_vertices_per_cell[NUM_SRC_CELLS] =
972 {3, 3, 3, 3, 3, 3, 3, 3, 3};
973 int src_cell_to_vertex[NUM_SRC_CELLS * 3] =
974 {0,3,2, 1,5,0, 6,3,0, 5,6,0, 6,5,3, 4,5,1, 5,4,3, 4,1,2, 3,4,2};
975
976 // target points inside the source cells
977 enum {NUM_TGT_POINTS = 4};
978 double tgt_coordinates_x[NUM_TGT_POINTS] =
979 {10.903688, 10.266667, 11.603688, 9.466667};
980 double tgt_coordinates_y[NUM_TGT_POINTS] =
981 {20.153383, 20.366667, 19.853383, 19.200641};
982
983 int is_tgt = split_comm_size == 1;
984
985 struct yac_basic_grid_data grid_data =
986 is_tgt?
987 yac_generate_basic_grid_data_cloud_deg(
988 NUM_TGT_POINTS, tgt_coordinates_x, tgt_coordinates_y):
989 yac_generate_basic_grid_data_unstruct_deg(
990 NUM_SRC_POINTS, NUM_SRC_CELLS, src_num_vertices_per_cell,
991 src_coordinates_x, src_coordinates_y, src_cell_to_vertex);
992
993 struct yac_basic_grid * grids[2] =
994 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
995 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
996
997 struct yac_dist_grid_pair * grid_pair =
998 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
999
1000 struct yac_interp_field src_fields[] =
1001 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX,
1002 .masks_idx = SIZE_MAX}};
1003 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1004 struct yac_interp_field tgt_field =
1005 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX,
1006 .masks_idx = SIZE_MAX};
1007
1008 struct yac_interp_grid * interp_grid =
1009 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
1010 num_src_fields, src_fields, tgt_field);
1011
1012 struct interp_method * method_stack[] =
1013 {yac_interp_method_hcsbb_new(),
1014 yac_interp_method_fixed_new(-1.0), NULL};
1015
1016 struct yac_interp_weights * weights =
1017 yac_interp_method_do_search(method_stack, interp_grid);
1018
1019 for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
1020
1021 struct yac_interpolation * interpolation =
1022 yac_interp_weights_get_interpolation(
1023 weights, reorder_types[i], 1,
1024 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1025
1026 // a constant source field has to be reproduced exactly, because the
1027 // generated weights add up to one
1028 double * src_field = NULL;
1029 double ** src_fields_ = &src_field;
1030 double * tgt_field_data = NULL;
1031
1032 if (is_tgt) {
1033 tgt_field_data =
1034 xmalloc(grid_data.num_vertices * sizeof(*tgt_field_data));
1035 for (size_t j = 0; j < grid_data.num_vertices; ++j)
1036 tgt_field_data[j] = -1.0;
1037 } else {
1038 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
1039 for (size_t j = 0; j < grid_data.num_vertices; ++j)
1040 src_field[j] = 3.0;
1041 }
1042
1043 yac_interpolation_execute(interpolation, &src_fields_, &tgt_field_data);
1044
1045 if (is_tgt)
1046 for (size_t j = 0; j < grid_data.num_vertices; ++j)
1047 if (fabs(tgt_field_data[j] - 3.0) > 1e-5)
1048 PUT_ERR("wrong interpolation result");
1049
1050 free(tgt_field_data);
1051 free(src_field);
1052
1053 yac_interpolation_delete(interpolation);
1054 }
1055
1056 yac_interp_weights_delete(weights);
1057 yac_interp_method_delete(method_stack);
1058 yac_interp_grid_delete(interp_grid);
1059 yac_dist_grid_pair_delete(grid_pair);
1060 yac_basic_grid_delete(grids[1]);
1061 yac_basic_grid_delete(grids[0]);
1062 }
1063
1064 MPI_Comm_free(&split_comm);
1065
1066 xt_finalize();
1067
1068 MPI_Finalize();
1069
1070 return TEST_EXIT_CODE;
1071}
1072
1073#if defined __INTEL_COMPILER
1074#pragma intel optimization_level 0
1075#elif defined _CRAYC
1076#pragma _CRI noopt
1077#endif
1078static void utest_compute_field_data_XYZ(
1079 double (*vertices)[3], size_t num_vertices, double * f_v) {
1080
1081 for (size_t i = 0; i < num_vertices; ++i) {
1082
1083 double x = vertices[i][0];
1084 double y = vertices[i][1];
1085 double z = vertices[i][2];
1086
1087 // f = 1 + x^8 + e^(2y^3) + e^(2x^2) + 10xyz
1088 f_v[i] = 1.0 + pow(x, 8.0);
1089 f_v[i] += exp(2.0*y*y*y);
1090 f_v[i] += exp(2.0*x*x) + 10.0*x*y*z;
1091
1092 // f = x^3 + 3.0 * x*y^2 + 5 * x*y*z
1093 // f_v[i] = x*x*x + 3.0*x*y*y + 5.0*x*y*z + 1.0;
1094
1095 // f_v[i] = 1.0;
1096 }
1097}
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:57
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:202
yac_basic_grid_empty_new
struct yac_basic_grid * yac_basic_grid_empty_new(char const *name)
Definition basic_grid.c:70
yac_basic_grid_delete
void yac_basic_grid_delete(struct yac_basic_grid *grid)
Definition basic_grid.c:77
yac_generate_basic_grid_data_unstruct_deg
struct yac_basic_grid_data yac_generate_basic_grid_data_unstruct_deg(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell, double *x_vertices, double *y_vertices, int *cell_to_vertex)
Definition grid_unstruct.c:270
yac_generate_basic_grid_data_cloud_deg
struct yac_basic_grid_data yac_generate_basic_grid_data_cloud_deg(size_t nbr_points, double *x_points, double *y_points)
Definition grid_cloud.c:58
yac_dist_grid_pair_delete
void yac_dist_grid_pair_delete(struct yac_dist_grid_pair *grid_pair)
Definition dist_grid.c:2377
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:2089
dist_grid_utils.h
YAC_RAD
#define YAC_RAD
Definition dynamic_config.c:12
normalise_vector
static void normalise_vector(double v[])
Definition geometry.h:754
yac_interp_grid_delete
void yac_interp_grid_delete(struct yac_interp_grid *interp_grid)
Definition interp_grid.c:556
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)
Delete an interpolation stack and free its resources (but not the pointer array).
Definition interp_method.c:62
yac_interp_method_do_search
struct yac_interp_weights * yac_interp_method_do_search(struct interp_method **method, struct yac_interp_grid *interp_grid)
Perform weight computation using given interpolation stack and grid.
Definition interp_method.c:16
interp_method.h
Defines the interface of the interpolation method "base class" in YAC.
yac_interp_method_fixed_new
struct interp_method * yac_interp_method_fixed_new(double value)
Definition interp_method_fixed.c:63
interp_method_fixed.h
yac_interp_method_hcsbb_new
struct interp_method * yac_interp_method_hcsbb_new()
Definition interp_method_hcsbb.c:2264
interp_method_hcsbb.h
yac_interp_weights_get_interpolation
struct yac_interpolation * yac_interp_weights_get_interpolation(struct yac_interp_weights *weights, enum yac_interp_weights_reorder_type reorder, size_t collection_size, double frac_mask_fallback_value, double scaling_factor, double scaling_summand, char const *yaxt_exchanger_name, int is_source, int is_target)
Definition interp_weights.c:4885
yac_interp_weights_delete
void yac_interp_weights_delete(struct yac_interp_weights *weights)
Definition interp_weights.c:6459
yac_interp_weights_reorder_type
yac_interp_weights_reorder_type
Definition interp_weights.h:44
YAC_MAPPING_ON_TGT
@ YAC_MAPPING_ON_TGT
weights will be applied at target processes
Definition interp_weights.h:46
YAC_MAPPING_ON_SRC
@ YAC_MAPPING_ON_SRC
weights will be applied at source processes
Definition interp_weights.h:45
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:247
yac_interpolation_delete
void yac_interpolation_delete(struct yac_interpolation *interp)
Free an interpolation object and release all resources.
Definition interpolation.c:364
YAC_FRAC_MASK_NO_VALUE
double const YAC_FRAC_MASK_NO_VALUE
Definition interpolation.c:27
YAC_LOC_CORNER
@ YAC_LOC_CORNER
Definition location.h:15
YAC_LOC_CELL
@ YAC_LOC_CELL
Definition location.h:14
test_core.weights
weights
Definition test_core.py:69
test_core.src_field
src_field
Definition test_core.py:61
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
toy_multi_curve2d.y
y
Definition toy_multi_curve2d.py:62
toy_multi_curve2d.x
x
Definition toy_multi_curve2d.py:61
xmalloc
#define xmalloc(size)
Definition ppm_xfuncs.h:66
grid_data
Definition test_dummy_coupling3_c.c:36
grid_data::num_vertices_per_cell
int * num_vertices_per_cell
Definition test_dummy_coupling3_c.c:42
grid_data::cell_to_vertex
int * cell_to_vertex
Definition test_dummy_coupling3_c.c:43
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:24
yac_dist_grid_pair
Definition dist_grid.c:128
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:122
yac_interpolation
Definition interpolation.c:30
test_common.h
split_comm
static MPI_Comm split_comm
Definition test_interp_method_file_parallel.c:55
split_comm_rank
int split_comm_rank
Definition test_interp_method_fixed_parallel.c:28
split_comm_size
int split_comm_size
Definition test_interp_method_fixed_parallel.c:28
grid_names
static char const * grid_names[2]
Definition test_interp_method_hcsbb_parallel.c:31
reorder_types
enum yac_interp_weights_reorder_type reorder_types[]
Definition test_interp_method_hcsbb_parallel.c:25
NUM_REORDER_TYPES
@ NUM_REORDER_TYPES
Definition test_interp_method_hcsbb_parallel.c:27
NUM_SRC_POINTS
@ NUM_SRC_POINTS
Definition test_interp_operator_direct.c:37
NUM_TGT_POINTS
@ NUM_TGT_POINTS
Definition test_interp_operator_direct.c:38
tgt_field_data
static double tgt_field_data[MAX_COLLECTION_SIZE][NUM_TGT_POINTS]
Definition test_interp_operator_fixed.c:31
coordinates_x
double coordinates_x[]
Definition test_interpolation_parallel3.c:105
num_cells
size_t num_cells[2]
Definition test_interpolation_parallel3.c:104
coordinates_y
double coordinates_y[]
Definition test_interpolation_parallel3.c:106
with_halo
int with_halo
Definition test_interpolation_parallel3.c:110
TEST_EXIT_CODE
#define TEST_EXIT_CODE
Definition tests.h:15
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:21