YAC 3.12.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_interp_method_spmap_parallel.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#include <stdlib.h>
6#include <unistd.h>
7
8#include "tests.h"
9#include "test_common.h"
10#include "geometry.h"
14#include "dist_grid_utils.h"
15#include "yac_mpi.h"
16#include "io_utils.h"
17#include "point_selection.h"
18
19#include <mpi.h>
20#include <yaxt.h>
21#include <netcdf.h>
22
32size_t num_weight_types = sizeof(weight_types) / sizeof(weight_types[0]);
33
40size_t num_scale_types = sizeof(scale_types) / sizeof(scale_types[0]);
41
46size_t num_reorder_types = sizeof(reorder_types) / sizeof(reorder_types[0]);
47
48static char const * grid_names[2] = {"grid1", "grid2"};
49
50static double utest_get_sq_sphere_radius(
51 struct yac_spmap_cell_area_config const * cell_area_confi);
52static double utest_compute_scale(
53 enum yac_interp_spmap_scale_type scale_type,
54 double src_cell_area, double tgt_cell_area);
55static void utest_write_area_file(
56 char const * filename, double * cell_areas, size_t dim_x, size_t dim_y);
57
58int main(void) {
59
60 MPI_Init(NULL, NULL);
61
62 xt_initialize(MPI_COMM_WORLD);
63
64 int comm_rank, comm_size;
65 MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
66 MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
67 MPI_Barrier(MPI_COMM_WORLD);
68
69 if (comm_size != 6) {
70 PUT_ERR("ERROR: wrong number of processes");
71 xt_finalize();
72 MPI_Finalize();
73 return TEST_EXIT_CODE;
74 }
75
76 MPI_Comm split_comm;
77 MPI_Comm_split(MPI_COMM_WORLD, comm_rank < 1, 0, &split_comm);
78
79 int split_comm_rank, split_comm_size;
80 MPI_Comm_rank(split_comm, &split_comm_rank);
81 MPI_Comm_size(split_comm, &split_comm_size);
82
83 { // parallel interpolation process
84 // corner and cell ids for a 7 x 7 grid
85 // 56-----57-----58-----59-----60-----61-----62-----63
86 // | | | | | | | |
87 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
88 // | | | | | | | |
89 // 48-----49-----50-----51-----52-----53-----54-----55
90 // | | | | | | | |
91 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
92 // | | | | | | | |
93 // 40-----41-----42-----43-----44-----45-----46-----47
94 // | | | | | | | |
95 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
96 // | | | | | | | |
97 // 32-----33-----34-----35-----36-----37-----38-----39
98 // | | | | | | | |
99 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
100 // | | | | | | | |
101 // 24-----25-----26-----27-----28-----29-----30-----31
102 // | | | | | | | |
103 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
104 // | | | | | | | |
105 // 16-----17-----18-----19-----20-----21-----22-----23
106 // | | | | | | | |
107 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
108 // | | | | | | | |
109 // 08-----09-----10-----11-----12-----13-----14-----15
110 // | | | | | | | |
111 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
112 // | | | | | | | |
113 // 00-----01-----02-----03-----04-----05-----06-----07
114 //
115 // the grid is distributed among the processes as follows:
116 // (index == process)
117 //
118 // 4---4---4---4---4---4---4---4
119 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
120 // 4---4---4---4---4---4---4---4
121 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
122 // 3---3---3---3---4---4---4---4
123 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
124 // 2---2---2---2---3---3---3---3
125 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
126 // 1---1---1---1---2---2---2---2
127 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
128 // 0---0---0---0---1---1---1---1
129 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
130 // 0---0---0---0---0---0---0---0
131 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
132 // 0---0---0---0---0---0---0---0
133 //
134 // mask
135 // land = 0
136 // coast = 1
137 // ocean = 2
138 // +---+---+---+---+---+---+---+
139 // | 1 | 2 | 2 | 2 | 2 | 2 | 2 |
140 // +---+---+---+---+---+---+---+
141 // | 1 | 2 | 2 | 1 | 1 | 1 | 2 |
142 // +---+---+---+---+---+---+---+
143 // | 1 | 2 | 2 | 1 | 0 | 1 | 2 |
144 // +---+---+---+---+---+---+---+
145 // | 1 | 2 | 2 | 1 | 1 | 1 | 2 |
146 // +---+---+---+---+---+---+---+
147 // | 1 | 1 | 2 | 2 | 2 | 2 | 2 |
148 // +---+---+---+---+---+---+---+
149 // | 0 | 1 | 1 | 2 | 2 | 2 | 2 |
150 // +---+---+---+---+---+---+---+
151 // | 0 | 0 | 1 | 1 | 1 | 1 | 1 |
152 // +---+---+---+---+---+---+---+
153 //
154 //---------------
155 // setup
156 //---------------
157
158 int is_tgt = split_comm_size == 1;
159 double coordinates_x[8] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
160 double coordinates_y[8] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
161 size_t const num_cells[2] = {7,7};
162 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
163 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{7,7}}};
164 int global_mask[7*7] = {
165 0,0,1,1,1,1,1,
166 0,1,1,2,2,2,2,
167 1,1,2,2,2,2,2,
168 1,2,2,1,1,1,2,
169 1,2,2,1,0,1,2,
170 1,2,2,1,1,1,2,
171 1,2,2,2,2,2,2};
172 int with_halo = 0;
173 for (size_t i = 0; i <= num_cells[0]; ++i)
175 for (size_t i = 0; i <= num_cells[1]; ++i)
177
181 local_start[is_tgt][split_comm_rank],
182 local_count[is_tgt][split_comm_rank], with_halo);
183
184 struct yac_basic_grid * grids[2] =
187
188 {
189 int valid_mask_value = (is_tgt)?2:1;
190 yac_coordinate_pointer point_coordinates =
191 xmalloc(grid_data.num_cells * sizeof(*point_coordinates));
192 int * mask = xmalloc(grid_data.num_cells * sizeof(*mask));
193 for (size_t i = 0; i < grid_data.num_cells; ++i) {
194 double * middle_point = point_coordinates[i];
195 for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
196 size_t * curr_vertices =
197 grid_data.cell_to_vertex + grid_data.cell_to_vertex_offsets[i];
198 size_t curr_num_vertices = grid_data.num_vertices_per_cell[i];
199 for (size_t j = 0; j < curr_num_vertices; ++j) {
200 double * curr_vertex_coord =
201 grid_data.vertex_coordinates[curr_vertices[j]];
202 for (size_t k = 0; k < 3; ++k)
203 middle_point[k] += curr_vertex_coord[k];
204 }
205 normalise_vector(middle_point);
206 mask[i] = global_mask[grid_data.cell_ids[i]] == valid_mask_value;
207 }
209 grids[0], YAC_LOC_CELL, point_coordinates);
211 grids[0], YAC_LOC_CELL, mask, NULL);
212 }
213
214 struct yac_dist_grid_pair * grid_pair =
215 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
216
217 struct yac_interp_field src_fields[] =
218 {{.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0}};
219 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
221 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
222
223 struct yac_interp_grid * interp_grid =
226
227 struct interp_method * method_stack[] =
230 yac_interp_method_fixed_new(-2.0), NULL};
231
232 struct yac_interp_weights * weights =
233 yac_interp_method_do_search(method_stack, interp_grid);
234
235 for (size_t i = 0; i < num_reorder_types; ++i) {
236
237 struct yac_interpolation * interpolation =
239 weights, reorder_types[i], 1,
240 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
241
242 // check generated interpolation
243 {
244 double * src_field = NULL;
245 double ** src_fields = &src_field;
246 double * tgt_field = NULL;
247
248 if (is_tgt) {
249 tgt_field = xmalloc(grid_data.num_cells * sizeof(*tgt_field));
250 for (size_t i = 0; i < grid_data.num_cells; ++i) tgt_field[i] = -1;
251 } else {
252 src_field = xmalloc(grid_data.num_cells * sizeof(*src_field));
253 for (size_t i = 0; i < grid_data.num_cells; ++i)
254 src_field[i] = (double)(grid_data.cell_ids[i]);
255 }
256
257 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
258
259 double ref_tgt_field[7*7] = {
260 -1,-1,-1,-1,-1,-1,-1,
261 -1,-1,-1,2+3+9,4,5,6,
262 -1,-1,8+15,-2,25,-2,-2,
263 -1,14+21,24,-1,-1,-1,26,
264 -1,28,31,-1,-1,-1,33,
265 -1,35,38,-1,-1,-1,40,
266 -1,42,-2,-2,39,-2,-2};
267
268 if (is_tgt)
269 for (size_t i = 0; i < grid_data.num_cells; ++i)
270 if (ref_tgt_field[grid_data.cell_ids[i]] != tgt_field[i])
271 PUT_ERR("wrong interpolation result");
272
273 free(tgt_field);
274 free(src_field);
275 }
276
277 yac_interpolation_delete(interpolation);
278 }
279
281 yac_interp_method_delete(method_stack);
282 yac_interp_grid_delete(interp_grid);
283 yac_dist_grid_pair_delete(grid_pair);
286 }
287
288 { // parallel interpolation process
289 // corner and cell ids for a 7 x 7 grid
290 // 56-----57-----58-----59-----60-----61-----62-----63
291 // | | | | | | | |
292 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
293 // | | | | | | | |
294 // 48-----49-----50-----51-----52-----53-----54-----55
295 // | | | | | | | |
296 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
297 // | | | | | | | |
298 // 40-----41-----42-----43-----44-----45-----46-----47
299 // | | | | | | | |
300 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
301 // | | | | | | | |
302 // 32-----33-----34-----35-----36-----37-----38-----39
303 // | | | | | | | |
304 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
305 // | | | | | | | |
306 // 24-----25-----26-----27-----28-----29-----30-----31
307 // | | | | | | | |
308 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
309 // | | | | | | | |
310 // 16-----17-----18-----19-----20-----21-----22-----23
311 // | | | | | | | |
312 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
313 // | | | | | | | |
314 // 08-----09-----10-----11-----12-----13-----14-----15
315 // | | | | | | | |
316 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
317 // | | | | | | | |
318 // 00-----01-----02-----03-----04-----05-----06-----07
319 //
320 // the grid is distributed among the processes as follows:
321 // (index == process)
322 //
323 // 4---4---4---4---4---4---4---4
324 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
325 // 4---4---4---4---4---4---4---4
326 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
327 // 3---3---3---3---4---4---4---4
328 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
329 // 2---2---2---2---3---3---3---3
330 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
331 // 1---1---1---1---2---2---2---2
332 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
333 // 0---0---0---0---1---1---1---1
334 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
335 // 0---0---0---0---0---0---0---0
336 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
337 // 0---0---0---0---0---0---0---0
338 //
339 // mask
340 // land = 0
341 // coast = 1
342 // ocean = 2
343 // +---+---+---+---+---+---+---+
344 // | 1 | 2 | 2 | 2 | 2 | 2 | 2 |
345 // +---+---+---+---+---+---+---+
346 // | 1 | 2 | 2 | 1 | 1 | 1 | 2 |
347 // +---+---+---+---+---+---+---+
348 // | 1 | 2 | 2 | 1 | 0 | 1 | 2 |
349 // +---+---+---+---+---+---+---+
350 // | 1 | 2 | 2 | 1 | 1 | 1 | 2 |
351 // +---+---+---+---+---+---+---+
352 // | 1 | 1 | 2 | 2 | 2 | 2 | 2 |
353 // +---+---+---+---+---+---+---+
354 // | 0 | 1 | 1 | 2 | 2 | 2 | 2 |
355 // +---+---+---+---+---+---+---+
356 // | 0 | 0 | 1 | 1 | 1 | 1 | 1 |
357 // +---+---+---+---+---+---+---+
358 //
359 //---------------
360 // setup
361 //---------------
362
363 int is_tgt = split_comm_size == 1;
364 double coordinates_x[8] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
365 double coordinates_y[8] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
366 size_t const num_cells[2] = {7,7};
367 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
368 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{7,7}}};
369 int global_mask[7*7] = {
370 0,0,1,1,1,1,1,
371 0,1,1,2,2,2,2,
372 1,1,2,2,2,2,2,
373 1,2,2,1,1,1,2,
374 1,2,2,1,0,1,2,
375 1,2,2,1,1,1,2,
376 1,2,2,2,2,2,2};
377 int with_halo = 0;
378 for (size_t i = 0; i <= num_cells[0]; ++i)
380 for (size_t i = 0; i <= num_cells[1]; ++i)
382
386 local_start[is_tgt][split_comm_rank],
387 local_count[is_tgt][split_comm_rank], with_halo);
388
389 struct yac_basic_grid * grids[2] =
392
393 {
394 int valid_mask_value = (is_tgt)?2:1;
395 yac_coordinate_pointer point_coordinates =
396 xmalloc(grid_data.num_cells * sizeof(*point_coordinates));
397 int * mask = xmalloc(grid_data.num_cells * sizeof(*mask));
398 for (size_t i = 0; i < grid_data.num_cells; ++i) {
399 double * middle_point = point_coordinates[i];
400 for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
401 size_t * curr_vertices =
402 grid_data.cell_to_vertex + grid_data.cell_to_vertex_offsets[i];
403 size_t curr_num_vertices = grid_data.num_vertices_per_cell[i];
404 for (size_t j = 0; j < curr_num_vertices; ++j) {
405 double * curr_vertex_coord =
406 grid_data.vertex_coordinates[curr_vertices[j]];
407 for (size_t k = 0; k < 3; ++k)
408 middle_point[k] += curr_vertex_coord[k];
409 }
410 normalise_vector(middle_point);
411 mask[i] = global_mask[grid_data.cell_ids[i]] == valid_mask_value;
412 }
414 grids[0], YAC_LOC_CELL, point_coordinates);
416 grids[0], YAC_LOC_CELL, mask, NULL);
417 }
418
419 struct yac_dist_grid_pair * grid_pair =
420 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
421
422 struct yac_interp_field src_fields[] =
423 {{.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0}};
424 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
426 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
427
428 struct yac_interp_grid * interp_grid =
431
432 struct yac_interp_spmap_config * spmap_config =
434 YAC_RAD * 1.1,
438 struct interp_method * method_stack[] =
441 yac_interp_method_fixed_new(-2.0), NULL};
442 yac_interp_spmap_config_delete(spmap_config);
443
444 struct yac_interp_weights * weights =
445 yac_interp_method_do_search(method_stack, interp_grid);
446
447 for (size_t i = 0; i < num_reorder_types; ++i) {
448
449 struct yac_interpolation * interpolation =
451 weights, reorder_types[i], 1,
452 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
453
454 // check generated interpolation
455 {
456 double * src_field = NULL;
457 double ** src_fields = &src_field;
458 double * tgt_field = NULL;
459
460 if (is_tgt) {
461 tgt_field = xmalloc(grid_data.num_cells * sizeof(*tgt_field));
462 for (size_t i = 0; i < grid_data.num_cells; ++i) tgt_field[i] = -1;
463 } else {
464 src_field = xmalloc(grid_data.num_cells * sizeof(*src_field));
465 for (size_t i = 0; i < grid_data.num_cells; ++i)
466 src_field[i] = (double)(grid_data.cell_ids[i]);
467 }
468
469 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
470
471 double ref_tgt_field[7*7] = {
472 -1,-1,-1,-1,-1,-1,-1,
473 -1,-1,-1, 0, 0, 0, 0,
474 -1,-1, 0, 0, 0, 0, 0,
475 -1, 0, 0,-1,-1,-1, 0,
476 -1, 0, 0,-1,-1,-1, 0,
477 -1, 0, 0,-1,-1,-1, 0,
478 -1, 0, 0, 0, 0, 0, 0};
479 size_t coast_point[] = {2,3,4,5,6,
480 8,9,
481 14,15,
482 21,24,25,26,
483 28,31,33,
484 35,38,39,40,
485 42};
486 size_t num_coast_points = sizeof(coast_point)/sizeof(coast_point[0]);
487 size_t num_tgt_per_coast[] = {3,3,4,4,3,
488 3,3,
489 3,3,
490 3,4,4,3,
491 4,4,3,
492 4,4,3,3,
493 3};
494 size_t tgts[] = {10,11,17, 10,11,17, 10,11,12,18, 11,12,13,19, 12,13,20,
495 16,17,23, 10,11,17,
496 22,23,29, 16,17,23,
497 22,23,29, 16,22,23,30, 11,17,18,19, 20,27,34,
498 22,29,30,36, 23,29,30,37, 27,34,41,
499 29,36,37,43, 30,36,37,44, 45,46,47, 34,41,48,
500 36,43,44};
501
502 for (size_t i = 0, k = 0; i < num_coast_points; ++i) {
503 size_t curr_num_tgt = num_tgt_per_coast[i];
504 double curr_data = (double)(coast_point[i]) / (double)curr_num_tgt;
505 for (size_t j = 0; j < curr_num_tgt; ++j, ++k)
506 ref_tgt_field[tgts[k]] += curr_data;
507 }
508
509 if (is_tgt)
510 for (size_t i = 0; i < grid_data.num_cells; ++i)
511 if (fabs(ref_tgt_field[grid_data.cell_ids[i]] -
512 tgt_field[i]) > 1e-6)
513 PUT_ERR("wrong interpolation result");
514
515 free(tgt_field);
516 free(src_field);
517 }
518
519 yac_interpolation_delete(interpolation);
520 }
521
523 yac_interp_method_delete(method_stack);
524 yac_interp_grid_delete(interp_grid);
525 yac_dist_grid_pair_delete(grid_pair);
528 }
529
530 { // parallel interpolation process
531 // corner and cell ids for a 7 x 9 grid
532 // 72-----73-----74-----75-----76-----77-----78-----79
533 // | | | | | | | |
534 // | 56 | 57 | 58 | 59 | 60 | 61 | 62 |
535 // | | | | | | | |
536 // 64-----65-----66-----67-----68-----69-----70-----71
537 // | | | | | | | |
538 // | 49 | 50 | 51 | 52 | 53 | 54 | 55 |
539 // | | | | | | | |
540 // 56-----57-----58-----59-----60-----61-----62-----63
541 // | | | | | | | |
542 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
543 // | | | | | | | |
544 // 48-----49-----50-----51-----52-----53-----54-----55
545 // | | | | | | | |
546 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
547 // | | | | | | | |
548 // 40-----41-----42-----43-----44-----45-----46-----47
549 // | | | | | | | |
550 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
551 // | | | | | | | |
552 // 32-----33-----34-----35-----36-----37-----38-----39
553 // | | | | | | | |
554 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
555 // | | | | | | | |
556 // 24-----25-----26-----27-----28-----29-----30-----31
557 // | | | | | | | |
558 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
559 // | | | | | | | |
560 // 16-----17-----18-----19-----20-----21-----22-----23
561 // | | | | | | | |
562 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
563 // | | | | | | | |
564 // 08-----09-----10-----11-----12-----13-----14-----15
565 // | | | | | | | |
566 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
567 // | | | | | | | |
568 // 00-----01-----02-----03-----04-----05-----06-----07
569 //
570 // the grid is distributed among the processes as follows:
571 // (index == process)
572 //
573 // 4---4---4---4---4---4---4---4
574 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
575 // 4---4---4---4---4---4---4---4
576 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
577 // 4---4---4---4---4---4---4---4
578 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
579 // 3---3---3---3---4---4---4---4
580 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
581 // 2---2---2---2---3---3---3---3
582 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
583 // 1---1---1---1---2---2---2---2
584 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
585 // 0---0---0---0---1---1---1---1
586 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
587 // 0---0---0---0---0---0---0---0
588 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
589 // 0---0---0---0---0---0---0---0
590 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
591 // 0---0---0---0---0---0---0---0
592 //
593 // mask
594 // land = 0
595 // coast = 1
596 // ocean = 2
597 // +---+---+---+---+---+---+---+
598 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
599 // +---+---+---+---+---+---+---+
600 // | 1 | 0 | 0 | 0 | 0 | 1 | 2 |
601 // +---+---+---+---+---+---+---+
602 // | 2 | 1 | 0 | 0 | 1 | 2 | 2 |
603 // +---+---+---+---+---+---+---+
604 // | 2 | 2 | 1 | 1 | 2 | 2 | 2 |
605 // +---+---+---+---+---+---+---+
606 // | 2 | 2 | 1 | 1 | 2 | 2 | 2 |
607 // +---+---+---+---+---+---+---+
608 // | 2 | 1 | 0 | 0 | 1 | 2 | 2 |
609 // +---+---+---+---+---+---+---+
610 // | 1 | 0 | 0 | 0 | 0 | 1 | 2 |
611 // +---+---+---+---+---+---+---+
612 // | 0 | 0 | 0 | 0 | 0 | 0 | 1 |
613 // +---+---+---+---+---+---+---+
614 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
615 // +---+---+---+---+---+---+---+
616 //
617 //---------------
618 // setup
619 //---------------
620
621 enum {NUM_X = 7, NUM_Y = 9};
622 int is_tgt = split_comm_size == 1;
623 double coordinates_x[NUM_X+1] = {0.0,1.01,2.03,3.06,4.1,5.15,6.21,7.28};
624 double coordinates_y[NUM_Y+1] =
625 {-1.0,0.0,1.01,2.03,3.06,4.1,5.15,6.21,7.28,8.36};
626 size_t const num_cells[2] = {NUM_X,NUM_Y};
627 size_t local_start[2][5][2] = {{{0,0},{0,3},{0,4},{4,3},{0,6}}, {{0,0}}};
628 size_t local_count[2][5][2] = {{{7,3},{4,1},{4,2},{3,3},{7,3}},
629 {{NUM_X,NUM_Y}}};
630 int global_mask[NUM_X * NUM_Y] = {
631 0,0,0,0,0,0,0,
632 0,0,0,0,0,0,1,
633 1,0,0,0,0,1,2,
634 2,1,0,0,1,2,2,
635 2,2,1,1,2,2,2,
636 2,2,1,1,2,2,2,
637 2,1,0,0,1,2,2,
638 1,0,0,0,0,1,2,
639 0,0,0,0,0,0,0};
640 int with_halo = 0;
641 double const grid_scale = 10.0;
642 for (size_t i = 0; i <= NUM_X; ++i)
643 coordinates_x[i] *= YAC_RAD * grid_scale;
644 for (size_t i = 0; i <= NUM_Y; ++i)
645 coordinates_y[i] *= YAC_RAD * grid_scale;
646
650 local_start[is_tgt][split_comm_rank],
651 local_count[is_tgt][split_comm_rank], with_halo);
652
653 struct yac_basic_grid * grids[2] =
656
657 yac_coordinate_pointer point_coordinates =
658 xmalloc(grid_data.num_cells * sizeof(*point_coordinates));
659 {
660 int valid_mask_value = (is_tgt)?2:1;
661 int * mask = xmalloc(grid_data.num_cells * sizeof(*mask));
662 for (size_t i = 0; i < grid_data.num_cells; ++i) {
663 double * middle_point = point_coordinates[i];
664 for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
665 size_t * curr_vertices =
666 grid_data.cell_to_vertex + grid_data.cell_to_vertex_offsets[i];
667 size_t curr_num_vertices = grid_data.num_vertices_per_cell[i];
668 for (size_t j = 0; j < curr_num_vertices; ++j) {
669 double * curr_vertex_coord =
670 grid_data.vertex_coordinates[curr_vertices[j]];
671 for (size_t k = 0; k < 3; ++k)
672 middle_point[k] += curr_vertex_coord[k];
673 }
674 normalise_vector(middle_point);
675 mask[i] = global_mask[grid_data.cell_ids[i]] == valid_mask_value;
676 }
678 grids[0], YAC_LOC_CELL, point_coordinates);
680 grids[0], YAC_LOC_CELL, mask, NULL);
681 }
682
683 struct yac_dist_grid_pair * grid_pair =
684 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
685
686 double grid_cell_areas[NUM_X*NUM_Y];
687 for (int i = 0, k = 0; i < NUM_Y; ++i)
688 for (int j = 0; j < NUM_X; ++j, ++k)
689 grid_cell_areas[k] =
690 fabs(
691 (coordinates_x[j+1] - coordinates_x[j+0]) *
692 (sin(coordinates_y[i+0]) - sin(coordinates_y[i+1])));
693
694 char const * area_filename = "test_interp_method_spmap_area.nc";
695 if (comm_rank == 0)
696 utest_write_area_file(area_filename, grid_cell_areas, NUM_X, NUM_Y);
697 MPI_Barrier(MPI_COMM_WORLD);
698
699 struct yac_interp_field src_fields[] =
700 {{.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0}};
701 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
703 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
704
705 struct yac_interp_grid * interp_grid =
708
709 double const spread_distances[] = {3.7 * grid_scale, 0.0, 0.001};
710 enum {NUM_SPREAD_DISTANCES =
711 sizeof(spread_distances) / sizeof(spread_distances[0])};
712
713 struct {
716 } scale_configs[] =
720 .tgt =
721 yac_spmap_cell_area_config_file_new(area_filename, "area_1d", 0)},
722 {.src =
723 yac_spmap_cell_area_config_file_new(area_filename, "area_1d", 0),
725 {.src =
726 yac_spmap_cell_area_config_file_new(area_filename, "area_2d", 0),
727 .tgt =
728 yac_spmap_cell_area_config_file_new(area_filename, "area_2d", 0)},
729 };
730 enum {NUM_SCALE_CONFIGS =
731 sizeof(scale_configs) / sizeof(scale_configs[0])};
732
733 char const * io_ranks[] = {"0", "3,5", "1,3,5"};
734 enum {NUM_IO_CONFIGS = sizeof(io_ranks) / sizeof(io_ranks[0])};
735
736 for (size_t spread_distance_idx = 0;
737 spread_distance_idx < NUM_SPREAD_DISTANCES; ++spread_distance_idx) {
738
739 for (size_t weight_type_idx = 0; weight_type_idx < num_weight_types;
740 ++weight_type_idx) {
741
742 for (size_t scale_type_idx = 0; scale_type_idx < num_scale_types;
743 ++scale_type_idx) {
744
745 for (size_t scale_config_idx = 0;
746 scale_config_idx < NUM_SCALE_CONFIGS; ++scale_config_idx) {
747
748 struct yac_spmap_scale_config * scale_config =
750 scale_types[scale_type_idx],
751 scale_configs[scale_config_idx].src,
752 scale_configs[scale_config_idx].tgt);
753
754 for (size_t io_config_idx = 0; io_config_idx < NUM_IO_CONFIGS;
755 ++io_config_idx) {
756
757 // clear environment
759
760 setenv("YAC_IO_RANK_LIST", io_ranks[io_config_idx], 1);
761 setenv("YAC_IO_MAX_NUM_RANKS_PER_NODE", "12", 1);
762
763 struct yac_interp_spmap_config * spmap_config =
765 YAC_RAD * spread_distances[spread_distance_idx],
767 weight_types[weight_type_idx],
769
770 struct interp_method * method_stack[] =
773 yac_interp_method_fixed_new(-2.0), NULL};
774 yac_interp_spmap_config_delete(spmap_config);
775
776 struct yac_interp_weights * weights =
777 yac_interp_method_do_search(method_stack, interp_grid);
778
779 for (size_t reorder_type_idx = 0;
780 reorder_type_idx < num_reorder_types;
781 ++reorder_type_idx) {
782
783 struct yac_interpolation * interpolation =
785 weights, reorder_types[reorder_type_idx], 1,
786 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
787
788 // check generated interpolation
789 {
790 double * src_field = NULL;
791 double ** src_fields = &src_field;
792 double * tgt_field = NULL;
793
794 if (is_tgt) {
795 tgt_field =
796 xmalloc(grid_data.num_cells * sizeof(*tgt_field));
797 for (size_t i = 0; i < grid_data.num_cells; ++i)
798 tgt_field[i] = -1;
799 } else {
800 src_field =
801 xmalloc(grid_data.num_cells * sizeof(*src_field));
802 for (size_t i = 0; i < grid_data.num_cells; ++i)
803 src_field[i] = (double)(grid_data.cell_ids[i]);
804 }
805
807 interpolation, &src_fields, &tgt_field);
808
809 if (is_tgt) {
810
811 double ref_tgt_field[NUM_X * NUM_Y] = {
812 -1,-1,-1,-1,-1,-1,-1,
813 -1,-1,-1,-1,-1,-1,-1,
814 -1,-1,-1,-1,-1,-1, 0,
815 0,-1,-1,-1,-1, 0, 0,
816 0, 0,-1,-1, 0, 0, 0,
817 0, 0,-1,-1, 0, 0, 0,
818 0,-1,-1,-1,-1, 0, 0,
819 -1,-1,-1,-1,-1,-1, 0,
820 -1,-1,-1,-1,-1,-1,-1};
821 size_t coast_point[] = {14,22,30,37,43,49,
822 13,19,25,31,38,46,54};
823 enum {
824 NUM_COAST_POINTS =
825 sizeof(coast_point)/sizeof(coast_point[0]),
826 MAX_NUM_TGT_PER_COAST = 12};
827 size_t num_tgt_per_coast[NUM_SPREAD_DISTANCES]
828 [NUM_COAST_POINTS] =
829 {{6,6,6,6,6,6,
830 9,9,11,12,12,11,9},
831 {1,1,1,1,1,1,
832 1,1,1,1,1,1,1},
833 {1,1,1,1,1,1,
834 1,1,1,1,1,1,1}};
835 size_t tgts[NUM_SPREAD_DISTANCES]
836 [NUM_COAST_POINTS]
837 [MAX_NUM_TGT_PER_COAST] =
838 {{{21,28,29,35,36,42},
839 {21,28,29,35,36,42},
840 {21,28,29,35,36,42},
841 {21,28,29,35,36,42},
842 {21,28,29,35,36,42},
843 {21,28,29,35,36,42},
844
845 {20,26,27,32,33,34,39,40,41},
846 {20,26,27,32,33,34,39,40,41},
847 {20,26,27,32,33,34,39,40,41,47,48},
848 {20,26,27,32,33,34,39,40,41,47,48,55},
849 {20,26,27,32,33,34,39,40,41,47,48,55},
850 {26,27,32,33,34,39,40,41,47,48,55},
851 {32,33,34,39,40,41,47,48,55}},
852 {{21},
853 {21},
854 {29},
855 {36},
856 {42},
857 {42},
858
859 {20},
860 {20},
861 {26},
862 {32},
863 {39},
864 {47},
865 {55}},
866 {{21},
867 {21},
868 {29},
869 {36},
870 {42},
871 {42},
872
873 {20},
874 {20},
875 {26},
876 {32},
877 {39},
878 {47},
879 {55}}};
880
881 switch (weight_types[weight_type_idx]) {
882 default:
883 case(YAC_INTERP_SPMAP_AVG): {
884 for (size_t i = 0; i < NUM_COAST_POINTS; ++i) {
885 size_t curr_num_tgt =
886 num_tgt_per_coast[spread_distance_idx][i];
887 double curr_data =
888 (double)(coast_point[i]) / (double)curr_num_tgt;
889 for (size_t j = 0; j < curr_num_tgt; ++j)
890 ref_tgt_field[tgts[spread_distance_idx][i][j]] +=
891 curr_data *
892 utest_compute_scale(
893 scale_types[scale_type_idx],
894 grid_cell_areas[coast_point[i]] *
895 utest_get_sq_sphere_radius(
897 scale_config)),
898 grid_cell_areas[
899 tgts[spread_distance_idx][i][j]] *
900 utest_get_sq_sphere_radius(
902 scale_config)));
903 }
904 break;
905 }
906 case(YAC_INTERP_SPMAP_DIST): {
907 for (size_t i = 0, k = 0; i < NUM_COAST_POINTS; ++i) {
908 size_t * curr_tgts = tgts[spread_distance_idx][i];
909 size_t curr_num_tgt =
910 num_tgt_per_coast[spread_distance_idx][i];
911 size_t curr_coast_point = coast_point[i];
912 double curr_src_data = (double)(curr_coast_point);
913 double inv_distances[curr_num_tgt];
914 double inv_distances_sum = 0.0;
915 for (size_t j = 0; j < curr_num_tgt; ++j)
916 inv_distances_sum +=
917 ((inv_distances[j] =
918 1.0 / get_vector_angle(
919 point_coordinates[curr_coast_point],
920 point_coordinates[curr_tgts[j]])));
921 for (size_t j = 0; j < curr_num_tgt; ++j, ++k)
922 ref_tgt_field[curr_tgts[j]] +=
923 curr_src_data
924 * (inv_distances[j] / inv_distances_sum) *
925 utest_compute_scale(
926 scale_types[scale_type_idx],
927 grid_cell_areas[coast_point[i]] *
928 utest_get_sq_sphere_radius(
930 scale_config)),
931 grid_cell_areas[
932 tgts[spread_distance_idx][i][j]] *
933 utest_get_sq_sphere_radius(
935 scale_config)));
936 }
937 }
938 }
939
940 for (size_t i = 0; i < grid_data.num_cells; ++i)
941 if (((ref_tgt_field[grid_data.cell_ids[i]] == 0.0) &&
942 (tgt_field[i] != -2.0)) ||
943 ((ref_tgt_field[grid_data.cell_ids[i]] != 0.0) &&
944 (fabs(ref_tgt_field[grid_data.cell_ids[i]] -
945 tgt_field[i]) > 1.0e-6)))
946 PUT_ERR("wrong interpolation result");
947
948 double src_sum = 0.0;
949 for (size_t i = 0; i < NUM_COAST_POINTS; ++i) {
950 double curr_src_data = (double)(coast_point[i]);
951 if (
952 (scale_types[scale_type_idx] ==
954 (scale_types[scale_type_idx] ==
956 curr_src_data *=
957 grid_cell_areas[coast_point[i]] *
958 utest_get_sq_sphere_radius(
960 scale_config));
961 src_sum += curr_src_data;
962 }
963 double tgt_sum = 0.0;
964 for (size_t i = 0; i < grid_data.num_cells; ++i) {
965 if ((tgt_field[i] != -1.0) &&
966 (tgt_field[i] != -2.0)) {
967 double curr_tgt_data = tgt_field[i];
968 if (
969 (scale_types[scale_type_idx] ==
971 (scale_types[scale_type_idx] ==
973 curr_tgt_data *=
974 grid_cell_areas[grid_data.cell_ids[i]] *
975 utest_get_sq_sphere_radius(
977 scale_config));
978 tgt_sum += curr_tgt_data;
979 }
980 }
981 if (fabs(src_sum - tgt_sum) > 1.0e-6)
982 PUT_ERR("wrong interpolation result (not conservative)");
983 }
984
985 free(tgt_field);
986 free(src_field);
987 } // check
988
989 yac_interpolation_delete(interpolation);
990 } // reorder_type_idx
991
993 yac_interp_method_delete(method_stack);
994 } // io_config_idx
995 yac_spmap_scale_config_delete(scale_config);
996 } // scale_config_idx
997 } // scale_type_idx
998 } // weight_type_idx
999 } // spread_distance_idx
1000
1001 for (size_t i = 0; i < NUM_SCALE_CONFIGS; ++i) {
1002 yac_spmap_cell_area_config_delete(scale_configs[i].src);
1003 yac_spmap_cell_area_config_delete(scale_configs[i].tgt);
1004 }
1005 yac_interp_grid_delete(interp_grid);
1006 yac_dist_grid_pair_delete(grid_pair);
1009
1010 MPI_Barrier(MPI_COMM_WORLD);
1011 if (comm_rank == 0) unlink(area_filename);
1012 }
1013
1014 { // parallel interpolation process
1015 // corner and cell ids for a 7 x 7 grid
1016 // 56-----57-----58-----59-----60-----61-----62-----63
1017 // | | | | | | | |
1018 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
1019 // | | | | | | | |
1020 // 48-----49-----50-----51-----52-----53-----54-----55
1021 // | | | | | | | |
1022 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
1023 // | | | | | | | |
1024 // 40-----41-----42-----43-----44-----45-----46-----47
1025 // | | | | | | | |
1026 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
1027 // | | | | | | | |
1028 // 32-----33-----34-----35-----36-----37-----38-----39
1029 // | | | | | | | |
1030 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
1031 // | | | | | | | |
1032 // 24-----25-----26-----27-----28-----29-----30-----31
1033 // | | | | | | | |
1034 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
1035 // | | | | | | | |
1036 // 16-----17-----18-----19-----20-----21-----22-----23
1037 // | | | | | | | |
1038 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
1039 // | | | | | | | |
1040 // 08-----09-----10-----11-----12-----13-----14-----15
1041 // | | | | | | | |
1042 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
1043 // | | | | | | | |
1044 // 00-----01-----02-----03-----04-----05-----06-----07
1045 //
1046 // the grid is distributed among the processes as follows:
1047 // (index == process)
1048 //
1049 // 4---4---4---4---4---4---4---4
1050 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
1051 // 4---4---4---4---4---4---4---4
1052 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
1053 // 3---3---3---3---4---4---4---4
1054 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
1055 // 2---2---2---2---3---3---3---3
1056 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
1057 // 1---1---1---1---2---2---2---2
1058 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
1059 // 0---0---0---0---1---1---1---1
1060 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
1061 // 0---0---0---0---0---0---0---0
1062 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
1063 // 0---0---0---0---0---0---0---0
1064 //
1065 // mask
1066 // land = 0
1067 // coast = 1
1068 // ocean = 2
1069 // +---+---+---+---+---+---+---+
1070 // | 0 | 0 | 0 | 1 | 2 | 2 | 2 |
1071 // +---+---+---+---+---+---+---+
1072 // | 0 | 0 | 0 | 0 | 1 | 2 | 2 |
1073 // +---+---+---+---+---+---+---+
1074 // | 0 | 0 | 0 | 0 | 1 | 1 | 2 |
1075 // +---+---+---+---+---+---+---+
1076 // | 0 | 0 | 0 | 1 | 0 | 0 | 1 |
1077 // +---+---+---+---+---+---+---+
1078 // | 0 | 0 | 1 | 0 | 0 | 0 | 0 |
1079 // +---+---+---+---+---+---+---+
1080 // | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
1081 // +---+---+---+---+---+---+---+
1082 // | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
1083 // +---+---+---+---+---+---+---+
1084 //
1085 //---------------
1086 // setup
1087 //---------------
1088
1089 int is_tgt = split_comm_size == 1;
1090 double coordinates_x[8] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
1091 double coordinates_y[8] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
1092 size_t const num_cells[2] = {7,7};
1093 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
1094 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{7,7}}};
1095 int global_mask[7*7] = {
1096 1,0,0,0,0,0,0,
1097 0,1,0,0,0,0,0,
1098 0,0,1,0,0,0,0,
1099 0,0,0,1,0,0,1,
1100 0,0,0,0,1,1,2,
1101 0,0,0,0,1,2,2,
1102 0,0,0,1,2,2,2};
1103 int with_halo = 0;
1104 for (size_t i = 0; i <= num_cells[0]; ++i)
1105 coordinates_x[i] *= YAC_RAD;
1106 for (size_t i = 0; i <= num_cells[1]; ++i)
1107 coordinates_y[i] *= YAC_RAD;
1108
1112 local_start[is_tgt][split_comm_rank],
1113 local_count[is_tgt][split_comm_rank], with_halo);
1114
1115 struct yac_basic_grid * grids[2] =
1118
1119 yac_coordinate_pointer point_coordinates =
1120 xmalloc(grid_data.num_cells * sizeof(*point_coordinates));
1121 {
1122 int valid_mask_value = (is_tgt)?2:1;
1123 int * mask = xmalloc(grid_data.num_cells * sizeof(*mask));
1124 for (size_t i = 0; i < grid_data.num_cells; ++i) {
1125 double * middle_point = point_coordinates[i];
1126 for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
1127 size_t * curr_vertices =
1128 grid_data.cell_to_vertex + grid_data.cell_to_vertex_offsets[i];
1129 size_t curr_num_vertices = grid_data.num_vertices_per_cell[i];
1130 for (size_t j = 0; j < curr_num_vertices; ++j) {
1131 double * curr_vertex_coord =
1132 grid_data.vertex_coordinates[curr_vertices[j]];
1133 for (size_t k = 0; k < 3; ++k)
1134 middle_point[k] += curr_vertex_coord[k];
1135 }
1136 normalise_vector(middle_point);
1137 mask[i] = global_mask[grid_data.cell_ids[i]] == valid_mask_value;
1138 }
1140 grids[0], YAC_LOC_CELL, point_coordinates);
1142 grids[0], YAC_LOC_CELL, mask, NULL);
1143 }
1144
1145 struct yac_dist_grid_pair * grid_pair =
1146 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
1147
1148 struct yac_interp_field src_fields[] =
1149 {{.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0}};
1150 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1152 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
1153
1154 struct yac_interp_grid * interp_grid =
1157
1158 struct yac_interp_spmap_config * spmap_config =
1161 YAC_RAD * 3.0,
1164
1165 struct interp_method * method_stack[] =
1168 yac_interp_method_fixed_new(-2.0), NULL};
1169 yac_interp_spmap_config_delete(spmap_config);
1170
1171 struct yac_interp_weights * weights =
1172 yac_interp_method_do_search(method_stack, interp_grid);
1173
1174 for (size_t i = 0; i < num_reorder_types; ++i) {
1175
1176 struct yac_interpolation * interpolation =
1178 weights, reorder_types[i], 1,
1179 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1180
1181 // check generated interpolation
1182 {
1183 double * src_field = NULL;
1184 double ** src_fields = &src_field;
1185 double * tgt_field = NULL;
1186
1187 if (is_tgt) {
1188 tgt_field = xmalloc(grid_data.num_cells * sizeof(*tgt_field));
1189 for (size_t i = 0; i < grid_data.num_cells; ++i) tgt_field[i] = -1;
1190 } else {
1191 src_field = xmalloc(grid_data.num_cells * sizeof(*src_field));
1192 for (size_t i = 0; i < grid_data.num_cells; ++i)
1193 src_field[i] = (double)(grid_data.cell_ids[i]);
1194 }
1195
1196 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
1197
1198 if (is_tgt) {
1199
1200 double ref_tgt_field[7*7] = {
1201 -1,-1,-1,-1,-1,-1,-1,
1202 -1,-1,-1,-1,-1,-1,-1,
1203 -1,-1,-1,-1,-1,-1,-1,
1204 -1,-1,-1,-1,-1,-1,-1,
1205 -1,-1,-1,-1,-1,-1, 0,
1206 -1,-1,-1,-1,-1, 0,-2,
1207 -1,-1,-1,-1, 0,-2,-2};
1208 size_t coast_point[] = {0,8,16,24,27,32,33,39,45};
1209 size_t num_coast_points = sizeof(coast_point)/sizeof(coast_point[0]);
1210 size_t num_tgt_per_coast[] = {0,0,0,1,1,1,1,1,1};
1211 size_t tgts[] = {40,34,40,34,40,46};
1212
1213 for (size_t i = 0, k = 0; i < num_coast_points; ++i) {
1214 size_t curr_num_tgt = num_tgt_per_coast[i];
1215 if (curr_num_tgt == 0) continue;
1216 double curr_data = (double)(coast_point[i]) / (double)curr_num_tgt;
1217 for (size_t j = 0; j < curr_num_tgt; ++j, ++k)
1218 ref_tgt_field[tgts[k]] += curr_data;
1219 }
1220
1221 for (size_t i = 0; i < grid_data.num_cells; ++i)
1222 if (fabs(ref_tgt_field[grid_data.cell_ids[i]] -
1223 tgt_field[i]) > 1e-6)
1224 PUT_ERR("wrong interpolation result");
1225 }
1226
1227 free(tgt_field);
1228 free(src_field);
1229 }
1230
1231 yac_interpolation_delete(interpolation);
1232 }
1233
1235 yac_interp_method_delete(method_stack);
1236
1237 yac_interp_grid_delete(interp_grid);
1238 yac_dist_grid_pair_delete(grid_pair);
1241 }
1242
1243 { // testing YAC_INTERP_SPMAP_DIST
1244 //---------------
1245 // setup
1246 //---------------
1247
1248 int is_tgt = split_comm_size == 1;
1249 double coordinates_x[4] = {0.0,1.0,2.0,3.0};
1250 double coordinates_y[4] = {0.0,1.0,2.0,3.0};
1251 size_t const num_cells[2] = {3,3};
1252 size_t local_start[2][5][2] = {{{0,0},{0,0},{0,0},{0,0},{0,0}}, {{0,0}}};
1253 size_t local_count[2][5][2] = {{{3,3},{3,3},{3,3},{3,3},{3,3}}, {{3,3}}};
1254 int global_mask[3*3] = {
1255 1,1,2,
1256 1,3,2,
1257 2,2,2};
1258 int with_halo = 0;
1259
1260 for (size_t i = 0; i <= num_cells[0]; ++i)
1261 coordinates_x[i] *= YAC_RAD;
1262 for (size_t i = 0; i <= num_cells[1]; ++i)
1263 coordinates_y[i] *= YAC_RAD;
1264
1268 local_start[is_tgt][split_comm_rank],
1269 local_count[is_tgt][split_comm_rank], with_halo);
1270
1271 struct yac_basic_grid * grids[2] =
1274
1275 yac_coordinate_pointer point_coordinates =
1276 xmalloc(grid_data.num_cells * sizeof(*point_coordinates));
1277 {
1278 int valid_mask_value = (is_tgt)?2:1;
1279 int * mask = xmalloc(grid_data.num_cells * sizeof(*mask));
1280 for (size_t i = 0; i < grid_data.num_cells; ++i) {
1281 double * middle_point = point_coordinates[i];
1282 for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
1283 size_t * curr_vertices =
1284 grid_data.cell_to_vertex + grid_data.cell_to_vertex_offsets[i];
1285 size_t curr_num_vertices = grid_data.num_vertices_per_cell[i];
1286 for (size_t j = 0; j < curr_num_vertices; ++j) {
1287 double * curr_vertex_coord =
1288 grid_data.vertex_coordinates[curr_vertices[j]];
1289 for (size_t k = 0; k < 3; ++k)
1290 middle_point[k] += curr_vertex_coord[k];
1291 }
1292 normalise_vector(middle_point);
1293 mask[i] = (global_mask[grid_data.cell_ids[i]] & valid_mask_value) > 0;
1294 }
1296 grids[0], YAC_LOC_CELL, point_coordinates);
1298 grids[0], YAC_LOC_CELL, mask, NULL);
1299 }
1300
1301 struct yac_dist_grid_pair * grid_pair =
1302 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
1303
1304 struct yac_interp_field src_fields[] =
1305 {{.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0}};
1306 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1308 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
1309
1310 struct yac_interp_grid * interp_grid =
1313
1314 struct yac_interp_spmap_config * spmap_config =
1316 1.1 * YAC_RAD,
1320
1321 struct interp_method * method_stack[] =
1324 yac_interp_method_fixed_new(-2.0), NULL};
1325 yac_interp_spmap_config_delete(spmap_config);
1326
1327 struct yac_interp_weights * weights =
1328 yac_interp_method_do_search(method_stack, interp_grid);
1329
1330 struct yac_interpolation * interpolation =
1332 weights, YAC_MAPPING_ON_SRC, 1,
1333 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1334
1335 // check generated interpolation
1336 {
1337 double * src_field = NULL;
1338 double ** src_fields = &src_field;
1339 double * tgt_field = NULL;
1340
1341 if (is_tgt) {
1342 tgt_field = xmalloc(grid_data.num_cells * sizeof(*tgt_field));
1343 for (size_t i = 0; i < grid_data.num_cells; ++i) tgt_field[i] = -1;
1344 } else {
1345 src_field = xmalloc(grid_data.num_cells * sizeof(*src_field));
1346 for (size_t i = 0; i < grid_data.num_cells; ++i)
1347 src_field[i] = (double)(grid_data.cell_ids[i] + 1);
1348 }
1349
1350 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
1351
1352 if (is_tgt) {
1353 double ref_tgt_field[3*3] = {-1,-1,0, -1,0,0, -2,0,-2};
1354 size_t coast_point[] = {0,1,3,4};
1355 enum {
1356 NUM_COAST_POINTS = sizeof(coast_point)/sizeof(coast_point[0]),
1357 MAX_NUM_TGT_PER_COAST = 3};
1358 size_t num_tgt_per_coast[NUM_COAST_POINTS] = {3,2,3,1};
1359 size_t tgts[NUM_COAST_POINTS][MAX_NUM_TGT_PER_COAST] =
1360 {{4,5,7},{2,5},{4,5,7},{4}};
1361
1362
1363 for (size_t i = 0, k = 0; i < NUM_COAST_POINTS; ++i) {
1364 size_t * curr_tgts = tgts[i];
1365 size_t curr_num_tgt = num_tgt_per_coast[i];
1366 size_t curr_coast_point = coast_point[i];
1367 double curr_src_data = (double)(curr_coast_point + 1);
1368 double inv_distances[curr_num_tgt];
1369 double inv_distances_sum = 0.0;
1370 if ( curr_num_tgt == 1) {
1371 ref_tgt_field[curr_tgts[0]] += curr_src_data;
1372 } else {
1373 for (size_t j = 0; j < curr_num_tgt; ++j)
1374 inv_distances_sum +=
1375 ((inv_distances[j] =
1376 1.0 / get_vector_angle(
1377 point_coordinates[curr_coast_point],
1378 point_coordinates[curr_tgts[j]])));
1379 for (size_t j = 0; j < curr_num_tgt; ++j, ++k)
1380 ref_tgt_field[curr_tgts[j]] +=
1381 curr_src_data * (inv_distances[j] / inv_distances_sum);
1382 }
1383 }
1384
1385 for (size_t i = 0; i < grid_data.num_cells; ++i)
1386 if (fabs(ref_tgt_field[grid_data.cell_ids[i]] - tgt_field[i]) > 1e-9)
1387 PUT_ERR("wrong interpolation result");
1388 }
1389
1390 free(tgt_field);
1391 free(src_field);
1392 }
1393
1394 yac_interpolation_delete(interpolation);
1395
1397 yac_interp_method_delete(method_stack);
1398 yac_interp_grid_delete(interp_grid);
1399 yac_dist_grid_pair_delete(grid_pair);
1402 }
1403
1404 { // parallel interpolation process
1405 // corner and cell ids for a 7 x 7 grid
1406 // 56-----57-----58-----59-----60-----61-----62-----63
1407 // | | | | | | | |
1408 // | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
1409 // | | | | | | | |
1410 // 48-----49-----50-----51-----52-----53-----54-----55
1411 // | | | | | | | |
1412 // | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
1413 // | | | | | | | |
1414 // 40-----41-----42-----43-----44-----45-----46-----47
1415 // | | | | | | | |
1416 // | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
1417 // | | | | | | | |
1418 // 32-----33-----34-----35-----36-----37-----38-----39
1419 // | | | | | | | |
1420 // | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
1421 // | | | | | | | |
1422 // 24-----25-----26-----27-----28-----29-----30-----31
1423 // | | | | | | | |
1424 // | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
1425 // | | | | | | | |
1426 // 16-----17-----18-----19-----20-----21-----22-----23
1427 // | | | | | | | |
1428 // | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
1429 // | | | | | | | |
1430 // 08-----09-----10-----11-----12-----13-----14-----15
1431 // | | | | | | | |
1432 // | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
1433 // | | | | | | | |
1434 // 00-----01-----02-----03-----04-----05-----06-----07
1435 //
1436 // the grid is distributed among the processes as follows:
1437 // (index == process)
1438 //
1439 // 4---4---4---4---4---4---4---4
1440 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
1441 // 4---4---4---4---4---4---4---4
1442 // | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
1443 // 3---3---3---3---4---4---4---4
1444 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
1445 // 2---2---2---2---3---3---3---3
1446 // | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
1447 // 1---1---1---1---2---2---2---2
1448 // | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
1449 // 0---0---0---0---1---1---1---1
1450 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
1451 // 0---0---0---0---0---0---0---0
1452 // | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
1453 // 0---0---0---0---0---0---0---0
1454 //
1455 // global_mask
1456 // land = 0
1457 // coast = 1
1458 // ocean = 2
1459 // +---+---+---+---+---+---+---+
1460 // | 1 | 2 | 2 | 2 | 2 | 2 | 2 |
1461 // +---+---+---+---+---+---+---+
1462 // | 1 | 2 | 2 | 2 | 2 | 2 | 2 |
1463 // +---+---+---+---+---+---+---+
1464 // | 1 | 2 | 2 | 2 | 2 | 2 | 2 |
1465 // +---+---+---+---+---+---+---+
1466 // | 1 | 2 | 2 | 2 | 2 | 2 | 2 |
1467 // +---+---+---+---+---+---+---+
1468 // | 1 | 1 | 2 | 2 | 2 | 2 | 2 |
1469 // +---+---+---+---+---+---+---+
1470 // | 0 | 1 | 1 | 2 | 2 | 2 | 2 |
1471 // +---+---+---+---+---+---+---+
1472 // | 0 | 0 | 1 | 1 | 1 | 1 | 1 |
1473 // +---+---+---+---+---+---+---+
1474 //
1475 //---------------
1476 // setup
1477 //---------------
1478
1479 int is_tgt = split_comm_size == 1;
1480 double coordinates_x[8] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
1481 double coordinates_y[8] = {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0};
1482 size_t const num_cells[2] = {7,7};
1483 size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
1484 size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{7,7}}};
1485 int global_mask[7*7] = {
1486 0,0,1,1,1,1,1,
1487 0,1,1,2,2,2,2,
1488 1,1,2,2,2,2,2,
1489 1,2,2,2,2,2,2,
1490 1,2,2,2,2,2,2,
1491 1,2,2,2,2,2,2,
1492 1,2,2,2,2,2,2};
1493 int with_halo = 0;
1494 for (size_t i = 0; i <= num_cells[0]; ++i)
1495 coordinates_x[i] *= YAC_RAD;
1496 for (size_t i = 0; i <= num_cells[1]; ++i)
1497 coordinates_y[i] *= YAC_RAD;
1498
1502 local_start[is_tgt][split_comm_rank],
1503 local_count[is_tgt][split_comm_rank], with_halo);
1504
1505 struct yac_basic_grid * grids[2] =
1508
1509 {
1510 int valid_mask_value = (is_tgt)?2:1;
1511 yac_coordinate_pointer point_coordinates =
1512 xmalloc(grid_data.num_cells * sizeof(*point_coordinates));
1513 int * mask = xmalloc(grid_data.num_cells * sizeof(*mask));
1514 for (size_t i = 0; i < grid_data.num_cells; ++i) {
1515 double * middle_point = point_coordinates[i];
1516 for (size_t k = 0; k < 3; ++k) middle_point[k] = 0.0;
1517 size_t * curr_vertices =
1518 grid_data.cell_to_vertex + grid_data.cell_to_vertex_offsets[i];
1519 size_t curr_num_vertices = grid_data.num_vertices_per_cell[i];
1520 for (size_t j = 0; j < curr_num_vertices; ++j) {
1521 double * curr_vertex_coord =
1522 grid_data.vertex_coordinates[curr_vertices[j]];
1523 for (size_t k = 0; k < 3; ++k)
1524 middle_point[k] += curr_vertex_coord[k];
1525 }
1526 normalise_vector(middle_point);
1527 mask[i] = global_mask[grid_data.cell_ids[i]] == valid_mask_value;
1528 }
1530 grids[0], YAC_LOC_CELL, point_coordinates);
1532 grids[0], YAC_LOC_CELL, mask, NULL);
1533 }
1534
1535 struct yac_dist_grid_pair * grid_pair =
1536 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
1537
1538 struct yac_interp_field src_fields[] =
1539 {{.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0}};
1540 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
1542 {.location = YAC_LOC_CELL, .coordinates_idx = 0, .masks_idx = 0};
1543
1544 struct yac_interp_grid * interp_grid =
1547
1548 enum {OVERWRITE_CONFIG_COUNT = 3};
1550 overwrite_configs[OVERWRITE_CONFIG_COUNT+1];
1551 overwrite_configs[OVERWRITE_CONFIG_COUNT] = NULL;
1552
1553 // select bounding circle containing cell with global id 8 and spread its
1554 // value to all surrounding cell of the initial target
1555 {
1556 struct yac_point_selection * src_point_selection =
1558 1.5 * YAC_RAD, 1.5 * YAC_RAD, 0.1 * YAC_RAD);
1559 struct yac_interp_spmap_config * spmap_config =
1561 1.5 * YAC_RAD,
1565 overwrite_configs[0] =
1566 yac_spmap_overwrite_config_new(src_point_selection, spmap_config);
1567 yac_interp_spmap_config_delete(spmap_config);
1568 yac_point_selection_delete(src_point_selection);
1569 }
1570
1571 // select bounding circle containing cells with global id 8, 9, 15 and
1572 // spread its value to all orthogonally adjacent cell of the initial target
1573 // (src 8 will be dealt with by the first overwrite config)
1574 {
1575 struct yac_point_selection * src_point_selection =
1577 1.5 * YAC_RAD, 1.5 * YAC_RAD, 1.1 * YAC_RAD);
1578 struct yac_interp_spmap_config * spmap_config =
1580 1.1 * YAC_RAD,
1584 overwrite_configs[1] =
1585 yac_spmap_overwrite_config_new(src_point_selection, spmap_config);
1586 yac_interp_spmap_config_delete(spmap_config);
1587 yac_point_selection_delete(src_point_selection);
1588 }
1589
1590 // select bounding circle containing cell with global id 6 and
1591 // and set its maximum search distance so low, that no matching target
1592 // point is found
1593 {
1594 struct yac_point_selection * src_point_selection =
1596 6.5 * YAC_RAD, 0.5 * YAC_RAD, 0.1 * YAC_RAD);
1597 struct yac_interp_spmap_config * spmap_config =
1600 0.1 * YAC_RAD,
1603 overwrite_configs[2] =
1604 yac_spmap_overwrite_config_new(src_point_selection, spmap_config);
1605 yac_interp_spmap_config_delete(spmap_config);
1606 yac_point_selection_delete(src_point_selection);
1607 }
1608
1609 double const fixed_value = -2.0;
1610 struct interp_method * method_stack[] =
1613 (struct yac_spmap_overwrite_config const * const *)overwrite_configs),
1614 yac_interp_method_fixed_new(fixed_value), NULL};
1615
1616 struct yac_interp_weights * weights =
1617 yac_interp_method_do_search(method_stack, interp_grid);
1618
1619 for (size_t reorder_idx = 0; reorder_idx < num_reorder_types;
1620 ++reorder_idx) {
1621
1622 struct yac_interpolation * interpolation =
1624 weights, reorder_types[reorder_idx], 1,
1625 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
1626
1627 // check generated interpolation
1628 {
1629 double * src_field = NULL;
1630 double ** src_fields = &src_field;
1631 double * tgt_field = NULL;
1632
1633 if (is_tgt) {
1634 tgt_field = xmalloc(grid_data.num_cells * sizeof(*tgt_field));
1635 for (size_t i = 0; i < grid_data.num_cells; ++i) tgt_field[i] = -1;
1636 } else {
1637 src_field = xmalloc(grid_data.num_cells * sizeof(*src_field));
1638 for (size_t i = 0; i < grid_data.num_cells; ++i)
1639 src_field[i] = (double)(grid_data.cell_ids[i]);
1640 }
1641
1642 yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
1643
1644 double ref_tgt_field[7*7];
1645 for (size_t i = 0; i < 7*7; ++i) ref_tgt_field[i] = 0.0;
1646
1647 size_t unset_tgt[] = {0,1,2,3,4,5,6, 7,8,9, 14,15, 21, 28, 35, 42};
1648 for (size_t i = 0; i < sizeof(unset_tgt)/sizeof(unset_tgt[0]); ++i)
1649 ref_tgt_field[unset_tgt[i]] = -1.0;
1650
1651 struct {
1652 yac_int global_id;
1653 size_t tgt_idx;
1654 } direct[] =
1655 {{.global_id = 2, .tgt_idx = 10}, {.global_id = 3, .tgt_idx = 10},
1656 {.global_id = 4, .tgt_idx = 11}, {.global_id = 5, .tgt_idx = 12},
1657 {.global_id = 14, .tgt_idx = 22}, {.global_id = 21, .tgt_idx = 22},
1658 {.global_id = 28, .tgt_idx = 29}, {.global_id = 35, .tgt_idx = 36},
1659 {.global_id = 42, .tgt_idx = 43}};
1660 for (size_t i = 0; i < sizeof(direct)/sizeof(direct[0]); ++i)
1661 ref_tgt_field[direct[i].tgt_idx] += (double)(direct[i].global_id);
1662
1663 struct {
1664 yac_int global_id;
1665 size_t * tgt_idx;
1666 size_t tgt_count;
1667 } wsum[] =
1668 {{.global_id = 8,
1669 .tgt_idx = (size_t[]){10, 16,17, 22,23,24},
1670 .tgt_count = 6},
1671 {.global_id = 9,
1672 .tgt_idx = (size_t[]){10,11, 17},
1673 .tgt_count = 3},
1674 {.global_id = 15,
1675 .tgt_idx = (size_t[]){16,17, 23},
1676 .tgt_count = 3}};
1677 for (size_t i = 0; i < sizeof(wsum)/sizeof(wsum[0]); ++i)
1678 for (size_t j = 0; j < wsum[i].tgt_count; ++j)
1679 ref_tgt_field[wsum[i].tgt_idx[j]] +=
1680 (double)(wsum[i].global_id) / (double)(wsum[i].tgt_count);
1681
1682 size_t fixed_tgt[] =
1683 {13, 18,19,20, 25,26,27, 30,31,32,33,34,
1684 37,38,39,40,41, 44,45,46,47,48};
1685 for (size_t i = 0; i < sizeof(fixed_tgt)/sizeof(fixed_tgt[0]); ++i)
1686 ref_tgt_field[fixed_tgt[i]] = fixed_value;
1687
1688 if (is_tgt)
1689 for (size_t i = 0; i < grid_data.num_cells; ++i)
1690 if (fabs(
1691 ref_tgt_field[grid_data.cell_ids[i]] -
1692 tgt_field[i]) > 1e-9) PUT_ERR("wrong interpolation result");
1693
1694 free(tgt_field);
1695 free(src_field);
1696 }
1697
1698 yac_interpolation_delete(interpolation);
1699 }
1700
1702 yac_interp_method_delete(method_stack);
1703 for (size_t i = 0; i < OVERWRITE_CONFIG_COUNT; ++i)
1704 yac_spmap_overwrite_config_delete(overwrite_configs[i]);
1705 yac_interp_grid_delete(interp_grid);
1706 yac_dist_grid_pair_delete(grid_pair);
1709 }
1710
1711 MPI_Comm_free(&split_comm);
1712
1713 xt_finalize();
1714
1715 MPI_Finalize();
1716
1717 return TEST_EXIT_CODE;
1718}
1719
1720static double utest_get_sq_sphere_radius(
1721 struct yac_spmap_cell_area_config const * cell_area_config) {
1722
1723 double sq_sphere_radius;
1724 if (
1725 yac_spmap_cell_area_config_get_type(cell_area_config) ==
1727 double sphere_radius =
1729 sq_sphere_radius = sphere_radius * sphere_radius;
1730 } else sq_sphere_radius = 1.0;
1731 return sq_sphere_radius;
1732}
1733
1734static double utest_compute_scale(
1735 enum yac_interp_spmap_scale_type scale_type,
1736 double src_cell_area, double tgt_cell_area) {
1737
1738 double scale = 1.0;
1739
1740 if ((scale_type == YAC_INTERP_SPMAP_SRCAREA) ||
1741 (scale_type == YAC_INTERP_SPMAP_FRACAREA))
1742 scale *= src_cell_area;
1743 if ((scale_type == YAC_INTERP_SPMAP_INVTGTAREA) ||
1744 (scale_type == YAC_INTERP_SPMAP_FRACAREA))
1745 scale /= tgt_cell_area;
1746
1747 return scale;
1748}
1749
1750static void utest_write_area_file(
1751 char const * filename, double * cell_areas, size_t dim_x, size_t dim_y) {
1752
1753 // create file
1754 int ncid;
1755 yac_nc_create(filename, NC_CLOBBER, &ncid);
1756
1757 int dim_x_id;
1758 int dim_y_id;
1759 int dim_xy_id;
1760
1761 // define dimensions
1762 YAC_HANDLE_ERROR(nc_def_dim(ncid, "x", dim_x, &dim_x_id));
1763 YAC_HANDLE_ERROR(nc_def_dim(ncid, "y", dim_y, &dim_y_id));
1764 YAC_HANDLE_ERROR(nc_def_dim(ncid, "xy", dim_x * dim_y, &dim_xy_id));
1765
1766 char const * area_2d_varname = "area_2d";
1767 char const * area_1d_varname = "area_1d";
1768
1769 int area_2d_dim_ids[2] = {dim_x_id, dim_y_id};
1770 int area_1d_dim_ids[1] = {dim_xy_id};
1771
1772 int area_2d_varid;
1773 int area_1d_varid;
1774
1775 // define variable
1777 nc_def_var(
1778 ncid, area_2d_varname, NC_DOUBLE, 2, area_2d_dim_ids, &area_2d_varid));
1780 nc_def_var(
1781 ncid, area_1d_varname, NC_DOUBLE, 1, area_1d_dim_ids, &area_1d_varid));
1782
1783 // end definition
1784 YAC_HANDLE_ERROR(nc_enddef(ncid));
1785
1786 // write grid data
1787 YAC_HANDLE_ERROR(nc_put_var_double(ncid, area_2d_varid, cell_areas));
1788 YAC_HANDLE_ERROR(nc_put_var_double(ncid, area_1d_varid, cell_areas));
1789
1790 // close file
1791 YAC_HANDLE_ERROR(nc_close(ncid));
1792}
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)
#define YAC_RAD
static void normalise_vector(double v[])
Definition geometry.h:635
static double get_vector_angle(double const a[3], double const b[3])
Definition geometry.h:340
void yac_interp_grid_delete(struct yac_interp_grid *interp_grid)
struct yac_interp_grid * yac_interp_grid_new(struct yac_dist_grid_pair *grid_pair, char const *src_grid_name, char const *tgt_grid_name, size_t num_src_fields, struct yac_interp_field const *src_fields, struct yac_interp_field const tgt_field)
Definition interp_grid.c:31
void yac_interp_method_delete(struct interp_method **method)
struct yac_interp_weights * yac_interp_method_do_search(struct interp_method **method, struct yac_interp_grid *interp_grid)
struct interp_method * yac_interp_method_fixed_new(double value)
struct yac_spmap_cell_area_config const * yac_spmap_scale_config_get_tgt_cell_area_config(struct yac_spmap_scale_config const *scale_config)
double yac_spmap_cell_area_config_get_sphere_radius(struct yac_spmap_cell_area_config const *cell_area_config)
struct yac_spmap_cell_area_config * yac_spmap_cell_area_config_file_new(char const *filename, char const *varname, yac_int min_global_id)
struct yac_spmap_overwrite_config * yac_spmap_overwrite_config_new(struct yac_point_selection const *src_point_selection, struct yac_interp_spmap_config const *config)
enum yac_interp_spmap_cell_area_provider yac_spmap_cell_area_config_get_type(struct yac_spmap_cell_area_config const *cell_area_config)
struct yac_spmap_cell_area_config const * yac_spmap_scale_config_get_src_cell_area_config(struct yac_spmap_scale_config const *scale_config)
struct interp_method * yac_interp_method_spmap_new(struct yac_interp_spmap_config const *default_config, struct yac_spmap_overwrite_config const *const *overwrite_configs)
struct yac_interp_spmap_config * yac_interp_spmap_config_new(double spread_distance, double max_search_distance, enum yac_interp_spmap_weight_type weight_type, struct yac_spmap_scale_config const *scale_config)
void yac_spmap_scale_config_delete(struct yac_spmap_scale_config *scale_config)
void yac_interp_spmap_config_delete(struct yac_interp_spmap_config *config)
void yac_spmap_cell_area_config_delete(struct yac_spmap_cell_area_config *cell_area_config)
struct yac_spmap_cell_area_config * yac_spmap_cell_area_config_yac_new(double sphere_radius)
void yac_spmap_overwrite_config_delete(struct yac_spmap_overwrite_config *overwrite_config)
struct yac_spmap_scale_config * yac_spmap_scale_config_new(enum yac_interp_spmap_scale_type scale_type, struct yac_spmap_cell_area_config const *source_cell_area_config, struct yac_spmap_cell_area_config const *target_cell_area_config)
#define YAC_INTERP_SPMAP_MAX_SEARCH_DISTANCE_DEFAULT
yac_interp_spmap_scale_type
@ YAC_INTERP_SPMAP_NONE
weights are not scaled
@ YAC_INTERP_SPMAP_INVTGTAREA
@ YAC_INTERP_SPMAP_SRCAREA
@ YAC_INTERP_SPMAP_FRACAREA
#define YAC_INTERP_SPMAP_SCALE_CONFIG_DEFAULT
#define YAC_INTERP_SPMAP_WEIGHTED_DEFAULT
#define YAC_INTERP_SPMAP_SPREAD_DISTANCE_DEFAULT
@ YAC_INTERP_SPMAP_CELL_AREA_YAC
yac_interp_spmap_weight_type
@ YAC_INTERP_SPMAP_AVG
@ YAC_INTERP_SPMAP_DIST
#define YAC_INTERP_SPMAP_OVERWRITE_DEFAULT
#define YAC_INTERP_SPMAP_DEFAULT_CONFIG
struct yac_interpolation * yac_interp_weights_get_interpolation(struct yac_interp_weights *weights, enum yac_interp_weights_reorder_type reorder, size_t collection_size, double frac_mask_fallback_value, double scaling_factor, double scaling_summand, char const *yaxt_exchanger_name, int is_source, int is_target)
void yac_interp_weights_delete(struct yac_interp_weights *weights)
yac_interp_weights_reorder_type
@ YAC_MAPPING_ON_TGT
weights will be applied at target processes
@ YAC_MAPPING_ON_SRC
weights will be applied at source processes
void yac_interpolation_execute(struct yac_interpolation *interp, double ***src_fields, double **tgt_field)
Execute interpolation synchronously and write results to the target field.
void yac_interpolation_delete(struct yac_interpolation *interp)
Free an interpolation object and release all resources.
double const YAC_FRAC_MASK_NO_VALUE
void yac_nc_create(const char *path, int cmode, int *ncidp)
Definition io_utils.c:367
@ YAC_LOC_CELL
Definition location.h:14
Definition __init__.py:1
struct yac_point_selection * yac_point_selection_bnd_circle_new(double center_lon, double center_lat, double inc_angle)
void yac_point_selection_delete(struct yac_point_selection *point_select)
#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
struct yac_spmap_scale_config * scale_config
struct yac_spmap_cell_area_config * tgt
void clear_yac_io_env()
static MPI_Comm split_comm
static char const * grid_names[2]
enum yac_interp_spmap_scale_type scale_types[]
enum yac_interp_weights_reorder_type reorder_types[]
enum yac_interp_spmap_weight_type weight_types[]
double coordinates_x[]
size_t num_cells[2]
double coordinates_y[]
static int mask[16]
#define TEST_EXIT_CODE
Definition tests.h:14
#define PUT_ERR(string)
Definition tests.h:10
#define YAC_HANDLE_ERROR(exp)
Definition toy_output.c:13
struct yac_basic_grid ** grids
Definition yac.c:152
YAC_INT yac_int
Definition yac_types.h:15
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19