YAC 3.12.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_interp_method_ncc_parallel.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#include <stdlib.h>
6
7#include "tests.h"
8#include "test_common.h"
9#include "geometry.h"
10#include "interpolation/methods/interp_method.h"
11#include "interpolation/methods/interp_method_ncc.h"
12#include "interpolation/methods/interp_method_fixed.h"
13#include "dist_grid_utils.h"
14#include "clipping.h"
15#include "yac_mpi.h"
16
17#include <mpi.h>
18#include <yaxt.h>
19#include <netcdf.h>
20
28size_t num_reorder_types = sizeof(reorder_types) / sizeof(reorder_types[0]);
29static char const * grid_names[2] = {"src_grid", "tgt_grid"};
30
31static double utest_compute_ref_tgt_value(
32 size_t tgt_corner, int * tgt_global_core_mask, int * tgt_global_field_mask,
33 int * src_global_core_mask, int * src_global_field_maske, int with_src_mask,
34 int with_tgt_mask, enum yac_interp_ncc_weight_type weight_type,
35 int partial_coverage);
36
37int main(void) {
38
39 MPI_Init(NULL, NULL);
40
41 xt_initialize(MPI_COMM_WORLD);
42
43 int comm_rank, comm_size;
44 MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
45 MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
46 MPI_Barrier(MPI_COMM_WORLD);
47
48 if (comm_size != 3) {
49 PUT_ERR("ERROR: wrong number of processes");
50 xt_finalize();
51 MPI_Finalize();
52 return TEST_EXIT_CODE;
53 }
54
55 MPI_Comm split_comm;
56 MPI_Comm_split(
57 MPI_COMM_WORLD, comm_rank < 2, 0, &split_comm);
58
59 int split_comm_rank, split_comm_size;
60 MPI_Comm_rank(split_comm, &split_comm_rank);
61 MPI_Comm_size(split_comm, &split_comm_size);
62
63 {// test with various configuration options
64
65 // The global source grid has 5x4 cells:
66 //
67 // 24--44--25--45--26--46--27--47--28--48--29
68 // | | | | | |
69 // 34 15 36 16 38 17 40 18 42 19 43
70 // | | | | | |
71 // 18--33--19--35--20--37--21--39--22--41--23
72 // | | | | | |
73 // 23 10 25 11 27 12 29 13 31 14 32
74 // | | | | | |
75 // 12--22--13--24--14--26--15--28--16--30--17
76 // | | | | | |
77 // 12 05 14 06 16 07 18 08 20 09 21
78 // | | | | | |
79 // 06--11--07--13--08--15--09--17--10--19--11
80 // | | | | | |
81 // 01 00 03 01 05 02 07 03 09 04 10
82 // | | | | | |
83 // 00--01--01--02--02--04--03--06--04--08--05
84 //
85 //---------------
86 // setup
87 //---------------
88
89 int is_tgt = split_comm_size == 1;
90 double coordinates_x[2][10] =
91 {{0,1,2,3,4,5}, {0.25,0.75,1.25,1.75,2.25,2.75,3.25,3.75,4.25,4.75}};
92 double coordinates_y[2][8] =
93 {{-2,-1,0,1,2}, {-1.75,-1.25,-0.75,-0.25,0.25,0.75,1.25,1.75}};
94 double cell_coordiantes_x[5 * 4] =
95 {0.5,1.5,2.5,3.5,4.5,
96 0.5,1.5,2.5,3.5,4.5,
97 0.5,1.5,2.5,3.5,4.5,
98 0.5,1.5,2.5,3.5,4.5};
99 double cell_coordiantes_y[5 * 4] =
100 {-1.5,-1.5,-1.5,-1.5,-1.5,
101 -0.5,-0.5,-0.5,-0.5,-0.5,
102 0.5,0.5,0.5,0.5,0.5,
103 1.5,1.5,1.5,1.5,1.5};
104 size_t const num_cells[2][2] = {{5,4}, {9,7}};
105 size_t local_start[2][2][2] = {{{0,0},{3,0}}, {{0,0}}};
106 size_t local_count[2][2][2] = {{{3,4},{2,4}}, {{9,7}}};
107 int global_core_mask[2][10*8] =
108 {{1,0,0,0,0,
109 1,1,0,0,0,
110 1,1,1,1,0,
111 1,1,1,1,0},
112 {1,1,1,0,0,0,0,0,0,0,
113 1,1,1,0,0,0,0,0,0,0,
114 1,1,1,1,1,1,0,0,0,0,
115 1,1,1,1,1,1,0,0,0,0,
116 1,1,1,1,1,1,1,1,1,0,
117 1,1,1,1,1,1,1,1,1,0,
118 1,1,1,1,1,1,1,1,1,0,
119 1,1,1,1,1,1,1,1,1,0}};
120 int global_field_mask[2][10*8] =
121 {{0,1,1,1,1,
122 0,1,1,1,1,
123 0,1,1,1,1,
124 0,1,1,1,1},
125 {1,1,1,1,1,1,1,1,1,1,
126 1,1,1,1,1,1,1,1,1,1,
127 1,1,1,1,1,1,1,1,1,1,
128 1,1,1,1,1,1,1,1,1,1,
129 1,1,1,1,1,1,1,1,1,1,
130 1,1,1,1,1,1,1,1,1,1,
131 1,1,1,1,1,1,1,1,1,1,
132 0,0,0,0,0,0,0,0,0,0}};
133 int with_halo = 0;
134 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
135 coordinates_x[is_tgt][i] *= YAC_RAD;
136 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
137 coordinates_y[is_tgt][i] *= YAC_RAD;
138
139 struct yac_basic_grid_data grid_data =
140 yac_generate_basic_grid_data_reg2d(
141 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
142 local_start[is_tgt][split_comm_rank],
143 local_count[is_tgt][split_comm_rank], with_halo);
144
145 if (is_tgt) {
146 for (size_t i = 0; i < grid_data.num_vertices; ++i)
147 grid_data.core_vertex_mask[i] =
148 global_core_mask[is_tgt][grid_data.vertex_ids[i]];
149 free(grid_data.core_edge_mask);
150 free(grid_data.core_cell_mask);
151 grid_data.core_edge_mask = NULL;
152 grid_data.core_cell_mask = NULL;
153 } else {
154 for (size_t i = 0; i < grid_data.num_cells; ++i)
155 grid_data.core_cell_mask[i] =
156 global_core_mask[is_tgt][grid_data.cell_ids[i]];
157 }
158
159 struct yac_basic_grid * grids[2] =
160 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
161 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
162
163 int field_mask[10*8];
164 double field_coords[10*8][3];
165 if (is_tgt) {
166 for (size_t i = 0; i < grid_data.num_vertices; ++i)
167 field_mask[i] = global_field_mask[is_tgt][grid_data.vertex_ids[i]];
168 yac_basic_grid_add_mask(
169 grids[0], YAC_LOC_CORNER, field_mask, grid_data.num_vertices, NULL);
170 } else {
171 for (size_t i = 0; i < grid_data.num_cells; ++i) {
172 field_mask[i] = global_field_mask[is_tgt][grid_data.cell_ids[i]];
173 LLtoXYZ_deg(
174 cell_coordiantes_x[grid_data.cell_ids[i]],
175 cell_coordiantes_y[grid_data.cell_ids[i]], field_coords[i]);
176 }
177 yac_basic_grid_add_mask(
178 grids[0], YAC_LOC_CELL, field_mask, grid_data.num_cells, NULL);
179 yac_basic_grid_add_coordinates(
180 grids[0], YAC_LOC_CELL, field_coords, grid_data.num_cells);
181 }
182
183 struct yac_dist_grid_pair * grid_pair =
184 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
185
186 for (int with_src_mask = 0; with_src_mask < 2; ++with_src_mask) {
187 for (int with_tgt_mask = 0; with_tgt_mask < 2; ++with_tgt_mask) {
188
189 struct yac_interp_field src_fields[] =
190 {{.location = YAC_LOC_CELL,
191 .coordinates_idx = 0,
192 .masks_idx = (with_src_mask)?0:SIZE_MAX}};
193 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
194 struct yac_interp_field tgt_field =
195 {.location = YAC_LOC_CORNER,
196 .coordinates_idx = SIZE_MAX,
197 .masks_idx = (with_tgt_mask)?0:SIZE_MAX};
198
199 struct yac_interp_grid * interp_grid =
200 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
201 num_src_fields, src_fields, tgt_field);
202
203 enum yac_interp_ncc_weight_type weight_types[] =
204 {YAC_INTERP_NCC_AVG,
205 YAC_INTERP_NCC_DIST};
206 size_t const num_weight_types =
207 sizeof(weight_types) / sizeof(weight_types[0]);
208
209 for (size_t weight_types_idx = 0; weight_types_idx < num_weight_types;
210 ++weight_types_idx) {
211 for (int partial_coverage = 0; partial_coverage < 2;
212 ++partial_coverage) {
213
214 struct interp_method * method_stack[] =
215 {yac_interp_method_ncc_new(
216 weight_types[weight_types_idx], partial_coverage),
217 yac_interp_method_fixed_new(-1.0), NULL};
218
219 struct yac_interp_weights * weights =
220 yac_interp_method_do_search(method_stack, interp_grid);
221
222 for (size_t i = 0; i < num_reorder_types; ++i) {
223
224 struct yac_interpolation * interpolation =
225 yac_interp_weights_get_interpolation(
226 weights, reorder_types[i], 1,
227 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
228
229 // check generated interpolation
230 {
231 double src_field_data[5*4];
232 double * src_field = src_field_data;
233 double ** src_fields = &src_field;
234 double tgt_field_data[10*8] =
235 {-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,
236 -2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,
237 -2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,
238 -2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,
239 -2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,
240 -2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,
241 -2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,
242 -2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,-2.0};
243 double * tgt_field = tgt_field_data;
244
245 if (!is_tgt)
246 for (size_t i = 0; i < grid_data.num_cells; ++i)
247 src_field[i] = (double)(grid_data.cell_ids[i] + 1);
248
249 yac_interpolation_execute(
250 interpolation, &src_fields, &tgt_field);
251
252 if (is_tgt) {
253
254 for (size_t i = 0; i < grid_data.num_vertices; ++i) {
255
256 double ref_tgt_value =
257 utest_compute_ref_tgt_value(
258 i, global_core_mask[1], global_field_mask[1],
259 global_core_mask[0], global_field_mask[0],
260 with_src_mask, with_tgt_mask,
261 weight_types[weight_types_idx], partial_coverage);
262
263 if (fabs(ref_tgt_value - tgt_field[i]) > 1e-3)
264 PUT_ERR("wrong interpolation result");
265 }
266 }
267 }
268
269 yac_interpolation_delete(interpolation);
270 } // reorder type
271 yac_interp_weights_delete(weights);
272 yac_interp_method_delete(method_stack);
273 } // partial coverage
274 } // weight type
275 yac_interp_grid_delete(interp_grid);
276 } // with tgt mask
277 } // with src mask
278 yac_dist_grid_pair_delete(grid_pair);
279 yac_basic_grid_delete(grids[1]);
280 yac_basic_grid_delete(grids[0]);
281 }
282
283 {// small test
284
285 // The global source grid has 5x4 cells:
286 //
287 //---------------
288 // setup
289 //---------------
290
291 int is_tgt = split_comm_size == 1;
292 double coordinates_x[2][4] =
293 {{-2,-1,1,2}, {-1.0,0.0,1.0}};
294 double coordinates_y[2][4] =
295 {{-2,-1,1,2}, {-1.0,0.0,1.0}};
296 double field_coordiantes_x[2][3*3] =
297 {{-1.5,0.0,1.5,
298 -1.5,0.0,1.5,
299 -1.5,0.0,1.5},
300 {-1.0,0.0,1.0,
301 -1.0,0.0,1.0,
302 -1.0,0.0,1.0}};
303 double field_coordiantes_y[2][3*3] =
304 {{-1.5,-1.5,-1.5,
305 0.0,0.0,0.0,
306 1.5,1.5,1.5},
307 {-1.0,-1.0,-1.0,
308 0.0,0.0,0.0,
309 1.0,1.0,1.0}};
310 int field_mask[3*3] = {0,1,1, 1,1,1, 1,1,0};
311 size_t const num_cells[2][2] = {{3,3}, {2,2}};
312 size_t local_start[2][2][2] = {{{0,0},{0,1}}, {{0,0}}};
313 size_t local_count[2][2][2] = {{{3,2},{3,2}}, {{2,2}}};
314 int with_halo = 0;
315 for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
316 coordinates_x[is_tgt][i] *= YAC_RAD;
317 for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
318 coordinates_y[is_tgt][i] *= YAC_RAD;
319
320 struct yac_basic_grid_data grid_data =
321 yac_generate_basic_grid_data_reg2d(
322 coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
323 local_start[is_tgt][split_comm_rank],
324 local_count[is_tgt][split_comm_rank], with_halo);
325
326 struct yac_basic_grid * grids[2] =
327 {yac_basic_grid_new(grid_names[is_tgt], grid_data),
328 yac_basic_grid_empty_new(grid_names[is_tgt^1])};
329
330 double field_coords[3*3][3];
331 {
332 const yac_int * ids = (is_tgt)?grid_data.vertex_ids:grid_data.cell_ids;
333 size_t count = (is_tgt)?grid_data.num_vertices:grid_data.num_cells;
334 enum yac_location field_location[2] = {YAC_LOC_CELL, YAC_LOC_CORNER};
335 int src_field_mask[3*3];
336 for (size_t i = 0; i < count; ++i) {
337 LLtoXYZ_deg(
338 field_coordiantes_x[is_tgt][ids[i]],
339 field_coordiantes_y[is_tgt][ids[i]], field_coords[i]);
340 if (!is_tgt) src_field_mask[i] = field_mask[ids[i]];
341 }
342 yac_basic_grid_add_coordinates(
343 grids[0], field_location[is_tgt], field_coords, count);
344 if (!is_tgt)
345 yac_basic_grid_add_mask(
346 grids[0], field_location[is_tgt], src_field_mask, count, NULL);
347 }
348
349 struct yac_dist_grid_pair * grid_pair =
350 yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
351
352 for (int with_src_mask = 0; with_src_mask < 2; ++with_src_mask) {
353
354 struct yac_interp_field src_fields[] =
355 {{.location = YAC_LOC_CELL,
356 .coordinates_idx = 0,
357 .masks_idx = with_src_mask?0:SIZE_MAX}};
358 size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
359 struct yac_interp_field tgt_field =
360 {.location = YAC_LOC_CORNER,
361 .coordinates_idx = 0,
362 .masks_idx = SIZE_MAX};
363
364 struct yac_interp_grid * interp_grid =
365 yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
366 num_src_fields, src_fields, tgt_field);
367
368 for (int partial_coverage = 0; partial_coverage < 2; ++partial_coverage) {
369
370 struct interp_method * method_stack[] =
371 {yac_interp_method_ncc_new(
372 YAC_INTERP_NCC_DIST, partial_coverage),
373 yac_interp_method_fixed_new(-1.0), NULL};
374
375 struct yac_interp_weights * weights =
376 yac_interp_method_do_search(method_stack, interp_grid);
377
378 struct yac_interpolation * interpolation =
379 yac_interp_weights_get_interpolation(
380 weights, YAC_MAPPING_ON_SRC, 1,
381 YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0, NULL, 1, 1);
382
383 // check generated interpolation
384 {
385 double src_field_data[3*3];
386 double * src_field = src_field_data;
387 double ** src_fields = &src_field;
388 double tgt_field_data[3*3] =
389 {-2.0,-2.0,-2.0,
390 -2.0,-2.0,-2.0,
391 -2.0,-2.0,-2.0};
392 double * tgt_field = tgt_field_data;
393
394 if (!is_tgt)
395 for (size_t i = 0; i < grid_data.num_cells; ++i)
396 src_field[i] = (double)(grid_data.cell_ids[i] + 1);
397
398 yac_interpolation_execute(
399 interpolation, &src_fields, &tgt_field);
400
401 if (is_tgt) {
402
403 yac_int tgt_to_src[2][2][9][4] =
404 {{{{0,1,3,4},{0,1,3,4},{1,2,4,5},
405 {0,1,3,4},{4},{1,2,4,5},
406 {3,4,6,7},{3,4,6,7},{4,5,7,8}},
407 {{0,1,3,4},{1,2,4,5},{1,2,4,5},
408 {3,4,6,7},{4},{4,5,7,8},
409 {3,4,6,7},{4,5,7,8},{4,5,7,8}}},
410 {{{1,3,4},{1,3,4},{1,2,4,5},
411 {1,3,4},{4},{1,2,4,5},
412 {3,4,6,7},{3,4,6,7},{4,5,7}},
413 {{1,3,4},{1,2,4,5},{1,2,4,5},
414 {3,4,6,7},{4},{4,5,7},
415 {3,4,6,7},{4,5,7},{4,5,7}}}};
416 int num_src_per_tgt[2][2][2][9] =
417 {{{{4,4,4, 4,1,4, 4,4,4},{4,4,4, 4,1,4, 4,4,4}},
418 {{0,0,4, 0,1,4, 4,4,0},{0,4,4, 4,1,0, 4,0,0}}},
419 {{{4,4,4, 4,1,4, 4,4,4},{4,4,4, 4,1,4, 4,4,4}},
420 {{3,3,4, 3,1,4, 4,4,3},{3,4,4, 4,1,3, 4,3,3}}}};
421
422 for (size_t i = 0; i < grid_data.num_vertices; ++i) {
423
424 yac_int tgt_id = grid_data.vertex_ids[i];
425 double ref_tgt_value[2] = {0.0, 0.0};
426
427 double tgt_coord[3];
428 LLtoXYZ_deg(
429 field_coordiantes_x[1][tgt_id], field_coordiantes_y[1][tgt_id],
430 tgt_coord);
431
432 for (int j = 0; j < 2; ++j) {
433
434 double inv_distance_sum = 0.0;
435 int num_src =
436 num_src_per_tgt[partial_coverage][with_src_mask][j][tgt_id];
437
438 if (num_src > 0) {
439 for (int k = 0; k < num_src; ++k) {
440 yac_int src_id = tgt_to_src[with_src_mask][j][tgt_id][k];
441 double src_coord[3];
442 LLtoXYZ_deg(
443 field_coordiantes_x[0][src_id],
444 field_coordiantes_y[0][src_id], src_coord);
445 double angle = get_vector_angle(src_coord, tgt_coord);
446 if (angle < yac_angle_tol) {
447 ref_tgt_value[j] = (double)(src_id + 1);
448 inv_distance_sum = 1.0;
449 break;
450 }
451 double inv_distance =
452 1.0 / get_vector_angle(src_coord, tgt_coord);
453 inv_distance_sum += inv_distance;
454 ref_tgt_value[j] += (double)(src_id + 1) * inv_distance;
455 }
456
457 ref_tgt_value[j] /= inv_distance_sum;
458 } else {
459 ref_tgt_value[j] = -1.0;
460 }
461 }
462
463 if ((fabs(ref_tgt_value[0] - tgt_field[i]) > 1e-3) &&
464 (fabs(ref_tgt_value[1] - tgt_field[i]) > 1e-3))
465 PUT_ERR("wrong interpolation result");
466 }
467 }
468 }
469
470 yac_interpolation_delete(interpolation);
471 yac_interp_weights_delete(weights);
472 yac_interp_method_delete(method_stack);
473 } // partial_coverage
474 yac_interp_grid_delete(interp_grid);
475 } // with_src_mask
476 yac_dist_grid_pair_delete(grid_pair);
477 yac_basic_grid_delete(grids[1]);
478 yac_basic_grid_delete(grids[0]);
479 }
480
481 MPI_Comm_free(&split_comm);
482
483 xt_finalize();
484
485 MPI_Finalize();
486
487 return TEST_EXIT_CODE;
488}
489
490static double utest_compute_ref_tgt_value(
491 size_t tgt_corner, int * tgt_global_core_mask, int * tgt_global_field_mask,
492 int * src_global_core_mask, int * src_global_field_mask, int with_src_mask,
493 int with_tgt_mask, enum yac_interp_ncc_weight_type weight_type,
494 int partial_coverage) {
495
496 // if target corner is masked
497 if (!tgt_global_core_mask[tgt_corner] ||
498 (with_tgt_mask && !tgt_global_field_mask[tgt_corner])) return -2.0;
499
500 // determine matching source cell
501 int source_grid_x_idx = ((int)tgt_corner % 10) / 2;
502 int source_grid_y_idx = ((int)tgt_corner / 10) / 2;
503 int source_cell = source_grid_x_idx + source_grid_y_idx * 5;
504
505 // check source cell core mask
506 if (!src_global_core_mask[source_cell]) return -1.0;
507
508 // determine closest source cell corner
509 int is_right_corner = (int)tgt_corner & 1;
510 int is_upper_corner = ((int)tgt_corner / 10) & 1;
511 int source_corner =
512 source_grid_x_idx + is_right_corner +
513 6 * (source_grid_y_idx + is_upper_corner);
514
515 // determine source cells surrounding the source vertex
516 int source_cells[4];
517 source_cells[0] = (5 * (source_corner / 6) + source_corner % 6) - 5 - 1;
518 source_cells[1] = (5 * (source_corner / 6) + source_corner % 6) - 5;
519 source_cells[2] = (5 * (source_corner / 6) + source_corner % 6) - 1;
520 source_cells[3] = (5 * (source_corner / 6) + source_corner % 6);
521
522 // special handling for grid edge corners
523 if (source_corner <= 5) {
524 source_cells[0] = -1;
525 source_cells[1] = -1;
526 }
527 if (source_corner >= 24) {
528 source_cells[2] = -1;
529 source_cells[3] = -1;
530 }
531 if ((source_corner % 6) == 0) {
532 source_cells[0] = -1;
533 source_cells[2] = -1;
534 }
535 if (((source_corner + 1) % 6) == 0) {
536 source_cells[3] = -1;
537 source_cells[1] = -1;
538 }
539
540 // apply core mask to source cells
541 for (int i = 0; i < 4; ++i)
542 if ((source_cells[i] >= 0) &&
543 (!src_global_core_mask[source_cells[i]]))
544 source_cells[i] = -1;
545
546 // apply field mask to source cells
547 if (with_src_mask) {
548 for (int i = 0; i < 4; ++i) {
549 if ((source_cells[i] >= 0) &&
550 (!src_global_field_mask[source_cells[i]])) {
551 if (!partial_coverage) return -1.0;
552 source_cells[i] = -1;
553 }
554 }
555 }
556
557 double const ref_distances[2][2] =
558 {{sqrt(0.25*0.25 + 0.25*0.25),
559 sqrt(0.75*0.75 + 0.25*0.25)},
560 {sqrt(0.25*0.25 + 0.75*0.75),
561 sqrt(0.75*0.75 + 0.75*0.75)}};
562
563 double distances[4];
564 int num_src_cells = 0;
565 for (int i = 0, k = 0; i < 2; ++i) {
566 for (int j = 0; j < 2; ++j, ++k) {
567 if (source_cells[k] >= 0) {
568 distances[num_src_cells] =
569 ref_distances[i^is_upper_corner^1][j^is_right_corner^1];
570 source_cells[num_src_cells] = source_cells[k];
571 ++num_src_cells;
572 }
573 }
574 }
575
576 if (num_src_cells == 0) return -1.0;
577
578 double tgt_value = 0.0;
579 switch(weight_type) {
580 default:
581 case (YAC_INTERP_NCC_AVG): {
582 double weight = 1.0 / (double)num_src_cells;
583 for (int i = 0; i < num_src_cells; ++i)
584 tgt_value += (double)(source_cells[i] + 1) * weight;
585 break;
586 }
587 case (YAC_INTERP_NCC_DIST): {
588 double inv_distance_sum = 0.0;
589 for (int i = 0; i < num_src_cells; ++i)
590 inv_distance_sum += 1.0 / distances[i];
591 for (int i = 0; i < num_src_cells; ++i)
592 tgt_value +=
593 (double)(source_cells[i] + 1) / (distances[i] * inv_distance_sum);
594 break;
595 }
596 };
597 return tgt_value;
598}
yac_basic_grid_new
struct yac_basic_grid * yac_basic_grid_new(char const *name, struct yac_basic_grid_data grid_data)
Definition basic_grid.c:50
yac_basic_grid_add_mask
size_t yac_basic_grid_add_mask(struct yac_basic_grid *grid, enum yac_location location, int const *mask, size_t count, char const *mask_name)
Definition basic_grid.c:284
yac_basic_grid_add_coordinates
size_t yac_basic_grid_add_coordinates(struct yac_basic_grid *grid, enum yac_location location, yac_coordinate_pointer coordinates, size_t count)
Definition basic_grid.c:232
yac_basic_grid_empty_new
struct yac_basic_grid * yac_basic_grid_empty_new(char const *name)
Definition basic_grid.c:63
yac_basic_grid_delete
void yac_basic_grid_delete(struct yac_basic_grid *grid)
Definition basic_grid.c:70
yac_dist_grid_pair_delete
void yac_dist_grid_pair_delete(struct yac_dist_grid_pair *grid_pair)
Definition dist_grid.c:2313
yac_dist_grid_pair_new
struct yac_dist_grid_pair * yac_dist_grid_pair_new(struct yac_basic_grid *grid_a, struct yac_basic_grid *grid_b, MPI_Comm comm)
Definition dist_grid.c:2061
yac_generate_basic_grid_data_reg2d
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)
Definition dist_grid_utils.c:200
dist_grid_utils.h
YAC_RAD
#define YAC_RAD
Definition dynamic_config.c:12
LLtoXYZ_deg
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:269
yac_angle_tol
#define yac_angle_tol
Definition geometry.h:26
get_vector_angle
static double get_vector_angle(double const a[3], double const b[3])
Definition geometry.h:340
yac_interp_grid_delete
void yac_interp_grid_delete(struct yac_interp_grid *interp_grid)
Definition interp_grid.c:575
yac_interp_grid_new
struct yac_interp_grid * yac_interp_grid_new(struct yac_dist_grid_pair *grid_pair, char const *src_grid_name, char const *tgt_grid_name, size_t num_src_fields, struct yac_interp_field const *src_fields, struct yac_interp_field const tgt_field)
Definition interp_grid.c:31
yac_interp_method_delete
void yac_interp_method_delete(struct interp_method **method)
Definition interp_method.c:54
yac_interp_method_do_search
struct yac_interp_weights * yac_interp_method_do_search(struct interp_method **method, struct yac_interp_grid *interp_grid)
Definition interp_method.c:15
interp_method.h
yac_interp_method_fixed_new
struct interp_method * yac_interp_method_fixed_new(double value)
Definition interp_method_fixed.c:52
interp_method_fixed.h
yac_interp_method_ncc_new
struct interp_method * yac_interp_method_ncc_new(enum yac_interp_ncc_weight_type weight_type, int partial_coverage)
Definition interp_method_ncc.c:426
interp_method_ncc.h
yac_interp_ncc_weight_type
yac_interp_ncc_weight_type
Definition interp_method_ncc.h:12
YAC_INTERP_NCC_DIST
@ YAC_INTERP_NCC_DIST
distance weighted average of n source points
Definition interp_method_ncc.h:14
YAC_INTERP_NCC_AVG
@ YAC_INTERP_NCC_AVG
average of n source points
Definition interp_method_ncc.h:13
yac_interp_weights_get_interpolation
struct yac_interpolation * yac_interp_weights_get_interpolation(struct yac_interp_weights *weights, enum yac_interp_weights_reorder_type reorder, size_t collection_size, double frac_mask_fallback_value, double scaling_factor, double scaling_summand, char const *yaxt_exchanger_name, int is_source, int is_target)
Definition interp_weights.c:4835
yac_interp_weights_delete
void yac_interp_weights_delete(struct yac_interp_weights *weights)
Definition interp_weights.c:6414
yac_interp_weights_reorder_type
yac_interp_weights_reorder_type
Definition interp_weights.h:42
YAC_MAPPING_ON_TGT
@ YAC_MAPPING_ON_TGT
weights will be applied at target processes
Definition interp_weights.h:44
YAC_MAPPING_ON_SRC
@ YAC_MAPPING_ON_SRC
weights will be applied at source processes
Definition interp_weights.h:43
yac_interpolation_execute
void yac_interpolation_execute(struct yac_interpolation *interp, double ***src_fields, double **tgt_field)
Execute interpolation synchronously and write results to the target field.
Definition interpolation.c:251
yac_interpolation_delete
void yac_interpolation_delete(struct yac_interpolation *interp)
Free an interpolation object and release all resources.
Definition interpolation.c:371
YAC_FRAC_MASK_NO_VALUE
double const YAC_FRAC_MASK_NO_VALUE
Definition interpolation.c:25
yac_location
yac_location
Definition location.h:12
YAC_LOC_CORNER
@ YAC_LOC_CORNER
Definition location.h:15
YAC_LOC_CELL
@ YAC_LOC_CELL
Definition location.h:14
test_core.weights
weights
Definition test_core.py:63
test_core.src_field
src_field
Definition test_core.py:55
test_raw_exchange.i
int i
Definition test_raw_exchange.py:89
test_raw_exchange.tgt_field
tgt_field
Definition test_raw_exchange.py:64
grid_data
Definition test_dummy_coupling3_c.c:34
interp_method
Definition interp_method_internal.h:36
yac_basic_grid_data
Definition basic_grid_data.h:13
yac_basic_grid
Definition basic_grid.c:17
yac_dist_grid_pair
Definition dist_grid.c:119
yac_interp_field
Definition basic_grid.h:15
yac_interp_field::location
enum yac_location location
Definition basic_grid.h:16
yac_interp_grid
Definition interp_grid.c:22
yac_interp_grid::tgt_field
struct yac_interp_field tgt_field
Definition interp_grid.c:26
yac_interp_grid::num_src_fields
size_t num_src_fields
Definition interp_grid.c:27
yac_interp_grid::grid_pair
struct yac_dist_grid_pair * grid_pair
Definition interp_grid.c:25
yac_interp_grid::src_fields
struct yac_interp_field src_fields[]
Definition interp_grid.c:28
yac_interp_weights
Definition interp_weights.c:116
yac_interpolation
Definition interpolation.c:28
test_common.h
split_comm
static MPI_Comm split_comm
Definition test_interp_method_file_parallel.c:56
grid_names
static char const * grid_names[2]
Definition test_interp_method_ncc_parallel.c:29
num_reorder_types
size_t num_reorder_types
Definition test_interp_method_ncc_parallel.c:28
reorder_types
enum yac_interp_weights_reorder_type reorder_types[]
Definition test_interp_method_ncc_parallel.c:26
num_weight_types
size_t num_weight_types
Definition test_interp_method_spmap_parallel.c:32
weight_types
enum yac_interp_spmap_weight_type weight_types[]
Definition test_interp_method_spmap_parallel.c:28
tgt_field_data
static double tgt_field_data[MAX_COLLECTION_SIZE][NUM_TGT_POINTS]
Definition test_interp_operator_fixed.c:30
weight_type
weight_type
Definition test_interp_weights_parallel.c:30
coordinates_x
double coordinates_x[]
Definition test_interpolation_parallel3.c:104
num_cells
size_t num_cells[2]
Definition test_interpolation_parallel3.c:103
coordinates_y
double coordinates_y[]
Definition test_interpolation_parallel3.c:105
with_halo
int with_halo
Definition test_interpolation_parallel3.c:109
TEST_EXIT_CODE
#define TEST_EXIT_CODE
Definition tests.h:14
PUT_ERR
#define PUT_ERR(string)
Definition tests.h:10
grids
struct yac_basic_grid ** grids
Definition yac.c:152
yac_mpi.h
yac_int
YAC_INT yac_int
Definition yac_types.h:15