YAC 3.21.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_interpolation_parallel2.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#include <stdlib.h>
6#include <stdio.h>
7#include <math.h>
8#include <string.h>
9
10#include <mpi.h>
11#include <yaxt.h>
12
13#include "tests.h"
14#include "test_common.h"
15#include "dist_grid_utils.h"
16#include "geometry.h"
21
27#define YAC_RAD (0.01745329251994329576923690768489) // M_PI / 180
28
29static void utest_submain_1(
30 MPI_Comm comp_comm, enum yac_interp_weights_reorder_type reorder_type);
31
32static void utest_submain_2(
33 MPI_Comm comp_comm, enum yac_interp_weights_reorder_type reorder_type);
34
35int main(int argc, char *argv[]) {
36
37 if (argc != 2) {
38 PUT_ERR("wrong number of arguments\n");
39 return TEST_EXIT_CODE;
40 }
41
42 enum yac_interp_weights_reorder_type reorder_type =
43 (strcmp(argv[1], "src") == 0)?YAC_MAPPING_ON_SRC:YAC_MAPPING_ON_TGT;
44
45 if ((reorder_type != YAC_MAPPING_ON_SRC) && strcmp(argv[1], "tgt")) {
46 PUT_ERR("invalid argument (has to be either \"src\" or \"tgt\")\n");
47 return TEST_EXIT_CODE;
48 }
49
50 MPI_Init(NULL, NULL);
51
52 xt_initialize(MPI_COMM_WORLD);
53
54 int comm_rank, comm_size;
55 MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
56 MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
57 MPI_Barrier(MPI_COMM_WORLD);
58
59 if (comm_size != 4) {
60 PUT_ERR("ERROR: wrong number of processes");
61 xt_finalize();
62 MPI_Finalize();
63 return TEST_EXIT_CODE;
64 }
65
66 // split processes into source an target
67
68 int comp_flag = comm_rank < 2;
69
70 MPI_Comm split_comm;
71 MPI_Comm_split(MPI_COMM_WORLD, comp_flag, 0, &split_comm);
72
73 if (comp_flag) utest_submain_1(split_comm, reorder_type);
74 else utest_submain_2(split_comm, reorder_type);
75
76 MPI_Comm_free(&split_comm);
77 xt_finalize();
78 MPI_Finalize();
79
80 return TEST_EXIT_CODE;
81}
82
83/*
84 * The grid is distributed among 2 processes.
85 *
86 * The global grid has 4x4 cells:
87 *
88 * 20--36--21--37--22--38--23--39--24
89 * | | | | |
90 * 28 12 30 13 32 14 34 15 35
91 * | | | | |
92 * 15--27--16--29--17--31--18--33--19
93 * | | | | |
94 * 19 08 21 09 23 10 25 11 26
95 * | | | | |
96 * 10--18--11--20--12--22--13--24--14
97 * | | | | |
98 * 10 04 12 05 14 06 16 07 17
99 * | | | | |
100 * 05--09--06--11--07--13--08--15--09
101 * | | | | |
102 * 01 00 03 01 05 02 07 03 08
103 * | | | | |
104 * 00--00--01--02--02--04--03--06--04
105 *
106 * The mask looks as follows (# = masked points)
107 *
108 * +-------+-------+-------+-------+
109 * | | | | |
110 * | | | | |
111 * | | | | |
112 * #---#---#---#---#-------+-------+
113 * | | | | |
114 * # # # # # | |
115 * | | | | |
116 * #---#---#---#---#-------+-------+
117 * | | | | |
118 * | | | | |
119 * | | | | |
120 * +-------+-------+-------+-------+
121 * | | | | |
122 * | | | | |
123 * | | | | |
124 * +-------+-------+-------+-------+
125 */
126static void utest_submain_2(
127 MPI_Comm comp_comm, enum yac_interp_weights_reorder_type reorder_type) {
128
129 char const local_grid_name[] = "grid2";
130 char const remote_grid_name[] = "grid1";
131
132 int my_rank;
133 MPI_Comm_rank(comp_comm, &my_rank);
134
135 double coordinates_x[] = {0.0, 1.0, 2.0, 3.0, 4.0};
136 double coordinates_y[] = {0.0, 1.0, 2.0, 3.0, 4.0};
137 double edge_coordinates_x[] = {
138 0.25,0.0,1.25,1.0,2.25,2.0,3.25,3.0,4.0,
139 0.25,0.0,1.25,1.0,2.25,2.0,3.25,3.0,4.0,
140 0.25,0.0,1.25,1.0,2.25,2.0,3.25,3.0,4.0,
141 0.25,0.0,1.25,1.0,2.25,2.0,3.25,3.0,4.0,
142 0.25,1.25,2.25,3.25};
143 double edge_coordinates_y[] = {
144 0.0,0.25,0.0,0.25,0.0,0.25,0.0,0.25,0.25,
145 1.0,1.25,1.0,1.25,1.0,1.25,1.0,1.25,1.25,
146 2.0,2.25,2.0,2.25,2.0,2.25,2.0,2.25,2.25,
147 3.0,3.25,3.0,3.25,3.0,3.25,3.0,3.25,3.25,
148 4.0,4.0,4.0,4.0};
149 double edge_coords[40][3];
150 size_t const num_cells[2] = {4,4};
151 size_t local_start[2][2] = {{0,0},{2,0}};
152 size_t local_count[2][2] = {{2,4},{2,4}};
153 int with_halo = 1;
154 for (size_t i = 0; i <= num_cells[0]; ++i) coordinates_x[i] *= YAC_RAD;
155 for (size_t i = 0; i <= num_cells[1]; ++i) coordinates_y[i] *= YAC_RAD;
156 for (size_t i = 0; i < 40; ++i)
157 LLtoXYZ_deg(edge_coordinates_x[i], edge_coordinates_y[i], edge_coords[i]);
158
160 utest_generate_basic_grid_data_reg2d(
162 local_start[my_rank], local_count[my_rank], with_halo);
163 size_t num_vertices = grid_data.num_vertices;
164 size_t num_edges = grid_data.num_edges;
165
166 struct yac_basic_grid * local_grid =
167 yac_basic_grid_new(local_grid_name, grid_data);
168 struct yac_basic_grid * remote_grid =
169 yac_basic_grid_empty_new(remote_grid_name);
170
171 yac_int masked_corner_ids[] = {10,11,12,15,16,17};
172 size_t num_masked_corners =
173 sizeof(masked_corner_ids)/sizeof(masked_corner_ids[0]);
174
175 int * corner_mask = xmalloc(num_vertices * sizeof(*corner_mask));
176 for (size_t i = 0; i < num_vertices; ++i) {
177 corner_mask[i] = 1;
178 for (size_t j = 0; j < num_masked_corners; ++j)
179 if (grid_data.vertex_ids[i] == masked_corner_ids[j]) corner_mask[i] = 0;
180 }
181
182 yac_int masked_edge_ids[] = {18,19,20,21,23,27,29};
183 size_t num_masked_edges =
184 sizeof(masked_edge_ids)/sizeof(masked_edge_ids[0]);
185
186 int * edge_mask = xmalloc(num_edges * sizeof(*edge_mask));
187 for (size_t i = 0; i < num_edges; ++i) {
188 edge_mask[i] = 1;
189 for (size_t j = 0; j < num_masked_edges; ++j)
190 if (grid_data.edge_ids[i] == masked_edge_ids[j]) edge_mask[i] = 0;
191 }
192
193 yac_coordinate_pointer edge_field_coords =
194 xmalloc(num_edges * sizeof(*edge_field_coords));
195 for (size_t i = 0; i < num_edges; ++i)
196 memcpy(edge_field_coords[i], edge_coords[grid_data.edge_ids[i]],
197 3 * sizeof(double));
198
200 local_grid, YAC_LOC_CORNER, corner_mask, NULL);
202 local_grid, YAC_LOC_EDGE, edge_mask, NULL);
204 local_grid, YAC_LOC_EDGE, edge_field_coords);
205
206 struct yac_dist_grid_pair * grid_pair =
207 yac_dist_grid_pair_new(local_grid, remote_grid, MPI_COMM_WORLD);
208
209 { // test interpolation without using a mask
210
211 struct yac_interp_field src_fields[] =
212 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
213 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
215 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
216
217 struct yac_interp_grid * interp_grid_out =
218 yac_interp_grid_new(grid_pair, local_grid_name, remote_grid_name,
220 struct yac_interp_grid * interp_grid_in =
221 yac_interp_grid_new(grid_pair, remote_grid_name, local_grid_name,
223
224 struct interp_method * method_stack_out[] =
227 NULL};
228 struct interp_method * method_stack_in[] =
231 NULL};
232
233 struct yac_interp_weights * weights_out =
234 yac_interp_method_do_search(method_stack_out, interp_grid_out);
235 struct yac_interp_weights * weights_in =
236 yac_interp_method_do_search(method_stack_in, interp_grid_in);
237
238 yac_interp_method_delete(method_stack_in);
239 yac_interp_method_delete(method_stack_out);
240
241 struct yac_interpolation * interpolation_out =
243 weights_out, reorder_type, 1,
244 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
245 struct yac_interpolation * interpolation_in =
247 weights_in, reorder_type, 1,
248 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
249
250 //---------------------
251 // do the interpolation
252 //---------------------
253
254 // source_data dimensions [collection_idx]
255 // [pointset_idx]
256 // [local_idx]
257 double * source_data_field =
258 xmalloc(num_vertices * sizeof(*source_data_field));
259 double * source_data_pointset[1] = {source_data_field}; // num_pointset == 1
260 double ** source_data[1] = {source_data_pointset}; // collection_size == 1
261 for (size_t i = 0; i < num_vertices; ++i)
262 source_data_field[i] =
263 (grid_data.core_vertex_mask[i])?
264 ((double)(grid_data.vertex_ids[i])):(-1.0);
265 // target_data dimensions [collection_idx]
266 // [local_idx]
267 double * target_data_field =
268 xmalloc(num_vertices * sizeof(*target_data_field));
269 double * target_data[1] = {target_data_field}; // collection_size == 1
270 for (size_t i = 0; i < num_vertices; ++i) target_data_field[i] = -1.0;
271
272 yac_interpolation_execute_put(interpolation_out, source_data);
273 yac_interpolation_execute_get(interpolation_in, target_data);
274
275 //----------------------------
276 // check interpolation results
277 //----------------------------
278
279 double ref_global_target_data[] = {1338,1338,1338,1338,1338,
280 1338, 3, 4, 5, 6,
281 1338, 8, 9, 10, 11,
282 1338, 13, 14, 15, 16,
283 1338, 18, 19, 20, 21};
284
285 for (size_t i = 0; i < num_vertices; ++i) {
286 if (grid_data.core_vertex_mask[i]) {
287 if (utest_double_are_unequal(
288 target_data[0][i],
289 ref_global_target_data[grid_data.vertex_ids[i]]))
290 PUT_ERR("error in interpolated data on target side\n");
291 } else {
292 if (target_data[0][i] != -1.0)
293 PUT_ERR("error in interpolated data on target side\n");
294 }
295 }
296
297 //--------
298 // cleanup
299 //--------
300
301 free(target_data_field);
302 free(source_data_field);
303 yac_interpolation_delete(interpolation_out);
304 yac_interpolation_delete(interpolation_in);
305 yac_interp_weights_delete(weights_out);
306 yac_interp_weights_delete(weights_in);
307 yac_interp_grid_delete(interp_grid_out);
308 yac_interp_grid_delete(interp_grid_in);
309 }
310
311 { // test interpolation using a source mask
312
313 struct yac_interp_field src_fields[] =
314 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
315 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
317 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
318
319 struct yac_interp_grid * interp_grid_out =
320 yac_interp_grid_new(grid_pair, local_grid_name, remote_grid_name,
322 struct yac_interp_grid * interp_grid_in =
323 yac_interp_grid_new(grid_pair, remote_grid_name, local_grid_name,
325
326 struct interp_method * method_stack_out[] =
329 NULL};
330 struct interp_method * method_stack_in[] =
333 NULL};
334
335 struct yac_interp_weights * weights_out =
336 yac_interp_method_do_search(method_stack_out, interp_grid_out);
337 struct yac_interp_weights * weights_in =
338 yac_interp_method_do_search(method_stack_in, interp_grid_in);
339
340 yac_interp_method_delete(method_stack_in);
341 yac_interp_method_delete(method_stack_out);
342
343 struct yac_interpolation * interpolation_out =
345 weights_out, reorder_type, 1,
346 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
347 struct yac_interpolation * interpolation_in =
349 weights_in, reorder_type, 1,
350 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
351
352 //---------------------
353 // do the interpolation
354 //---------------------
355
356 // source_data dimensions [collection_idx]
357 // [pointset_idx]
358 // [local_idx]
359 double * source_data_field =
360 xmalloc(num_vertices * sizeof(*source_data_field));
361 double * source_data_pointset[1] = {source_data_field}; // num_pointset == 1
362 double ** source_data[1] = {source_data_pointset}; // collection_size == 1
363 for (size_t i = 0; i < num_vertices; ++i)
364 source_data_field[i] =
365 (grid_data.core_vertex_mask[i])?
366 ((double)(grid_data.vertex_ids[i])):(-1.0);
367 // target_data dimensions [collection_idx]
368 // [local_idx]
369 double * target_data_field =
370 xmalloc(num_vertices * sizeof(*target_data_field));
371 double * target_data[1] = {target_data_field}; // collection_size == 1
372 for (size_t i = 0; i < num_vertices; ++i) target_data_field[i] = -1.0;
373
374 yac_interpolation_execute_put(interpolation_out, source_data);
375 yac_interpolation_execute_get(interpolation_in, target_data);
376
377 //----------------------------
378 // check interpolation results
379 //----------------------------
380
381 double ref_global_target_data[] = {1338,1338,1338,1338,1338,
382 1338, 3,1338,1338,1338,
383 1338, 8,1338,1338,1338,
384 1338, 13,1338,1338,1338,
385 1338, 18, 19, 20, 21};
386
387 for (size_t i = 0; i < num_vertices; ++i) {
388 if (grid_data.core_vertex_mask[i]) {
389 if (utest_double_are_unequal(
390 target_data[0][i],
391 ref_global_target_data[grid_data.vertex_ids[i]]))
392 PUT_ERR("error in interpolated data on target side\n");
393 } else {
394 if (target_data[0][i] != -1.0)
395 PUT_ERR("error in interpolated data on target side\n");
396 }
397 }
398
399 //--------
400 // cleanup
401 //--------
402
403 free(target_data_field);
404 free(source_data_field);
405 yac_interpolation_delete(interpolation_out);
406 yac_interpolation_delete(interpolation_in);
407 yac_interp_weights_delete(weights_out);
408 yac_interp_weights_delete(weights_in);
409 yac_interp_grid_delete(interp_grid_out);
410 yac_interp_grid_delete(interp_grid_in);
411 }
412
413 { // test interpolation using a source mask and allowing partial coverage
414
415 struct yac_interp_field src_fields[] =
416 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
417 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
419 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
420
421 struct yac_interp_grid * interp_grid_out =
422 yac_interp_grid_new(grid_pair, local_grid_name, remote_grid_name,
424 struct yac_interp_grid * interp_grid_in =
425 yac_interp_grid_new(grid_pair, remote_grid_name, local_grid_name,
427
428 struct interp_method * method_stack_out[] =
431 NULL};
432 struct interp_method * method_stack_in[] =
435 NULL};
436
437 struct yac_interp_weights * weights_out =
438 yac_interp_method_do_search(method_stack_out, interp_grid_out);
439 struct yac_interp_weights * weights_in =
440 yac_interp_method_do_search(method_stack_in, interp_grid_in);
441
442 yac_interp_method_delete(method_stack_in);
443 yac_interp_method_delete(method_stack_out);
444
445 struct yac_interpolation * interpolation_out =
447 weights_out, reorder_type, 1,
448 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
449 struct yac_interpolation * interpolation_in =
451 weights_in, reorder_type, 1,
452 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
453
454 //---------------------
455 // do the interpolation
456 //---------------------
457
458 // source_data dimensions [collection_idx]
459 // [pointset_idx]
460 // [local_idx]
461 double * source_data_field =
462 xmalloc(num_vertices * sizeof(*source_data_field));
463 double * source_data_pointset[1] = {source_data_field}; // num_pointset == 1
464 double ** source_data[1] = {source_data_pointset}; // collection_size == 1
465 for (size_t i = 0; i < num_vertices; ++i)
466 source_data_field[i] =
467 (grid_data.core_vertex_mask[i])?
468 ((double)(grid_data.vertex_ids[i])):(-1.0);
469 // target_data dimensions [collection_idx]
470 // [local_idx]
471 double * target_data_field =
472 xmalloc(num_vertices * sizeof(*target_data_field));
473 double * target_data[1] = {target_data_field}; // collection_size == 1
474 for (size_t i = 0; i < num_vertices; ++i) target_data_field[i] = -1.0;
475
476 yac_interpolation_execute_put(interpolation_out, source_data);
477 yac_interpolation_execute_get(interpolation_in, target_data);
478
479 //----------------------------
480 // check interpolation results
481 //----------------------------
482
483 double ref_global_target_data[] = {1338,1338, 1338,1338,1338,
484 1338, 3, 3, 2.5, 3.5,
485 1338, 8, 8.5,1338,1338,
486 1338, 13,44.0/3.0,17.5,18.5,
487 1338, 18, 19, 20, 21};
488
489 for (size_t i = 0; i < num_vertices; ++i) {
490 if (grid_data.core_vertex_mask[i]) {
491 if (fabs(target_data[0][i] -
492 ref_global_target_data[grid_data.vertex_ids[i]]) > 1e-9)
493 PUT_ERR("error in interpolated data on target side\n");
494 } else {
495 if (target_data[0][i] != -1.0)
496 PUT_ERR("error in interpolated data on target side\n");
497 }
498 }
499
500 //--------
501 // cleanup
502 //--------
503
504 free(target_data_field);
505 free(source_data_field);
506 yac_interpolation_delete(interpolation_out);
507 yac_interpolation_delete(interpolation_in);
508 yac_interp_weights_delete(weights_out);
509 yac_interp_weights_delete(weights_in);
510 yac_interp_grid_delete(interp_grid_out);
511 yac_interp_grid_delete(interp_grid_in);
512 }
513
514 { // test interpolation using a source mask and allowing partial coverage and
515 // use a target mask
516
517 struct yac_interp_grid * interp_grid_in, * interp_grid_out;
518
519 {
520 struct yac_interp_field src_fields[] =
521 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
522 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
524 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0};
525
526 interp_grid_out =
527 yac_interp_grid_new(grid_pair, local_grid_name, remote_grid_name,
528 num_src_fields, src_fields, tgt_field);
529 }
530
531 {
532 struct yac_interp_field src_fields[] =
533 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
534 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
536 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0};
537
538 interp_grid_in =
539 yac_interp_grid_new(grid_pair, remote_grid_name, local_grid_name,
540 num_src_fields, src_fields, tgt_field);
541 }
542
543 struct interp_method * method_stack_out[] =
546 NULL};
547 struct interp_method * method_stack_in[] =
550 NULL};
551
552 struct yac_interp_weights * weights_out =
553 yac_interp_method_do_search(method_stack_out, interp_grid_out);
554 struct yac_interp_weights * weights_in =
555 yac_interp_method_do_search(method_stack_in, interp_grid_in);
556
557 yac_interp_method_delete(method_stack_in);
558 yac_interp_method_delete(method_stack_out);
559
560 struct yac_interpolation * interpolation_out =
562 weights_out, reorder_type, 1,
563 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
564 struct yac_interpolation * interpolation_in =
566 weights_in, reorder_type, 1,
567 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
568
569 //---------------------
570 // do the interpolation
571 //---------------------
572
573 // source_data dimensions [collection_idx]
574 // [pointset_idx]
575 // [local_idx]
576 double * source_data_field =
577 xmalloc(num_vertices * sizeof(*source_data_field));
578 double * source_data_pointset[1] = {source_data_field}; // num_pointset == 1
579 double ** source_data[1] = {source_data_pointset}; // collection_size == 1
580 for (size_t i = 0; i < num_vertices; ++i)
581 source_data_field[i] =
582 (grid_data.core_vertex_mask[i])?
583 ((double)(grid_data.vertex_ids[i])):(-1.0);
584 // target_data dimensions [collection_idx]
585 // [local_idx]
586 double * target_data_field =
587 xmalloc(num_vertices * sizeof(*target_data_field));
588 double * target_data[1] = {target_data_field}; // collection_size == 1
589 for (size_t i = 0; i < num_vertices; ++i) target_data_field[i] = -1.0;
590
591 yac_interpolation_execute_put(interpolation_out, source_data);
592 yac_interpolation_execute_get(interpolation_in, target_data);
593
594 //----------------------------
595 // check interpolation results
596 //----------------------------
597
598 double ref_global_target_data[] = {1338,1338,1338,1338,1338,
599 1338, 3, 3, 2.5, 3.5,
600 -1, -1, -1,1338,1338,
601 -1, -1, -1,17.5,18.5,
602 1338, 18, 19, 20, 21};
603
604 for (size_t i = 0; i < num_vertices; ++i) {
605 if ((grid_data.core_vertex_mask[i]) && corner_mask[i]) {
606 if (fabs(target_data[0][i] -
607 ref_global_target_data[grid_data.vertex_ids[i]]) > 1e-9)
608 PUT_ERR("error in interpolated data on target side\n");
609 } else {
610 if (target_data[0][i] != -1.0)
611 PUT_ERR("error in interpolated data on target side\n");
612 }
613 }
614
615 //--------
616 // cleanup
617 //--------
618
619 free(target_data_field);
620 free(source_data_field);
621 yac_interpolation_delete(interpolation_out);
622 yac_interpolation_delete(interpolation_in);
623 yac_interp_weights_delete(weights_out);
624 yac_interp_weights_delete(weights_in);
625 yac_interp_grid_delete(interp_grid_out);
626 yac_interp_grid_delete(interp_grid_in);
627 }
628
629 { // test interpolation using a source mask and allowing partial coverage and
630 // use a target mask (target points on edges)
631
632 struct yac_interp_grid * interp_grid_out, * interp_grid_in;
633 {
634 struct yac_interp_field src_fields[] =
635 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
636 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
638 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
639 interp_grid_out =
640 yac_interp_grid_new(grid_pair, local_grid_name, remote_grid_name,
641 num_src_fields, src_fields, tgt_field);
642 }
643 {
644 struct yac_interp_field src_fields[] =
645 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
646 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
648 {.location = YAC_LOC_EDGE, .coordinates_idx = 0, .masks_idx = 0};
649 interp_grid_in =
650 yac_interp_grid_new(grid_pair, remote_grid_name, local_grid_name,
651 num_src_fields, src_fields, tgt_field);
652 }
653
654 struct interp_method * method_stack_out[] =
657 NULL};
658 struct interp_method * method_stack_in[] =
661 NULL};
662
663 struct yac_interp_weights * weights_out =
664 yac_interp_method_do_search(method_stack_out, interp_grid_out);
665 struct yac_interp_weights * weights_in =
666 yac_interp_method_do_search(method_stack_in, interp_grid_in);
667
668 yac_interp_method_delete(method_stack_in);
669 yac_interp_method_delete(method_stack_out);
670
671 struct yac_interpolation * interpolation_out =
673 weights_out, reorder_type, 1,
674 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
675 struct yac_interpolation * interpolation_in =
677 weights_in, reorder_type, 1,
678 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
679
680 //---------------------
681 // do the interpolation
682 //---------------------
683
684 // source_data dimensions [collection_idx]
685 // [pointset_idx]
686 // [local_idx]
687 double * source_data_field =
688 xmalloc(num_vertices * sizeof(*source_data_field));
689 double * source_data_pointset[1] = {source_data_field}; // num_pointset == 1
690 double ** source_data[1] = {source_data_pointset}; // collection_size == 1
691 for (size_t i = 0; i < num_vertices; ++i)
692 source_data_field[i] =
693 (grid_data.core_vertex_mask[i])?
694 ((double)(grid_data.vertex_ids[i])):(-1.0);
695 // target_data dimensions [collection_idx]
696 // [local_idx]
697 double * target_data_field =
698 xmalloc(num_edges * sizeof(*target_data_field));
699 double * target_data[1] = {target_data_field}; // collection_size == 1
700 for (size_t i = 0; i < num_edges; ++i) target_data_field[i] = -1.0;
701
702 yac_interpolation_execute_put(interpolation_out, source_data);
703 yac_interpolation_execute_get(interpolation_in, target_data);
704
705 //----------------------------
706 // check interpolation results
707 //----------------------------
708
709 double ref_global_target_data[] = {
710 1338,1338,1338,1338, 1338, 1338,1338,1338,1338,
711 1338,1338, 3, 3, 3, 3, 2.5, 2.5, 3.5,
712 -1, -1, -1, -1, 8.5, -1,1338,1338,1338,
713 -1,1338, -1, 13,44.0/3.0,44.0/3.0,17.5,17.5,18.5,
714 1338, 18, 19, 20};
715
716 for (size_t i = 0; i < num_edges; ++i) {
717 if ((grid_data.core_edge_mask[i]) && edge_mask[i]) {
718 if (fabs(target_data[0][i] -
719 ref_global_target_data[grid_data.edge_ids[i]]) > 1e-9)
720 PUT_ERR("error in interpolated data on target side\n");
721 } else {
722 if (target_data[0][i] != -1.0)
723 PUT_ERR("error in interpolated data on target side\n");
724 }
725 }
726
727 //--------
728 // cleanup
729 //--------
730
731 free(target_data_field);
732 free(source_data_field);
733 yac_interpolation_delete(interpolation_out);
734 yac_interpolation_delete(interpolation_in);
735 yac_interp_weights_delete(weights_out);
736 yac_interp_weights_delete(weights_in);
737 yac_interp_grid_delete(interp_grid_out);
738 yac_interp_grid_delete(interp_grid_in);
739 }
740
741 { // test interpolation using a fractional mask
742
743 struct yac_interp_grid * interp_grid;
744 {
745 struct yac_interp_field src_fields[] =
746 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
747 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
749 {.location = YAC_LOC_EDGE, .coordinates_idx = 0, .masks_idx = 0};
750 interp_grid =
751 yac_interp_grid_new(grid_pair, remote_grid_name, local_grid_name,
752 num_src_fields, src_fields, tgt_field);
753 }
754
755 struct interp_method * method_stack[] =
758 NULL};
759
760 struct yac_interp_weights * weights =
761 yac_interp_method_do_search(method_stack, interp_grid);
762
763 yac_interp_method_delete(method_stack);
764
765 double frac_mask_value = -1337.0;
767 interpolations[0] =
769 weights, reorder_type, 1, frac_mask_value, 1.0, 0.0, NULL, 1, 1);
770 interpolations[1] =
772
773 double * target_data_field =
774 xmalloc(num_edges * sizeof(*target_data_field));
775 double * target_data[1] = {target_data_field}; // collection_size == 1
776
777 for (int interp_idx = 0; interp_idx < 2; ++interp_idx) {
778
779 //---------------------
780 // do the interpolation
781 //---------------------
782
783 for (size_t i = 0; i < num_edges; ++i) target_data_field[i] = -1.0;
784
785 yac_interpolation_execute_get(interpolations[interp_idx], target_data);
786
787 //----------------------------
788 // check interpolation results
789 //----------------------------
790
791 double ref_global_target_data[] = {
792 1337,1337,1337,1337, 1337, 1337,1337,1337,1337,
793 1337,1337, 3, 3, 3, 3, 2.5, 2.5, 3.5,
794 -1, -1, -1, -1, 8.5, -1,1337,1337,1337,
795 -1,1337, -1, 13,44.0/3.0,44.0/3.0,17.5,17.5,18.5,
796 1337, 18, 19, 20};
797
798 for (size_t i = 0; i < num_edges; ++i) {
799 if ((grid_data.core_edge_mask[i]) && edge_mask[i]) {
800 if (fabs(target_data[0][i] -
801 ref_global_target_data[grid_data.edge_ids[i]]) > 1e-9)
802 PUT_ERR("error in interpolated data on target side\n");
803 } else {
804 if (target_data[0][i] != -1.0)
805 PUT_ERR("error in interpolated data on target side\n");
806 }
807 }
808 }
809
810 //--------
811 // cleanup
812 //--------
813
814 free(target_data_field);
815 for (int interp_idx = 0; interp_idx < 2; ++interp_idx)
818 yac_interp_grid_delete(interp_grid);
819 }
820
821 //--------
822 // cleanup
823 //--------
824
825 yac_dist_grid_pair_delete(grid_pair);
826 yac_basic_grid_delete(remote_grid);
827 yac_basic_grid_delete(local_grid);
828}
829
830/*
831 * The grid is distributed among 2 processes.
832 *
833 * The global grid has 4x4 cells:
834 *
835 * 20--36--21--37--22--38--23--39--24
836 * | | | | |
837 * 28 12 30 13 32 14 34 15 35
838 * | | | | |
839 * 15--27--16--29--17--31--18--33--19
840 * | | | | |
841 * 19 08 21 09 23 10 25 11 26
842 * | | | | |
843 * 10--18--11--20--12--22--13--24--14
844 * | | | | |
845 * 10 04 12 05 14 06 16 07 17
846 * | | | | |
847 * 05--09--06--11--07--13--08--15--09
848 * | | | | |
849 * 01 00 03 01 05 02 07 03 08
850 * | | | | |
851 * 00--00--01--02--02--04--03--06--04
852 *
853 * The mask looks as follows (# = masked points)
854 *
855 * +-------+-------+-------+-------+
856 * | | | | |
857 * | | | | |
858 * | | | | |
859 * +-------+-------+-------+-------+
860 * | | | | |
861 * | | | | |
862 * | | | | |
863 * +-------+-------#---#---#---#---#
864 * | | | | |
865 * | | # # # # #
866 * | | | | |
867 * +-------+-------#---#---#---#---#
868 * | | | | |
869 * | | | | |
870 * | | | | |
871 * +-------+-------+-------+-------+
872 */
873static void utest_submain_1(
874 MPI_Comm comp_comm, enum yac_interp_weights_reorder_type reorder_type) {
875
876 char const local_grid_name[] = "grid1";
877 char const remote_grid_name[] = "grid2";
878
879 int my_rank;
880 MPI_Comm_rank(comp_comm, &my_rank);
881
882 double vertex_coordinates_x[] = {0.5, 1.5, 2.5, 3.5, 4.5};
883 double vertex_coordinates_y[] = {0.5, 1.5, 2.5, 3.5, 4.5};
884 double cell_coordinates_x[] = {0.75, 1.75, 2.75, 3.75};
885 double cell_coordinates_y[] = {0.75, 1.75, 2.75, 3.75};
886 double cell_coords[16][3];
887 size_t const num_global_cells[2] = {4,4};
888 size_t local_start[2][2] = {{0,0},{0,2}};
889 size_t local_count[2][2] = {{4,2},{4,2}};
890 int with_halo = 1;
891 for (size_t i = 0; i <= num_global_cells[0]; ++i)
892 vertex_coordinates_x[i] *= YAC_RAD;
893 for (size_t i = 0; i <= num_global_cells[1]; ++i)
894 vertex_coordinates_y[i] *= YAC_RAD;
895 for (size_t i = 0, k = 0; i < num_global_cells[1]; ++i)
896 for (size_t j = 0; j < num_global_cells[0]; ++j, ++k)
898
900 utest_generate_basic_grid_data_reg2d(
901 vertex_coordinates_x, vertex_coordinates_y, num_global_cells,
902 local_start[my_rank], local_count[my_rank], with_halo);
903 size_t num_vertices = grid_data.num_vertices;
904 size_t num_cells = grid_data.num_cells;
905
906 struct yac_basic_grid * local_grid =
907 yac_basic_grid_new(local_grid_name, grid_data);
908 struct yac_basic_grid * remote_grid =
909 yac_basic_grid_empty_new(remote_grid_name);
910
911 yac_int masked_corner_ids[] = {7,8,9,12,13,14};
912 size_t num_masked_corners =
913 sizeof(masked_corner_ids)/sizeof(masked_corner_ids[0]);
914
915 int * corner_mask = xmalloc(num_vertices * sizeof(*corner_mask));
916 for (size_t i = 0; i < num_vertices; ++i) {
917 corner_mask[i] = 1;
918 for (size_t j = 0; j < num_masked_corners; ++j)
919 if (grid_data.vertex_ids[i] == masked_corner_ids[j]) corner_mask[i] = 0;
920 }
921
922 yac_int masked_cell_ids[] = {6,7};
923 size_t num_masked_cells =
924 sizeof(masked_cell_ids)/sizeof(masked_cell_ids[0]);
925
926 int * cell_mask = xmalloc(num_cells * sizeof(*cell_mask));
927 for (size_t i = 0; i < num_cells; ++i) {
928 cell_mask[i] = 1;
929 for (size_t j = 0; j < num_masked_cells; ++j)
930 if (grid_data.cell_ids[i] == masked_cell_ids[j]) cell_mask[i] = 0;
931 }
932
933 yac_coordinate_pointer cell_field_coords =
934 xmalloc(num_cells * sizeof(*cell_field_coords));
935 for (size_t i = 0; i < num_cells; ++i)
936 memcpy(cell_field_coords[i], cell_coords[grid_data.cell_ids[i]],
937 3 * sizeof(double));
938
940 local_grid, YAC_LOC_CORNER, corner_mask, NULL);
942 local_grid, YAC_LOC_CELL, cell_mask, NULL);
944 local_grid, YAC_LOC_CELL, cell_field_coords);
945
946 struct yac_dist_grid_pair * grid_pair =
947 yac_dist_grid_pair_new(local_grid, remote_grid, MPI_COMM_WORLD);
948
949 { // test interpolation without using a mask
950
951 struct yac_interp_field src_fields[] =
952 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX}};
953 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
955 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
956
957 struct yac_interp_grid * interp_grid_in =
958 yac_interp_grid_new(grid_pair, remote_grid_name, local_grid_name,
960 struct yac_interp_grid * interp_grid_out =
961 yac_interp_grid_new(grid_pair, local_grid_name, remote_grid_name,
963
964 struct interp_method * method_stack_in[] =
967 NULL};
968 struct interp_method * method_stack_out[] =
971 NULL};
972
973 struct yac_interp_weights * weights_in =
974 yac_interp_method_do_search(method_stack_in, interp_grid_in);
975 struct yac_interp_weights * weights_out =
976 yac_interp_method_do_search(method_stack_out, interp_grid_out);
977
978 yac_interp_method_delete(method_stack_out);
979 yac_interp_method_delete(method_stack_in);
980
981 struct yac_interpolation * interpolation_in =
983 weights_in, reorder_type, 1,
984 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
985 struct yac_interpolation * interpolation_out =
987 weights_out, reorder_type, 1,
988 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
989
990 //---------------------
991 // do the interpolation
992 //---------------------
993
994 // source_data dimensions [collection_idx]
995 // [pointset_idx]
996 // [local_idx]
997 double * source_data_field =
998 xmalloc(num_vertices * sizeof(*source_data_field));
999 double * source_data_pointset[1] = {source_data_field}; // num_pointset == 1
1000 double ** source_data[1] = {source_data_pointset}; // collection_size == 1
1001 for (size_t i = 0; i < num_vertices; ++i)
1002 source_data_field[i] =
1003 (grid_data.core_vertex_mask[i])?
1004 ((double)(grid_data.vertex_ids[i])):(-1.0);
1005 // target_data dimensions [collection_idx]
1006 // [local_idx]
1007 double * target_data_field =
1008 xmalloc(num_vertices * sizeof(*target_data_field));
1009 double * target_data[1] = {target_data_field}; // collection_size == 1
1010 for (size_t i = 0; i < num_vertices; ++i) target_data_field[i] = -1.0;
1011
1012 yac_interpolation_execute_put(interpolation_out, source_data);
1013 yac_interpolation_execute_get(interpolation_in, target_data);
1014
1015 //----------------------------
1016 // check interpolation results
1017 //----------------------------
1018
1019 double ref_global_target_data[] = { 3, 4, 5, 6,1337,
1020 8, 9, 10, 11,1337,
1021 13, 14, 15, 16,1337,
1022 18, 19, 20, 21,1337,
1023 1337,1337,1337,1337,1337};
1024
1025 for (size_t i = 0; i < num_vertices; ++i) {
1026 if (grid_data.core_vertex_mask[i]) {
1027 if (utest_double_are_unequal(
1028 target_data[0][i],
1029 ref_global_target_data[grid_data.vertex_ids[i]]))
1030 PUT_ERR("error in interpolated data on target side\n");
1031 } else {
1032 if (target_data[0][i] != -1.0)
1033 PUT_ERR("error in interpolated data on target side\n");
1034 }
1035 }
1036
1037 //--------
1038 // cleanup
1039 //--------
1040
1041 free(target_data_field);
1042 free(source_data_field);
1043 yac_interpolation_delete(interpolation_out);
1044 yac_interpolation_delete(interpolation_in);
1045 yac_interp_weights_delete(weights_out);
1046 yac_interp_weights_delete(weights_in);
1047 yac_interp_grid_delete(interp_grid_out);
1048 yac_interp_grid_delete(interp_grid_in);
1049 }
1050
1051 { // test interpolation using a source mask
1052
1053 struct yac_interp_field src_fields[] =
1054 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
1055 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1057 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
1058
1059 struct yac_interp_grid * interp_grid_in =
1060 yac_interp_grid_new(grid_pair, remote_grid_name, local_grid_name,
1062 struct yac_interp_grid * interp_grid_out =
1063 yac_interp_grid_new(grid_pair, local_grid_name, remote_grid_name,
1065
1066 struct interp_method * method_stack_in[] =
1069 NULL};
1070 struct interp_method * method_stack_out[] =
1073 NULL};
1074
1075 struct yac_interp_weights * weights_in =
1076 yac_interp_method_do_search(method_stack_in, interp_grid_in);
1077 struct yac_interp_weights * weights_out =
1078 yac_interp_method_do_search(method_stack_out, interp_grid_out);
1079
1080 yac_interp_method_delete(method_stack_out);
1081 yac_interp_method_delete(method_stack_in);
1082
1083 struct yac_interpolation * interpolation_in =
1085 weights_in, reorder_type, 1,
1086 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1087 struct yac_interpolation * interpolation_out =
1089 weights_out, reorder_type, 1,
1090 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1091
1092 //---------------------
1093 // do the interpolation
1094 //---------------------
1095
1096 // source_data dimensions [collection_idx]
1097 // [pointset_idx]
1098 // [local_idx]
1099 double * source_data_field =
1100 xmalloc(num_vertices * sizeof(*source_data_field));
1101 double * source_data_pointset[1] = {source_data_field}; // num_pointset == 1
1102 double ** source_data[1] = {source_data_pointset}; // collection_size == 1
1103 for (size_t i = 0; i < num_vertices; ++i)
1104 source_data_field[i] =
1105 (grid_data.core_vertex_mask[i])?
1106 ((double)(grid_data.vertex_ids[i])):(-1.0);
1107 // target_data dimensions [collection_idx]
1108 // [local_idx]
1109 double * target_data_field =
1110 xmalloc(num_vertices * sizeof(*target_data_field));
1111 double * target_data[1] = {target_data_field}; // collection_size == 1
1112 for (size_t i = 0; i < num_vertices; ++i) target_data_field[i] = -1.0;
1113
1114 yac_interpolation_execute_put(interpolation_out, source_data);
1115 yac_interpolation_execute_get(interpolation_in, target_data);
1116
1117 //----------------------------
1118 // check interpolation results
1119 //----------------------------
1120
1121 double ref_global_target_data[] = { 3, 4, 5, 6,1337,
1122 1337,1337,1337, 11,1337,
1123 1337,1337,1337, 16,1337,
1124 1337,1337,1337, 21,1337,
1125 1337,1337,1337,1337,1337};
1126
1127 for (size_t i = 0; i < num_vertices; ++i) {
1128 if (grid_data.core_vertex_mask[i]) {
1129 if (utest_double_are_unequal(
1130 target_data[0][i],
1131 ref_global_target_data[grid_data.vertex_ids[i]]))
1132 PUT_ERR("error in interpolated data on target side\n");
1133 } else {
1134 if (target_data[0][i] != -1.0)
1135 PUT_ERR("error in interpolated data on target side\n");
1136 }
1137 }
1138
1139 //--------
1140 // cleanup
1141 //--------
1142
1143 free(target_data_field);
1144 free(source_data_field);
1145 yac_interpolation_delete(interpolation_out);
1146 yac_interpolation_delete(interpolation_in);
1147 yac_interp_weights_delete(weights_out);
1148 yac_interp_weights_delete(weights_in);
1149 yac_interp_grid_delete(interp_grid_out);
1150 yac_interp_grid_delete(interp_grid_in);
1151 }
1152
1153 { // test interpolation using a source mask and allowing partial coverage
1154
1155 struct yac_interp_field src_fields[] =
1156 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
1157 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1159 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
1160
1161 struct yac_interp_grid * interp_grid_in =
1162 yac_interp_grid_new(grid_pair, remote_grid_name, local_grid_name,
1164 struct yac_interp_grid * interp_grid_out =
1165 yac_interp_grid_new(grid_pair, local_grid_name, remote_grid_name,
1167
1168 struct interp_method * method_stack_in[] =
1171 NULL};
1172 struct interp_method * method_stack_out[] =
1175 NULL};
1176
1177 struct yac_interp_weights * weights_in =
1178 yac_interp_method_do_search(method_stack_in, interp_grid_in);
1179 struct yac_interp_weights * weights_out =
1180 yac_interp_method_do_search(method_stack_out, interp_grid_out);
1181
1182 yac_interp_method_delete(method_stack_out);
1183 yac_interp_method_delete(method_stack_in);
1184
1185 struct yac_interpolation * interpolation_in =
1187 weights_in, reorder_type, 1,
1188 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1189 struct yac_interpolation * interpolation_out =
1191 weights_out, reorder_type, 1,
1192 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1193
1194 //---------------------
1195 // do the interpolation
1196 //---------------------
1197
1198 // source_data dimensions [collection_idx]
1199 // [pointset_idx]
1200 // [local_idx]
1201 double * source_data_field =
1202 xmalloc(num_vertices * sizeof(*source_data_field));
1203 double * source_data_pointset[1] = {source_data_field}; // num_pointset == 1
1204 double ** source_data[1] = {source_data_pointset}; // collection_size == 1
1205 for (size_t i = 0; i < num_vertices; ++i)
1206 source_data_field[i] =
1207 (grid_data.core_vertex_mask[i])?
1208 ((double)(grid_data.vertex_ids[i])):(-1.0);
1209 // target_data dimensions [collection_idx]
1210 // [local_idx]
1211 double * target_data_field =
1212 xmalloc(num_vertices * sizeof(*target_data_field));
1213 double * target_data[1] = {target_data_field}; // collection_size == 1
1214 for (size_t i = 0; i < num_vertices; ++i) target_data_field[i] = -1.0;
1215
1216 yac_interpolation_execute_put(interpolation_out, source_data);
1217 yac_interpolation_execute_get(interpolation_in, target_data);
1218
1219 //----------------------------
1220 // check interpolation results
1221 //----------------------------
1222
1223 double ref_global_target_data[] = { 3, 4, 5, 6,1337,
1224 5.5, 6.5,28.0/3.0, 11,1337,
1225 1337,1337, 15.5, 16,1337,
1226 20.5,21.5, 21, 21,1337,
1227 1337,1337, 1337,1337,1337};
1228
1229 for (size_t i = 0; i < num_vertices; ++i) {
1230 if (grid_data.core_vertex_mask[i]) {
1231 if (fabs(target_data[0][i] -
1232 ref_global_target_data[grid_data.vertex_ids[i]]) > 1e-9)
1233 PUT_ERR("error in interpolated data on target side\n");
1234 } else {
1235 if (target_data[0][i] != -1.0)
1236 PUT_ERR("error in interpolated data on target side\n");
1237 }
1238 }
1239
1240 //--------
1241 // cleanup
1242 //--------
1243
1244 free(target_data_field);
1245 free(source_data_field);
1246 yac_interpolation_delete(interpolation_out);
1247 yac_interpolation_delete(interpolation_in);
1248 yac_interp_weights_delete(weights_out);
1249 yac_interp_weights_delete(weights_in);
1250 yac_interp_grid_delete(interp_grid_out);
1251 yac_interp_grid_delete(interp_grid_in);
1252 }
1253
1254 { // test interpolation using a source mask and allowing partial coverage and
1255 // use a target mask
1256
1257 struct yac_interp_grid * interp_grid_in, * interp_grid_out;
1258
1259 {
1260 struct yac_interp_field src_fields[] =
1261 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
1262 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1264 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0};
1265
1266 interp_grid_in =
1267 yac_interp_grid_new(grid_pair, remote_grid_name, local_grid_name,
1268 num_src_fields, src_fields, tgt_field);
1269 }
1270
1271 {
1272 struct yac_interp_field src_fields[] =
1273 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
1274 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1276 {.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0};
1277
1278 interp_grid_out =
1279 yac_interp_grid_new(grid_pair, local_grid_name, remote_grid_name,
1280 num_src_fields, src_fields, tgt_field);
1281 }
1282
1283 struct interp_method * method_stack_in[] =
1286 NULL};
1287 struct interp_method * method_stack_out[] =
1290 NULL};
1291
1292 struct yac_interp_weights * weights_in =
1293 yac_interp_method_do_search(method_stack_in, interp_grid_in);
1294 struct yac_interp_weights * weights_out =
1295 yac_interp_method_do_search(method_stack_out, interp_grid_out);
1296
1297 yac_interp_method_delete(method_stack_out);
1298 yac_interp_method_delete(method_stack_in);
1299
1300 struct yac_interpolation * interpolation_in =
1302 weights_in, reorder_type, 1,
1303 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1304 struct yac_interpolation * interpolation_out =
1306 weights_out, reorder_type, 1,
1307 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1308
1309 //---------------------
1310 // do the interpolation
1311 //---------------------
1312
1313 // source_data dimensions [collection_idx]
1314 // [pointset_idx]
1315 // [local_idx]
1316 double * source_data_field =
1317 xmalloc(num_vertices * sizeof(*source_data_field));
1318 double * source_data_pointset[1] = {source_data_field}; // num_pointset == 1
1319 double ** source_data[1] = {source_data_pointset}; // collection_size == 1
1320 for (size_t i = 0; i < num_vertices; ++i)
1321 source_data_field[i] =
1322 (grid_data.core_vertex_mask[i])?
1323 ((double)(grid_data.vertex_ids[i])):(-1.0);
1324 // target_data dimensions [collection_idx]
1325 // [local_idx]
1326 double * target_data_field =
1327 xmalloc(num_vertices * sizeof(*target_data_field));
1328 double * target_data[1] = {target_data_field}; // collection_size == 1
1329 for (size_t i = 0; i < num_vertices; ++i) target_data_field[i] = -1.0;
1330
1331 yac_interpolation_execute_put(interpolation_out, source_data);
1332 yac_interpolation_execute_get(interpolation_in, target_data);
1333
1334 //----------------------------
1335 // check interpolation results
1336 //----------------------------
1337
1338 double ref_global_target_data[] = { 3, 4, 5, 6,1337,
1339 5.5, 6.5, -1, -1, -1,
1340 1337,1337, -1, -1, -1,
1341 20.5,21.5, 21, 21,1337,
1342 1337,1337,1337,1337,1337};
1343
1344 for (size_t i = 0; i < num_vertices; ++i) {
1345 if ((grid_data.core_vertex_mask[i]) && corner_mask[i]) {
1346 if (fabs(target_data[0][i] -
1347 ref_global_target_data[grid_data.vertex_ids[i]]) > 1e-9)
1348 PUT_ERR("error in interpolated data on target side\n");
1349 } else {
1350 if (target_data[0][i] != -1.0)
1351 PUT_ERR("error in interpolated data on target side\n");
1352 }
1353 }
1354
1355 //--------
1356 // cleanup
1357 //--------
1358
1359 free(target_data_field);
1360 free(source_data_field);
1361 yac_interpolation_delete(interpolation_out);
1362 yac_interpolation_delete(interpolation_in);
1363 yac_interp_weights_delete(weights_out);
1364 yac_interp_weights_delete(weights_in);
1365 yac_interp_grid_delete(interp_grid_out);
1366 yac_interp_grid_delete(interp_grid_in);
1367 }
1368
1369 { // test interpolation using a source mask and allowing partial coverage and
1370 // use a target mask (target points in cell centers)
1371
1372 struct yac_interp_grid * interp_grid_in, * interp_grid_out;
1373 {
1374 struct yac_interp_field src_fields[] =
1375 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
1376 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1378 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
1379 interp_grid_in =
1380 yac_interp_grid_new(grid_pair, remote_grid_name, local_grid_name,
1381 num_src_fields, src_fields, tgt_field);
1382 }
1383 {
1384 struct yac_interp_field src_fields[] =
1385 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
1386 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1388 {.location = YAC_LOC_EDGE, .coordinates_idx = 0, .masks_idx = 0};
1389 interp_grid_out =
1390 yac_interp_grid_new(grid_pair, local_grid_name, remote_grid_name,
1391 num_src_fields, src_fields, tgt_field);
1392 }
1393
1394 struct interp_method * method_stack_in[] =
1397 NULL};
1398 struct interp_method * method_stack_out[] =
1401 NULL};
1402
1403 struct yac_interp_weights * weights_in =
1404 yac_interp_method_do_search(method_stack_in, interp_grid_in);
1405 struct yac_interp_weights * weights_out =
1406 yac_interp_method_do_search(method_stack_out, interp_grid_out);
1407
1408 yac_interp_method_delete(method_stack_out);
1409 yac_interp_method_delete(method_stack_in);
1410
1411 struct yac_interpolation * interpolation_in =
1413 weights_in, reorder_type, 1,
1414 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1415 struct yac_interpolation * interpolation_out =
1417 weights_out, reorder_type, 1,
1418 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1419
1420 //---------------------
1421 // do the interpolation
1422 //---------------------
1423
1424 // source_data dimensions [collection_idx]
1425 // [pointset_idx]
1426 // [local_idx]
1427 double * source_data_field =
1428 xmalloc(num_vertices * sizeof(*source_data_field));
1429 double * source_data_pointset[1] = {source_data_field}; // num_pointset == 1
1430 double ** source_data[1] = {source_data_pointset}; // collection_size == 1
1431 for (size_t i = 0; i < num_vertices; ++i)
1432 source_data_field[i] =
1433 (grid_data.core_vertex_mask[i])?
1434 ((double)(grid_data.vertex_ids[i])):(-1.0);
1435 // target_data dimensions [collection_idx]
1436 // [local_idx]
1437 double * target_data_field =
1438 xmalloc(num_cells * sizeof(*target_data_field));
1439 double * target_data[1] = {target_data_field}; // collection_size == 1
1440 for (size_t i = 0; i < num_cells; ++i) target_data_field[i] = -1.0;
1441
1442 yac_interpolation_execute_put(interpolation_out, source_data);
1443 yac_interpolation_execute_get(interpolation_in, target_data);
1444
1445 //----------------------------
1446 // check interpolation results
1447 //----------------------------
1448
1449 double ref_global_target_data[] = { 3, 4, 5, 6,
1450 5.5, 6.5, -1, -1,
1451 1337,1337,15.5, 16,
1452 20.5,21.5, 21, 21};
1453
1454 for (size_t i = 0; i < num_cells; ++i) {
1455 if ((grid_data.core_cell_mask[i]) && cell_mask[i]) {
1456 if (fabs(target_data[0][i] -
1457 ref_global_target_data[grid_data.cell_ids[i]]) > 1e-9)
1458 PUT_ERR("error in interpolated data on target side\n");
1459 } else {
1460 if (target_data[0][i] != -1.0)
1461 PUT_ERR("error in interpolated data on target side\n");
1462 }
1463 }
1464
1465 //--------
1466 // cleanup
1467 //--------
1468
1469 free(target_data_field);
1470 free(source_data_field);
1471 yac_interpolation_delete(interpolation_out);
1472 yac_interpolation_delete(interpolation_in);
1473 yac_interp_weights_delete(weights_out);
1474 yac_interp_weights_delete(weights_in);
1475 yac_interp_grid_delete(interp_grid_out);
1476 yac_interp_grid_delete(interp_grid_in);
1477 }
1478
1479 { // test interpolation using a fractional mask
1480
1481 struct yac_interp_grid * interp_grid;
1482 {
1483 struct yac_interp_field src_fields[] =
1484 {{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = 0}};
1485 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1487 {.location = YAC_LOC_EDGE, .coordinates_idx = 0, .masks_idx = 0};
1488 interp_grid =
1489 yac_interp_grid_new(grid_pair, local_grid_name, remote_grid_name,
1490 num_src_fields, src_fields, tgt_field);
1491 }
1492
1493 struct interp_method * method_stack[] =
1496 NULL};
1497
1498 struct yac_interp_weights * weights =
1499 yac_interp_method_do_search(method_stack, interp_grid);
1500
1501 yac_interp_method_delete(method_stack);
1502
1503 double frac_mask_value = -1337.0;
1505 interpolations[0] =
1507 weights, reorder_type, 1, frac_mask_value, 1.0, 0.0, NULL, 1, 1);
1508 interpolations[1] =
1510
1511 //---------------------
1512 // do the interpolation
1513 //---------------------
1514
1515 // source_data dimensions [collection_idx]
1516 // [pointset_idx]
1517 // [local_idx]
1518 double * source_data_field =
1519 xmalloc(num_vertices * sizeof(*source_data_field));
1520 double * source_data_pointset[1] = {source_data_field}; // num_pointset == 1
1521 double ** source_data[1] = {source_data_pointset}; // collection_size == 1
1522 double * source_frac_mask_data =
1523 xmalloc(num_vertices * sizeof(*source_frac_mask_data));
1524 double * source_frac_mask_pointset[1] = {source_frac_mask_data};
1525 double ** source_frac_mask[1] = {source_frac_mask_pointset};
1526 for (size_t i = 0; i < num_vertices; ++i) source_frac_mask_data[i] = 0.5;
1527 for (size_t i = 0; i < num_vertices; ++i)
1528 source_data_field[i] =
1529 (grid_data.core_vertex_mask[i])?
1530 ((double)(grid_data.vertex_ids[i])*source_frac_mask_data[i]):(-1.0);
1531
1532 for (int interp_idx = 0; interp_idx < 2; ++interp_idx)
1534 interpolations[interp_idx], source_data, source_frac_mask);
1535
1536 //--------
1537 // cleanup
1538 //--------
1539
1540 free(source_frac_mask_data);
1541 free(source_data_field);
1542 for (int interp_idx = 0; interp_idx < 2; ++interp_idx)
1545 yac_interp_grid_delete(interp_grid);
1546 }
1547
1548 //---------
1549 // clean up
1550 //---------
1551
1552 yac_dist_grid_pair_delete(grid_pair);
1553 yac_basic_grid_delete(remote_grid);
1554 yac_basic_grid_delete(local_grid);
1555}
1556
struct yac_basic_grid * yac_basic_grid_new(char const *name, struct yac_basic_grid_data grid_data)
Definition basic_grid.c:57
size_t yac_basic_grid_add_coordinates_nocpy(struct yac_basic_grid *grid, enum yac_location location, yac_coordinate_pointer coordinates)
Definition basic_grid.c:202
struct yac_basic_grid * yac_basic_grid_empty_new(char const *name)
Definition basic_grid.c:70
void yac_basic_grid_delete(struct yac_basic_grid *grid)
Definition basic_grid.c:77
size_t yac_basic_grid_add_mask_nocpy(struct yac_basic_grid *grid, enum yac_location location, int const *mask, char const *mask_name)
Definition basic_grid.c:244
void yac_dist_grid_pair_delete(struct yac_dist_grid_pair *grid_pair)
Definition dist_grid.c:2377
struct yac_dist_grid_pair * yac_dist_grid_pair_new(struct yac_basic_grid *grid_a, struct yac_basic_grid *grid_b, MPI_Comm comm)
Definition dist_grid.c:2089
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:278
void yac_interp_grid_delete(struct yac_interp_grid *interp_grid)
struct yac_interp_grid * yac_interp_grid_new(struct yac_dist_grid_pair *grid_pair, char const *src_grid_name, char const *tgt_grid_name, size_t num_src_fields, struct yac_interp_field const *src_fields, struct yac_interp_field const tgt_field)
Definition interp_grid.c:31
void yac_interp_method_delete(struct interp_method **method)
Delete an interpolation stack and free its resources (but not the pointer array).
struct yac_interp_weights * yac_interp_method_do_search(struct interp_method **method, struct yac_interp_grid *interp_grid)
Perform weight computation using given interpolation stack and grid.
Defines the interface of the interpolation method "base class" in YAC.
struct interp_method * yac_interp_method_avg_new(enum yac_interp_avg_weight_type weight_type, int partial_coverage)
@ YAC_INTERP_AVG_ARITHMETIC
struct interp_method * yac_interp_method_fixed_new(double value)
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
struct yac_interpolation * yac_interpolation_copy(struct yac_interpolation *interp)
Create a deep copy of an interpolation object.
void yac_interpolation_delete(struct yac_interpolation *interp)
Free an interpolation object and release all resources.
void yac_interpolation_execute_get(struct yac_interpolation *interp, double **tgt_field)
Complete interpolation and write results to the target field (get phase).
void yac_interpolation_execute_put_frac(struct yac_interpolation *interp, double ***src_fields, double ***src_frac_masks)
Provide source field data with fractional masks and start asynchronous execution of interpolation (pu...
void yac_interpolation_execute_put(struct yac_interpolation *interp, double ***src_fields)
Provide source field data and start asynchronous execution of interpolation (put phase).
double const YAC_FRAC_MASK_NO_VALUE
@ YAC_LOC_CORNER
Definition location.h:15
@ YAC_LOC_EDGE
Definition location.h:16
@ YAC_LOC_CELL
Definition location.h:14
#define xmalloc(size)
Definition ppm_xfuncs.h:66
enum yac_location location
Definition basic_grid.h:16
struct yac_interp_field tgt_field
Definition interp_grid.c:26
size_t num_src_fields
Definition interp_grid.c:27
struct yac_dist_grid_pair * grid_pair
Definition interp_grid.c:25
struct yac_interp_field src_fields[]
Definition interp_grid.c:28
int * cell_mask
static MPI_Comm split_comm
double coordinates_x[]
size_t num_cells[2]
double cell_coordinates_y[]
double coordinates_y[]
double cell_coords[36][3]
double cell_coordinates_x[]
#define TEST_EXIT_CODE
Definition tests.h:15
#define PUT_ERR(string)
Definition tests.h:10
struct @6 interpolations[]
YAC_INT yac_int
Definition yac_types.h:15
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:21