YAC 3.21.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_interp_method_creep_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
7#include "tests.h"
8#include "test_common.h"
9#include "geometry.h"
10#include "interpolation/methods/interp_method.h"
11#include "interpolation/methods/interp_method_avg.h"
12#include "interpolation/methods/interp_method_conserv.h"
13#include "interpolation/methods/interp_method_creep.h"
14#include "interpolation/methods/interp_method_fixed.h"
15#include "dist_grid_utils.h"
16#include "yac_mpi.h"
17
18#include <mpi.h>
19#include <yaxt.h>
20
26static char const * grid_names[2] = {"src_grid", "tgt_grid"};
27
28int main(void) {
29
30 MPI_Init(NULL, NULL);
31
32 xt_initialize(MPI_COMM_WORLD);
33
34 int comm_rank, comm_size;
35 MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
36 MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
37 MPI_Barrier(MPI_COMM_WORLD);
38
39 if (comm_size != 3) {
40 PUT_ERR("ERROR: wrong number of processes");
41 xt_finalize();
42 MPI_Finalize();
43 return TEST_EXIT_CODE;
44 }
45
46 MPI_Comm split_comm;
47 MPI_Comm_split(MPI_COMM_WORLD, comm_rank < 2, 0, &split_comm);
48
49 int split_comm_rank, split_comm_size;
50 MPI_Comm_rank(split_comm, &split_comm_rank);
51 MPI_Comm_size(split_comm, &split_comm_size);
52
53 {// stack that starts with linear point interpolation, followed by creep
54 // fill and fixed as backup with two source processes and a single target
55 // process
56
57 // the global source grid is a 3x2 grid:
58 // 08--14--09--15--10--16--11
59 // | | | |
60 // 08 03 10 04 12 05 13
61 // | | | |
62 // 04--07--05--09--06--11--07
63 // | | | |
64 // 01 00 03 01 05 02 06
65 // | | | |
66 // 00--00--01--02--02--04--03
67 //
68 //---------------
69 // setup
70 //---------------
71
72 int is_tgt = split_comm_size == 1;
73 double coordinates_x[] = {0.0,1.0,2.0,3.0};
74 double coordinates_y[] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0};
75 double tgt_cell_coordinates_x[] = {0.5,1.5,2.5};
76 double tgt_cell_coordinates_y[] = {0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5};
77 double tgt_cell_coordinates[24][3];
78 int tgt_cell_mask[24] = {1, 1, 1,
79 1, 1, 1,
80 1, 1, 1,
81 1, 0, 1,
82 1, 0, 0,
83 1, 1, 1,
84 0, 0, 0,
85 1, 1, 1};
86 size_t const num_cells[2][2] = {{3,2}, {3,8}};
87 size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
88 size_t local_count[2][2][2] = {{{2,2},{2,2}}, {{3,8}}};
89 int with_halo = 0;
90 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
91 coordinates_x[i] *= YAC_RAD;
92 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
93 coordinates_y[i] *= YAC_RAD;
94
95 struct yac_basic_grid_data grid_data =
96 utest_generate_basic_grid_data_reg2d(
97 coordinates_x, coordinates_y, num_cells[is_tgt],
98 local_start[is_tgt][split_comm_rank],
99 local_count[is_tgt][split_comm_rank], with_halo);
100
101 struct yac_basic_grid * grids[2] =
102 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
103 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
104
105 if (is_tgt) {
106 for (size_t i = 0, k = 0; i < num_cells[is_tgt][1]; ++i)
107 for (size_t j = 0; j < num_cells[is_tgt][0]; ++j, ++k)
108 LLtoXYZ_deg(
109 tgt_cell_coordinates_x[j], tgt_cell_coordinates_y[i],
110 tgt_cell_coordinates[k]);
111 yac_basic_grid_add_coordinates(
112 grids[0], YAC_LOC_CELL, tgt_cell_coordinates, grid_data.num_cells);
113 yac_basic_grid_add_mask(
114 grids[0], YAC_LOC_CELL, tgt_cell_mask, grid_data.num_cells, NULL);
115 }
116
117 struct yac_dist_grid_pair * grid_pair =
118 yac_dist_grid_pair_new(grids[is_tgt], grids[is_tgt^1], MPI_COMM_WORLD);
119 struct yac_interp_field src_fields[] =
120 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
121 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
122 struct yac_interp_field tgt_field =
123 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
124
125 struct yac_interp_grid * interp_grid =
126 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
127 num_src_fields, src_fields, tgt_field);
128
129 for (int creep_distance = -1; creep_distance <= 7; ++creep_distance) {
130
131 struct interp_method * method_stack[] =
132 {yac_interp_method_avg_new(YAC_INTERP_AVG_ARITHMETIC, 0),
133 yac_interp_method_creep_new(creep_distance),
134 yac_interp_method_fixed_new(-1.0), NULL};
135
136 struct yac_interp_weights * weights =
137 yac_interp_method_do_search(method_stack, interp_grid);
138
139 enum yac_interp_weights_reorder_type reorder_type[2] =
140 {YAC_MAPPING_ON_SRC, YAC_MAPPING_ON_TGT};
141
142 for (size_t i = 0; i < 2; ++i) {
143
144 struct yac_interpolation * interpolation =
145 yac_interp_weights_get_interpolation(
146 weights, reorder_type[i], 1,
147 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
148
149 // check generated interpolation
150 {
151 double * src_field = NULL;
152 double ** src_fields = &src_field;
153 double * tgt_field = NULL;
154 double * ref_tgt_field = NULL;
155 double ref_global_tgt_field[9][24] =
156 {{ 2.5, 3.5, 4.5, // creep_distance = -1
157 6.5, 7.5, 8.5,
158 6.5, 7.5, 8.5,
159 6.5, -2.0, 8.5,
160 6.5, -2.0, -2.0,
161 6.5, 6.5, 6.5,
162 -2.0, -2.0, -2.0,
163 -1.0, -1.0, -1.0},
164 { 2.5, 3.5, 4.5, // creep_distance = 0
165 6.5, 7.5, 8.5,
166 -1.0, -1.0, -1.0,
167 -1.0, -2.0, -1.0,
168 -1.0, -2.0, -2.0,
169 -1.0, -1.0, -1.0,
170 -2.0, -2.0, -2.0,
171 -1.0, -1.0, -1.0},
172 { 2.5, 3.5, 4.5, // creep_distance = 1
173 6.5, 7.5, 8.5,
174 6.5, 7.5, 8.5,
175 -1.0, -2.0, -1.0,
176 -1.0, -2.0, -2.0,
177 -1.0, -1.0, -1.0,
178 -2.0, -2.0, -2.0,
179 -1.0, -1.0, -1.0},
180 { 2.5, 3.5, 4.5, // creep_distance = 2
181 6.5, 7.5, 8.5,
182 6.5, 7.5, 8.5,
183 6.5, -2.0, 8.5,
184 -1.0, -2.0, -2.0,
185 -1.0, -1.0, -1.0,
186 -2.0, -2.0, -2.0,
187 -1.0, -1.0, -1.0},
188 { 2.5, 3.5, 4.5, // creep_distance = 3
189 6.5, 7.5, 8.5,
190 6.5, 7.5, 8.5,
191 6.5, -2.0, 8.5,
192 6.5, -2.0, -2.0,
193 -1.0, -1.0, -1.0,
194 -2.0, -2.0, -2.0,
195 -1.0, -1.0, -1.0},
196 { 2.5, 3.5, 4.5, // creep_distance = 4
197 6.5, 7.5, 8.5,
198 6.5, 7.5, 8.5,
199 6.5, -2.0, 8.5,
200 6.5, -2.0, -2.0,
201 6.5, -1.0, -1.0,
202 -2.0, -2.0, -2.0,
203 -1.0, -1.0, -1.0},
204 { 2.5, 3.5, 4.5, // creep_distance = 5
205 6.5, 7.5, 8.5,
206 6.5, 7.5, 8.5,
207 6.5, -2.0, 8.5,
208 6.5, -2.0, -2.0,
209 6.5, 6.5, -1.0,
210 -2.0, -2.0, -2.0,
211 -1.0, -1.0, -1.0},
212 { 2.5, 3.5, 4.5, // creep_distance = 6
213 6.5, 7.5, 8.5,
214 6.5, 7.5, 8.5,
215 6.5, -2.0, 8.5,
216 6.5, -2.0, -2.0,
217 6.5, 6.5, 6.5,
218 -2.0, -2.0, -2.0,
219 -1.0, -1.0, -1.0},
220 { 2.5, 3.5, 4.5, // creep_distance = 7
221 6.5, 7.5, 8.5,
222 6.5, 7.5, 8.5,
223 6.5, -2.0, 8.5,
224 6.5, -2.0, -2.0,
225 6.5, 6.5, 6.5,
226 -2.0, -2.0, -2.0,
227 -1.0, -1.0, -1.0}};
228
229 if (is_tgt) {
230 tgt_field = xmalloc(grid_data.num_cells * sizeof(*tgt_field));
231 ref_tgt_field = xmalloc(grid_data.num_cells * sizeof(*ref_tgt_field));
232 for (size_t i = 0; i < grid_data.num_cells; ++i) {
233 tgt_field[i] = -2.0;
234 ref_tgt_field[i] =
235 ref_global_tgt_field[creep_distance+1][grid_data.cell_ids[i]];
236 }
237 } else {
238 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
239 for (size_t i = 0; i < grid_data.num_vertices; ++i)
240 src_field[i] = (double)(grid_data.vertex_ids[i]);
241 }
242
243 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
244
245 if (is_tgt)
246 for (size_t i = 0; i < grid_data.num_cells; ++i)
247 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-9)
248 PUT_ERR("wrong interpolation result");
249
250 free(src_field);
251 free(tgt_field);
252 free(ref_tgt_field);
253 }
254
255 yac_interpolation_delete(interpolation);
256 }
257
258 yac_interp_weights_delete(weights);
259 yac_interp_method_delete(method_stack);
260 }
261 yac_interp_grid_delete(interp_grid);
262 yac_dist_grid_pair_delete(grid_pair);
263 yac_basic_grid_delete(grids[1]);
264 yac_basic_grid_delete(grids[0]);
265 }
266
267 {// stack with average interpolation followed by unlimited creep
268 // from a single source cell to a 3x3 target with a masked center cell
269
270 // The interpolation stack consists of the average interpolation method
271 // followed by the creep method with an unlimited creep distance.
272 //
273 // The central cell in the 3x3 target grid is masked out. Only target cell 0
274 // is interpolated by the average method. The creep method then fills the
275 // remaining target cells by creeping from cell 0. When reaching cell 8
276 // (cell on the other size of the masked cell), creep will assigne values
277 // from two different path. Since both path originate from the same
278 // source cell, the same value will be assigned to cell 8 from both paths,
279 // resulting in a merge of the two stencils, which is the aim for this test.
280
281 // the global source grid is a 1x1 grid
282 // the global target grid is a 3x3 grid
283 //
284 //---------------
285 // setup
286 //---------------
287
288 int is_tgt = split_comm_size == 1;
289 double coordinates_x[2][4] = {{0.0,1.0,0.0,0.0},
290 {0.0,1.0,2.0,3.0}};
291 double coordinates_y[2][4] = {{0.0,1.0,0.0,0.0},
292 {0.0,1.0,2.0,3.0}};
293 double tgt_cell_coordinates_x[] = {0.5,1.5,2.5};
294 double tgt_cell_coordinates_y[] = {0.5,1.5,2.5};
295 double tgt_cell_coordinates[9][3];
296 int cell_mask[9] = {1,1,1,
297 1,0,1,
298 1,1,1};
299 size_t const num_cells[2][2] = {{1,1}, {3,3}};
300 size_t local_start[2][2][2] = {{{0,0},{1,0}}, {{0,0}}};
301 size_t local_count[2][2][2] = {{{1,1},{0,1}}, {{3,3}}};
302 int with_halo = 0;
303 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
304 coordinates_x[is_tgt][i] *= YAC_RAD;
305 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
306 coordinates_y[is_tgt][i] *= YAC_RAD;
307
308 struct yac_basic_grid_data grid_data =
309 utest_generate_basic_grid_data_reg2d(
310 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
311 local_start[is_tgt][split_comm_rank],
312 local_count[is_tgt][split_comm_rank], with_halo);
313
314 struct yac_basic_grid * grids[2] =
315 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
316 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
317
318 if (is_tgt) {
319 for (size_t j = 0, k = 0; j < num_cells[is_tgt][1]; ++j)
320 for (size_t l = 0; l < num_cells[is_tgt][0]; ++l, ++k)
321 LLtoXYZ_deg(
322 tgt_cell_coordinates_x[l], tgt_cell_coordinates_y[j],
323 tgt_cell_coordinates[k]);
324 yac_basic_grid_add_coordinates(
325 grids[0], YAC_LOC_CELL, tgt_cell_coordinates, grid_data.num_cells);
326 int temp_cell_mask[9];
327 for (size_t i = 0; i < grid_data.num_cells; ++i)
328 temp_cell_mask[i] = cell_mask[grid_data.cell_ids[i]];
329 yac_basic_grid_add_mask(
330 grids[0], YAC_LOC_CELL, temp_cell_mask, grid_data.num_cells, NULL);
331 }
332
333 struct yac_dist_grid_pair * grid_pair =
334 yac_dist_grid_pair_new(grids[is_tgt], grids[is_tgt^1], MPI_COMM_WORLD);
335 struct yac_interp_field src_fields[] =
336 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX,
337 .masks_idx = SIZE_MAX}};
338 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
339 struct yac_interp_field tgt_field =
340 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
341
342 struct yac_interp_grid * interp_grid =
343 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
344 num_src_fields, src_fields, tgt_field);
345
346 struct interp_method * method_stack[] =
347 {yac_interp_method_avg_new(YAC_INTERP_AVG_ARITHMETIC, 0),
348 yac_interp_method_creep_new(-1), NULL};
349
350 struct yac_interp_weights * weights =
351 yac_interp_method_do_search(method_stack, interp_grid);
352
353 enum yac_interp_weights_reorder_type reorder_type[2] =
354 {YAC_MAPPING_ON_SRC, YAC_MAPPING_ON_TGT};
355
356 for (size_t i = 0; i < 2; ++i) {
357
358 struct yac_interpolation * interpolation =
359 yac_interp_weights_get_interpolation(
360 weights, reorder_type[i], 1,
361 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
362
363 // check generated interpolation
364 {
365 double * src_field = NULL;
366 double ** src_fields = &src_field;
367 double * tgt_field = NULL;
368 double * ref_tgt_field = NULL;
369 double ref_global_tgt_field[9] =
370 {42.0, 42.0, 42.0,
371 42.0, -2.0, 42.0,
372 42.0, 42.0, 42.0};
373
374 if (is_tgt) {
375 tgt_field = xmalloc(grid_data.num_cells * sizeof(*tgt_field));
376 ref_tgt_field = xmalloc(grid_data.num_cells * sizeof(*ref_tgt_field));
377 for (size_t j = 0; j < grid_data.num_cells; ++j) {
378 tgt_field[j] = -2.0;
379 ref_tgt_field[j] = ref_global_tgt_field[grid_data.cell_ids[j]];
380 }
381 } else {
382 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
383 for (size_t j = 0; j < grid_data.num_vertices; ++j)
384 src_field[j] = 42.0;
385 }
386
387 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
388
389 if (is_tgt)
390 for (size_t j = 0; j < grid_data.num_cells; ++j)
391 if (fabs(ref_tgt_field[j] - tgt_field[j]) > 1e-9)
392 PUT_ERR("wrong interpolation result");
393
394 free(src_field);
395 free(tgt_field);
396 free(ref_tgt_field);
397 }
398
399 yac_interpolation_delete(interpolation);
400 }
401
402 yac_interp_weights_delete(weights);
403 yac_interp_method_delete(method_stack);
404 yac_interp_grid_delete(interp_grid);
405 yac_dist_grid_pair_delete(grid_pair);
406 yac_basic_grid_delete(grids[1]);
407 yac_basic_grid_delete(grids[0]);
408 }
409
410 {// stack that starts with linear point interpolation, followed by creep
411 // fill and fixed as backup with two source processes and a single target
412 // process
413
414 // the global source grid is a 7x7 grid
415 // the global target grid is a 6x6 grid
416 //
417 //---------------
418 // setup
419 //---------------
420
421 int is_tgt = split_comm_size == 1;
422 double coordinates_x[2][8] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},
423 {0.5,1.5,2.5,3.5,4.5,5.5,6.5}};
424 double coordinates_y[2][8] = {{0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0},
425 {0.5,1.5,2.5,3.5,4.5,5.5,6.5}};
426 int vertex_mask[2][8*8] =
427 {{1,1,1,1,1,1,1,1,
428 1,1,1,1,1,1,1,1,
429 1,1,0,0,0,0,1,1,
430 1,1,0,0,0,0,1,1,
431 1,1,0,0,0,0,1,1,
432 1,1,0,0,0,0,1,1,
433 1,1,1,1,1,1,1,1,
434 1,1,1,1,1,1,1,1},
435 {1,1,1,1,1,1,1,
436 1,1,1,1,1,1,1,
437 1,1,1,0,1,1,1,
438 1,1,0,1,0,1,1,
439 1,1,1,0,1,1,1,
440 1,1,1,1,1,1,1,
441 1,1,1,1,1,1,1}};
442 size_t const num_cells[2][2] = {{7,7}, {6,6}};
443 size_t local_start[2][2][2] = {{{0,0},{3,0}}, {{0,0}}};
444 size_t local_count[2][2][2] = {{{3,7},{4,7}}, {{6,6}}};
445 int with_halo = 0;
446 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
447 coordinates_x[is_tgt][i] *= YAC_RAD;
448 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
449 coordinates_y[is_tgt][i] *= YAC_RAD;
450
451 struct yac_basic_grid_data grid_data =
452 utest_generate_basic_grid_data_reg2d(
453 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
454 local_start[is_tgt][split_comm_rank],
455 local_count[is_tgt][split_comm_rank], with_halo);
456
457 struct yac_basic_grid * grids[2] =
458 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
459 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
460
461 int temp_vertex_mask[8*8];
462 for (size_t i = 0; i < grid_data.num_vertices; ++i)
463 temp_vertex_mask[i] = vertex_mask[is_tgt][grid_data.vertex_ids[i]];
464 yac_basic_grid_add_mask(
465 grids[0], YAC_LOC_CORNER, temp_vertex_mask, grid_data.num_vertices, NULL);
466
467 struct yac_dist_grid_pair * grid_pair =
468 yac_dist_grid_pair_new(grids[is_tgt], grids[is_tgt^1], MPI_COMM_WORLD);
469 struct yac_interp_field src_fields[] =
470 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
471 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
472 struct yac_interp_field tgt_field =
473 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0};
474
475 struct yac_interp_grid * interp_grid =
476 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
477 num_src_fields, src_fields, tgt_field);
478
479 for (int creep_distance = -1; creep_distance <= 3; ++creep_distance) {
480
481 struct interp_method * method_stack[] =
482 {yac_interp_method_avg_new(YAC_INTERP_AVG_ARITHMETIC, 0),
483 yac_interp_method_creep_new(creep_distance),
484 yac_interp_method_fixed_new(-1.0), NULL};
485
486 struct yac_interp_weights * weights =
487 yac_interp_method_do_search(method_stack, interp_grid);
488
489 enum yac_interp_weights_reorder_type reorder_type[2] =
490 {YAC_MAPPING_ON_SRC, YAC_MAPPING_ON_TGT};
491
492 for (size_t i = 0; i < 2; ++i) {
493
494 struct yac_interpolation * interpolation =
495 yac_interp_weights_get_interpolation(
496 weights, reorder_type[i], 1,
497 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
498
499 // check generated interpolation
500 {
501 double * src_field = NULL;
502 double ** src_fields = &src_field;
503 double * tgt_field = NULL;
504 double * ref_tgt_field = NULL;
505 double ref_global_tgt_field[5][7*7] =
506 {{ 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, // creep_distance = -1
507 12.5, 9.0, 6.5, 7.5, 8.5, 14.0, 18.5,
508 20.5, 20.5, 13.5, -2.0, 17.5, 26.5, 26.5,
509 28.5, 28.5, -2.0, -1.0, -2.0, 34.5, 34.5,
510 36.5, 36.5, 45.5, -2.0, 49.5, 42.5, 42.5,
511 44.5, 49.0, 54.5, 55.5, 56.5, 54.0, 50.5,
512 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5},
513 { 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, // creep_distance = 0
514 12.5, -1.0, -1.0, -1.0, -1.0, -1.0, 18.5,
515 20.5, -1.0, -1.0, -2.0, -1.0, -1.0, 26.5,
516 28.5, -1.0, -2.0, -1.0, -2.0, -1.0, 34.5,
517 36.5, -1.0, -1.0, -2.0, -1.0, -1.0, 42.5,
518 44.5, -1.0, -1.0, -1.0, -1.0, -1.0, 50.5,
519 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5},
520 { 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, // creep_distance = 1
521 12.5, 9.0, 6.5, 7.5, 8.5, 14.0, 18.5,
522 20.5, 20.5, -1.0, -2.0, -1.0, 26.5, 26.5,
523 28.5, 28.5, -2.0, -1.0, -2.0, 34.5, 34.5,
524 36.5, 36.5, -1.0, -2.0, -1.0, 42.5, 42.5,
525 44.5, 49.0, 54.5, 55.5, 56.5, 54.0, 50.5,
526 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5},
527 { 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, // creep_distance = 2
528 12.5, 9.0, 6.5, 7.5, 8.5, 14.0, 18.5,
529 20.5, 20.5, 13.5, -2.0, 17.5, 26.5, 26.5,
530 28.5, 28.5, -2.0, -1.0, -2.0, 34.5, 34.5,
531 36.5, 36.5, 45.5, -2.0, 49.5, 42.5, 42.5,
532 44.5, 49.0, 54.5, 55.5, 56.5, 54.0, 50.5,
533 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5},
534 { 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, // creep_distance = 3
535 12.5, 9.0, 6.5, 7.5, 8.5, 14.0, 18.5,
536 20.5, 20.5, 13.5, -2.0, 17.5, 26.5, 26.5,
537 28.5, 28.5, -2.0, -1.0, -2.0, 34.5, 34.5,
538 36.5, 36.5, 45.5, -2.0, 49.5, 42.5, 42.5,
539 44.5, 49.0, 54.5, 55.5, 56.5, 54.0, 50.5,
540 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5}};
541
542 if (is_tgt) {
543 tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
544 ref_tgt_field = xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
545 for (size_t i = 0; i < grid_data.num_vertices; ++i) {
546 tgt_field[i] = -2.0;
547 ref_tgt_field[i] =
548 ref_global_tgt_field[creep_distance+1][grid_data.vertex_ids[i]];
549 }
550 } else {
551 src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
552 for (size_t i = 0; i < grid_data.num_vertices; ++i)
553 src_field[i] = (double)(grid_data.vertex_ids[i]);
554 }
555
556 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
557
558 if (is_tgt)
559 for (size_t i = 0; i < grid_data.num_vertices; ++i)
560 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-9)
561 PUT_ERR("wrong interpolation result");
562
563 free(src_field);
564 free(tgt_field);
565 free(ref_tgt_field);
566 }
567
568 yac_interpolation_delete(interpolation);
569 }
570
571 yac_interp_weights_delete(weights);
572 yac_interp_method_delete(method_stack);
573 }
574 yac_interp_grid_delete(interp_grid);
575 yac_dist_grid_pair_delete(grid_pair);
576 yac_basic_grid_delete(grids[1]);
577 yac_basic_grid_delete(grids[0]);
578 }
579
580 {// stack that starts with linear point interpolation, followed by creep
581 // fill and fixed as backup with two source processes and a single target
582 // process
583
584 // the global source grid is a 5x1 grid
585 // the global target grid is a 5x3 grid
586 //
587 //---------------
588 // setup
589 //---------------
590
591 int is_tgt = split_comm_size == 1;
592 double coordinates_x[6] = {0.0,1.0,2.0,3.0,4.0,5.0};
593 double coordinates_y[4] = {0.0,1.0,2.0,3.0};
594 int cell_mask[2][5*3] =
595 {{1,0,1,0,1,},
596 {1,1,1,1,1,
597 1,1,0,1,1,
598 1,1,1,1,1}};
599 size_t const num_cells[2][2] = {{5,1}, {5,3}};
600 size_t local_start[2][2][2] = {{{0,0},{3,0}}, {{0,0}}};
601 size_t local_count[2][2][2] = {{{3,1},{2,1}}, {{5,3}}};
602 int with_halo = 0;
603 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
604 coordinates_x[i] *= YAC_RAD;
605 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
606 coordinates_y[i] *= YAC_RAD;
607
608 struct yac_basic_grid_data grid_data =
609 utest_generate_basic_grid_data_reg2d(
610 coordinates_x, coordinates_y, num_cells[is_tgt],
611 local_start[is_tgt][split_comm_rank],
612 local_count[is_tgt][split_comm_rank], with_halo);
613
614 struct yac_basic_grid * grids[2] =
615 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
616 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
617
618 int temp_cell_mask[5*3];
619 for (size_t i = 0; i < grid_data.num_cells; ++i)
620 temp_cell_mask[i] = cell_mask[is_tgt][grid_data.cell_ids[i]];
621 yac_basic_grid_add_mask(
622 grids[0], YAC_LOC_CELL, temp_cell_mask, grid_data.num_cells, NULL);
623
624 struct yac_dist_grid_pair * grid_pair =
625 yac_dist_grid_pair_new(grids[is_tgt], grids[is_tgt^1], MPI_COMM_WORLD);
626 struct yac_interp_field src_fields[] =
627 {{.location = YAC_LOC_CELL, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
628 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
629 struct yac_interp_field tgt_field =
630 {.location = YAC_LOC_CELL, .coordinates_idx = SIZE_MAX, .masks_idx = 0};
631
632 struct yac_interp_grid * interp_grid =
633 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
634 num_src_fields, src_fields, tgt_field);
635
636 for (int creep_distance = -1; creep_distance <= 5; ++creep_distance) {
637
638 struct interp_method * method_stack[] =
639 {yac_interp_method_conserv_new(1, 0, 1, YAC_INTERP_CONSERV_DESTAREA),
640 yac_interp_method_creep_new(creep_distance),
641 yac_interp_method_fixed_new(-1.0), NULL};
642
643 struct yac_interp_weights * weights =
644 yac_interp_method_do_search(method_stack, interp_grid);
645
646 enum yac_interp_weights_reorder_type reorder_type[2] =
647 {YAC_MAPPING_ON_SRC, YAC_MAPPING_ON_TGT};
648
649 for (size_t i = 0; i < 2; ++i) {
650
651 struct yac_interpolation * interpolation =
652 yac_interp_weights_get_interpolation(
653 weights, reorder_type[i], 1,
654 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
655
656 // check generated interpolation
657 {
658 double * src_field = NULL;
659 double ** src_fields = &src_field;
660 double * tgt_field = NULL;
661 double * ref_tgt_field = NULL;
662 double ref_global_tgt_field[7][7*7] =
663 {{ 0.0, 1.0 , 2.0 , 3.0 , 4.0, // creep_distance = -1
664 0.0, 0.5 , -2.0 , 3.5 , 4.0,
665 0.0, 0.25, 2.0 , 3.75, 4.0},
666 { 0.0, -1.0 , 2.0 , -1.0 , 4.0, // creep_distance = 0
667 -1.0, -1.0 , -2.0 , -1.0 , -1.0,
668 -1.0, -1.0 , -1.0 , -1.0 , -1.0},
669 { 0.0, 1.0 , 2.0 , 3.0 , 4.0, // creep_distance = 1
670 0.0, -1.0 , -2.0 , -1.0 , 4.0,
671 -1.0, -1.0 , -1.0 , -1.0 , -1.0},
672 { 0.0, 1.0 , 2.0 , 3.0 , 4.0, // creep_distance = 2
673 0.0, 0.5 , -2.0 , 3.5 , 4.0,
674 0.0, -1.0 , -1.0 , -1.0 , 4.0},
675 { 0.0, 1.0 , 2.0 , 3.0 , 4.0, // creep_distance = 3
676 0.0, 0.5 , -2.0 , 3.5 , 4.0,
677 0.0, 0.25, -1.0 , 3.75, 4.0},
678 { 0.0, 1.0 , 2.0 , 3.0 , 4.0, // creep_distance = 4
679 0.0, 0.5 , -2.0 , 3.5 , 4.0,
680 0.0, 0.25, 2.0 , 3.75, 4.0},
681 { 0.0, 1.0 , 2.0 , 3.0 , 4.0, // creep_distance = 5
682 0.0, 0.5 , -2.0 , 3.5 , 4.0,
683 0.0, 0.25, 2.0 , 3.75, 4.0}};
684
685 if (is_tgt) {
686 tgt_field = xmalloc(grid_data.num_cells* sizeof(*tgt_field));
687 ref_tgt_field = xmalloc(grid_data.num_cells* sizeof(*ref_tgt_field));
688 for (size_t i = 0; i < grid_data.num_cells; ++i) {
689 tgt_field[i] = -2.0;
690 ref_tgt_field[i] =
691 ref_global_tgt_field[creep_distance+1][grid_data.cell_ids[i]];
692 }
693 } else {
694 src_field = xmalloc(grid_data.num_cells* sizeof(*src_field));
695 for (size_t i = 0; i < grid_data.num_cells; ++i)
696 src_field[i] = (double)(grid_data.cell_ids[i]);
697 }
698
699 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
700
701 if (is_tgt)
702 for (size_t i = 0; i < grid_data.num_cells; ++i)
703 if (fabs(ref_tgt_field[i] - tgt_field[i]) > 1e-9)
704 PUT_ERR("wrong interpolation result");
705
706 free(src_field);
707 free(tgt_field);
708 free(ref_tgt_field);
709 }
710
711 yac_interpolation_delete(interpolation);
712 }
713
714 yac_interp_weights_delete(weights);
715 yac_interp_method_delete(method_stack);
716 }
717 yac_interp_grid_delete(interp_grid);
718 yac_dist_grid_pair_delete(grid_pair);
719 yac_basic_grid_delete(grids[1]);
720 yac_basic_grid_delete(grids[0]);
721 }
722
723 MPI_Comm_free(&split_comm);
724
725 xt_finalize();
726
727 MPI_Finalize();
728
729 return TEST_EXIT_CODE;
730}
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_mask
size_t yac_basic_grid_add_mask(struct yac_basic_grid *grid, enum yac_location location, int const *mask, size_t count, char const *mask_name)
Definition basic_grid.c:266
yac_basic_grid_add_coordinates
size_t yac_basic_grid_add_coordinates(struct yac_basic_grid *grid, enum yac_location location, yac_coordinate_pointer coordinates, size_t count)
Definition basic_grid.c:222
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_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
LLtoXYZ_deg
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:278
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_avg_new
struct interp_method * yac_interp_method_avg_new(enum yac_interp_avg_weight_type weight_type, int partial_coverage)
Definition interp_method_avg.c:1005
interp_method_avg.h
YAC_INTERP_AVG_ARITHMETIC
@ YAC_INTERP_AVG_ARITHMETIC
Definition interp_method_avg.h:13
yac_interp_method_conserv_new
struct interp_method * yac_interp_method_conserv_new(int order, int enforced_conserv, int partial_coverage, enum yac_interp_method_conserv_normalisation normalisation)
Definition interp_method_conserv.c:1237
interp_method_conserv.h
YAC_INTERP_CONSERV_DESTAREA
@ YAC_INTERP_CONSERV_DESTAREA
Definition interp_method_conserv.h:17
yac_interp_method_creep_new
struct interp_method * yac_interp_method_creep_new(int creep_distance)
Definition interp_method_creep.c:1240
interp_method_creep.h
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_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
yac_replay.vertex_mask
vertex_mask
Definition yac_replay.py:93
xmalloc
#define xmalloc(size)
Definition ppm_xfuncs.h:66
grid_data
Definition test_dummy_coupling3_c.c:36
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
cell_mask
int * cell_mask
Definition test_dummy_coupling_c.c:64
grid_names
static char const * grid_names[2]
Definition test_interp_method_creep_parallel.c:26
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
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