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