YAC 3.21.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_dummy_coupling8_c.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 <mpi.h>
6
7#include <stdlib.h>
8#include <stdio.h>
9#include <string.h>
10#include <math.h>
11#include <float.h>
12#include "tests.h"
13#include "test_common.h"
14#include "yac.h"
15
21#define YAC_RAD (0.01745329251994329576923690768489) // M_PI / 180
22
23enum {
29 TGT_DT = 2,
30 SRC_DT = 1,
31};
32
33#define FRAC_MASK_VALUE (1337.0)
34
35static void utest_aggregation();
36static void utest_changing_frac_mask();
37static void utest_self_coupling();
38
39int main (void) {
40
41
42 utest_aggregation();
43 utest_changing_frac_mask();
44 utest_self_coupling();
45
47
48 return TEST_EXIT_CODE;
49}
50
51static void utest_init_ref_recv_field(
52 double ref_recv_field[][COLLECTION_SIZE][NUM_POINTS]) {
53
54 for (int i = 0; i < 3; ++i)
55 for (int j = 0; j < COLLECTION_SIZE; ++j)
56 for (int k = 0; k < NUM_POINTS; ++k)
57 ref_recv_field[i][j][k] = 0.0;
58 for (int j = 0; j < COLLECTION_SIZE; ++j) {
59 for (int k = 0; k < NUM_POINTS; ++k) {
60 ref_recv_field[3][j][k] = DBL_MAX;
61 ref_recv_field[4][j][k] = -DBL_MAX;
62 }
63 }
64}
65
66static void utest_init_temp_frac_mask(
67 double temp_frac_mask[][COLLECTION_SIZE][NUM_POINTSETS][NUM_POINTS]) {
68
69 for (int i = 0; i < NUM_FIELDS; ++i)
70 for (int j = 0; j < COLLECTION_SIZE; ++j)
71 for (int k = 0; k < NUM_POINTS; ++k)
72 temp_frac_mask[i][j][0][k] = 0.0;
73}
74
75static void utest_aggregation() {
76
77 yac_cinit();
79 yac_cdef_datetime("1850-01-01T00:00:00", "1850-01-03T00:00:00");
80
81 int size, rank;
82 MPI_Comm_rank ( MPI_COMM_WORLD, &rank );
83 MPI_Comm_size ( MPI_COMM_WORLD, &size );
84
85 if (size != 2) {
86 fputs("wrong number of processes (has to be 2)\n", stderr);
87 exit(EXIT_FAILURE);
88 }
89
90 int is_target = rank == 1;
91
92 // define local component
93 int comp_id;
94 yac_cdef_comp((is_target)?"target_comp":"source_comp", &comp_id);
95
96 // define grid (both components use an identical grid)
97 int grid_id;
99 (is_target)?"target_grid":"source_grid", (int[2]){3,3}, (int[2]){0,0},
100 (double[]){0*YAC_RAD,1*YAC_RAD,2*YAC_RAD},
101 (double[]){0*YAC_RAD,1*YAC_RAD,2*YAC_RAD}, &grid_id);
102
103 // define points at the vertices of the grid
104 int point_id;
106 grid_id, (int[2]){3,3}, YAC_LOCATION_CORNER,
107 (double[]){0*YAC_RAD,1*YAC_RAD,2*YAC_RAD},
108 (double[]){0*YAC_RAD,1*YAC_RAD,2*YAC_RAD}, &point_id);
109
110 // define fields
111 int field_ids[NUM_FIELDS];
112 const char * fieldName[NUM_FIELDS] =
113 {"time_op_accu_field",
114 "time_op_avg_field",
115 "time_op_none_field",
116 "time_op_min_field",
117 "time_op_max_field"};
118 for (int field_idx = 0; field_idx < NUM_FIELDS; ++field_idx) {
120 fieldName[field_idx], comp_id, &point_id, NUM_POINTSETS,
121 COLLECTION_SIZE, (is_target)?"2":"1", YAC_TIME_UNIT_SECOND,
122 &field_ids[field_idx]);
124 (is_target)?"target_comp":"source_comp",
125 (is_target)?"target_grid":"source_grid",
126 fieldName[field_idx], FRAC_MASK_VALUE);
127 }
128
129 // define interpolation stacks
130 int interp_stack_nnn;
131 yac_cget_interp_stack_config(&interp_stack_nnn);
133 interp_stack_nnn, YAC_NNN_AVG, 1, 0.0, 0.0);
134
135 // define couplings
136 int reduction_type[NUM_FIELDS] =
142 for (int field_idx = 0; field_idx < NUM_FIELDS; ++field_idx)
144 "source_comp", "source_grid", fieldName[field_idx],
145 "target_comp", "target_grid", fieldName[field_idx],
146 "4", YAC_TIME_UNIT_SECOND, reduction_type[field_idx],
147 interp_stack_nnn, 0, 0);
148 yac_cfree_interp_stack_config(interp_stack_nnn);
149
150 yac_cenddef ( );
151
152 double send_field[COLLECTION_SIZE][NUM_POINTSETS][NUM_POINTS] =
153 {{{ 1, 2, 3, 4, 5, 6, 7, 8, 9}},
154 {{10,11,12,13,14,15,16,17,18}},
155 {{19,20,21,22,23,24,25,26,27}}};
157 {{{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}},
158 {{1.0, 1.0, 0.7, 0.7, 0.5, 0.3, 0.3, 0.0, 0.0}},
159 {{1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0}}},
160 {{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0}},
161 {{0.7, 0.7, 0.7, 0.7, 0.5, 0.3, 0.3, 0.3, 0.3}},
162 {{1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0}}},
163 {{{0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0}},
164 {{0.3, 0.3, 0.3, 0.3, 0.5, 0.7, 0.7, 0.7, 0.7}},
165 {{1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}}},
166 {{{0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}},
167 {{0.0, 0.0, 0.3, 0.3, 0.5, 0.7, 0.7, 1.0, 1.0}},
168 {{1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}}}};
169 double temp_frac_mask[NUM_FIELDS][COLLECTION_SIZE][NUM_POINTSETS][NUM_POINTS];
170 double ref_recv_field[NUM_FIELDS][COLLECTION_SIZE][NUM_POINTS];
171 utest_init_ref_recv_field(ref_recv_field);
172 utest_init_temp_frac_mask(temp_frac_mask);
173
174 // do time steps
175 for (int t = 0; t < 8 * COUPLING_DT; ++t) {
176
177 if (!is_target) {
178
179 for (int field_idx = 0; field_idx < NUM_FIELDS; ++field_idx) {
180
181 int info, ierror;
183 field_ids[field_idx], COLLECTION_SIZE, &send_field[0][0][0],
184 &frac_mask[t%COUPLING_DT][0][0][0], &info, &ierror);
185
186 int ref_send_info =
187 (t%COUPLING_DT)?
188 ((field_idx == 2)?YAC_ACTION_NONE:YAC_ACTION_REDUCTION):
190 if (info != ref_send_info) PUT_ERR("error in yac_cput_: wrong info");
191 if (ierror) PUT_ERR("error in yac_cput_: wrong ierror");
192 }
193 }
194
195 // the first timestep is coupled directly
196 double scale = ((t == 0)?1.0:0.25);
197 for (int i = 0; i < COLLECTION_SIZE; ++i) {
198 for (int j = 0; j < NUM_POINTS; ++j) {
199
200 double frac_mask_value = frac_mask[t%COUPLING_DT][i][0][j];
201 double frac_send_field_value =
202 frac_mask_value * send_field[i][0][j];
203
204 if (frac_mask_value != 0.0) {
205
206 // update ref_recv_field (accu)
207 ref_recv_field[0][i][j] += frac_send_field_value;
208 temp_frac_mask[0][i][0][j] += frac_mask_value * scale;
209
210 // update ref_recv_field (avg)
211 ref_recv_field[1][i][j] += frac_send_field_value * scale;
212 temp_frac_mask[1][i][0][j] += frac_mask_value * scale;
213
214 // update ref_recv_field (minimum)
215 if (ref_recv_field[3][i][j] > frac_send_field_value) {
216 ref_recv_field[3][i][j] = frac_send_field_value;
217 temp_frac_mask[3][i][0][j] = frac_mask_value;
218 }
219
220 // update ref_recv_field (maximum)
221 if (ref_recv_field[4][i][j] < frac_send_field_value) {
222 ref_recv_field[4][i][j] = frac_send_field_value;
223 temp_frac_mask[4][i][0][j] = frac_mask_value;
224 }
225
226 }
227
228 // update ref_recv_field (none)
229 ref_recv_field[2][i][j] = frac_send_field_value;
230 temp_frac_mask[2][i][0][j] = frac_mask_value;
231 }
232 }
233
234 if (is_target) {
235
236 // target calls get every second timestep
237 if ((t % TGT_DT) == 0) {
238 for (int field_idx = 0; field_idx < NUM_FIELDS; ++field_idx) {
239
240 // initialise recv_field
241 double recv_field[COLLECTION_SIZE][NUM_POINTS];
242 for (int j = 0; j < COLLECTION_SIZE; ++j)
243 for (int k = 0; k < NUM_POINTS; ++k)
244 recv_field[j][k] = -1;
245
246 int info, ierror;
247 yac_cget_(
248 field_ids[field_idx], COLLECTION_SIZE, &recv_field[0][0],
249 &info, &ierror);
250
251 int ref_recv_info =
253 if (info != ref_recv_info) PUT_ERR("error in yac_cget_: wrong info");
254 if (ierror) PUT_ERR("error in yac_cget_: wrong ierror");
255
256 if (info == YAC_ACTION_COUPLING) {
257 for (int j = 0; j < COLLECTION_SIZE; ++j) {
258 for (int k = 0; k < NUM_POINTS; ++k) {
259 if (temp_frac_mask[field_idx][j][0][k] != 0.0) {
260 if (fabs(recv_field[j][k] -
261 (ref_recv_field[field_idx][j][k] /
262 temp_frac_mask[field_idx][j][0][k])) > 1e-6)
263 PUT_ERR("error in yac_cget_: wrong recv_field (unmasked)");
264 } else {
265 if (recv_field[j][k] != FRAC_MASK_VALUE)
266 PUT_ERR("error in yac_cget_: wrong recv_field (masked)");
267 }
268 }
269 }
270 } else {
271 for (int j = 0; j < COLLECTION_SIZE; ++j)
272 for (int k = 0; k < NUM_POINTS; ++k)
273 if (recv_field[j][k] != -1)
274 PUT_ERR("error in yac_cget_: wrong recv_field");
275 }
276 }
277 }
278 }
279
280 // update send_field
281 for (int i = 0; i < COLLECTION_SIZE; ++i)
282 for (int j = 0; j < NUM_POINTS; ++j)
283 send_field[i][0][j] += i - 1;
284
285 // clean ref_recv_field at every coupling timestep
286 if ((t % COUPLING_DT) == 0) {
287 utest_init_ref_recv_field(ref_recv_field);
288 utest_init_temp_frac_mask(temp_frac_mask);
289 }
290 }
291
292 yac_ccleanup();
293}
294
295static void utest_changing_frac_mask() {
296
297 yac_cinit();
299 yac_cdef_datetime("1850-01-01T00:00:00", "1850-01-03T00:00:00");
300
301 int size, rank;
302 MPI_Comm_rank ( MPI_COMM_WORLD, &rank );
303 MPI_Comm_size ( MPI_COMM_WORLD, &size );
304
305 if (size != 2) {
306 fputs("wrong number of processes (has to be 2)\n", stderr);
307 exit(EXIT_FAILURE);
308 }
309
310 int is_target = rank == 1;
311
312 // define local component
313 int comp_id;
314 yac_cdef_comp((is_target)?"target_comp":"source_comp", &comp_id);
315
316 // define grid (source has a 2x2 grid, while the target has a 1x1 grid)
317 int grid_id;
318 if (is_target)
320 "target_grid_2", (int[2]){2,2}, (int[2]){0,0},
321 (double[]){0*YAC_RAD,2*YAC_RAD},
322 (double[]){0*YAC_RAD,2*YAC_RAD}, &grid_id);
323 else
325 "source_grid_2", (int[2]){3,3}, (int[2]){0,0},
326 (double[]){0*YAC_RAD,1*YAC_RAD,2*YAC_RAD},
327 (double[]){0*YAC_RAD,1*YAC_RAD,2*YAC_RAD}, &grid_id);
328
329 // define points at the vertices of the grid
330 int point_id;
331 if (is_target)
333 grid_id, (int[2]){1,1}, YAC_LOCATION_CELL,
334 (double[]){1*YAC_RAD}, (double[]){1*YAC_RAD}, &point_id);
335 else
337 grid_id, (int[2]){2,2}, YAC_LOCATION_CELL,
338 (double[]){0.5*YAC_RAD,1.5*YAC_RAD},
339 (double[]){0.5*YAC_RAD,1.5*YAC_RAD}, &point_id);
340
341 // define field
342 int field_id;
344 "field", comp_id, &point_id, 1, 1, "1", YAC_TIME_UNIT_SECOND, &field_id);
346 (is_target)?"target_comp":"source_comp",
347 (is_target)?"target_grid_2":"source_grid_2",
348 "field", FRAC_MASK_VALUE);
349
350 // define interpolation stacks
351 int interp_stack_nnn;
352 yac_cget_interp_stack_config(&interp_stack_nnn);
354 interp_stack_nnn, YAC_NNN_AVG, 4, 0.0, 0.0);
355
356 // define couplings
358 "source_comp", "source_grid_2", "field",
359 "target_comp", "target_grid_2", "field",
361 interp_stack_nnn, 0, 0);
362 yac_cfree_interp_stack_config(interp_stack_nnn);
363
364 yac_cenddef();
365
366 enum {NUM_TIMESTEP = 5, SRC_FIELD_SIZE = 4};
367
368 double send_field[NUM_TIMESTEP][SRC_FIELD_SIZE] =
369 {{ 1, 2, 3, 4}, {10,11,12,13}, {20,21,22,23},
370 {30,31,32,33}, {40,41,42,43}};
371 double frac_mask[NUM_TIMESTEP][SRC_FIELD_SIZE] =
372 {{1.0,1.0,0.0,0.0}, {1.0,0.0,1.0,0.0}, {1.0,1.0,0.0,0.0},
373 {1.0,1.0,1.0,1.0}, {0.0,0.0,0.0,0.0}};
374
375 for (int t = 0; t < NUM_TIMESTEP; ++t) {
376
377 if (!is_target) {
378
379 int info, ierror;
381 field_id, 1, &send_field[t][0], &frac_mask[t][0], &info, &ierror);
382 if (info != YAC_ACTION_COUPLING) PUT_ERR("error in yac_cput_: wrong info");
383 if (ierror) PUT_ERR("error in yac_cput_: wrong ierror");
384 }
385
386 if (is_target) {
387
388 double const weight = 0.25;
389 double ref_recv_field = 0.0;
390 double frac_weight_sum = 0.0;
391 for (int i = 0; i < SRC_FIELD_SIZE; ++i) {
392 ref_recv_field += send_field[t][i] * frac_mask[t][i] * weight;
393 frac_weight_sum += weight * frac_mask[t][i];
394 }
395 ref_recv_field =
396 (frac_weight_sum < 1.e-9)?
397 FRAC_MASK_VALUE:(ref_recv_field/frac_weight_sum);
398
399 double recv_field;
400 int info, ierror;
401 yac_cget_(field_id, 1, &recv_field, &info, &ierror);
402
404 PUT_ERR("error in yac_cget_: wrong info");
405 if (ierror) PUT_ERR("error in yac_cget_: wrong ierror");
406 if (fabs(recv_field - ref_recv_field) > 1.e-9)
407 PUT_ERR("error in yac_cget_: wrong recv_field");
408 }
409 }
410
411 yac_ccleanup();
412}
413
414static void utest_self_coupling() {
415
416 yac_cinit();
418 yac_cdef_datetime("1850-01-01T00:00:00", "1850-01-03T00:00:00");
419
420 int comm_rank, comm_size;
421 MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
422 MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
423
424 if (comm_size != 2) {
425 fputs("wrong number of processes (has to be 2)\n", stderr);
426 exit(EXIT_FAILURE);
427 }
428
429 int is_comp_b = (comm_rank == 1);
430
431 int comp_id;
432 yac_cdef_comp(is_comp_b ? "comp_b" : "comp_a", &comp_id);
433
434 enum { NX = 3, NY = 3, NUM_CELLS = NX * NY };
435 double x_verts[NX + 1], y_verts[NY + 1];
436 double x_cells[NX], y_cells[NY];
437 for (int i = 0; i <= NX; ++i) x_verts[i] = (double)i * YAC_RAD;
438 for (int j = 0; j <= NY; ++j) y_verts[j] = (double)j * YAC_RAD;
439 for (int i = 0; i < NX; ++i) x_cells[i] = (0.5 + i) * YAC_RAD;
440 for (int j = 0; j < NY; ++j) y_cells[j] = (0.5 + j) * YAC_RAD;
441
442 int grid_id, point_id;
444 is_comp_b ? "grid_b" : "grid_a",
445 (int[2]){NX + 1, NY + 1}, (int[2]){0, 0},
446 x_verts, y_verts, &grid_id);
448 grid_id, (int[2]){NX, NY}, YAC_LOCATION_CELL,
449 x_cells, y_cells, &point_id);
450
451 int src_field_id = -1, tgt_field_id = -1;
452 if (!is_comp_b) {
454 "src_field", comp_id, &point_id, 1, 1, "1",
455 YAC_TIME_UNIT_SECOND, &src_field_id);
457 "tgt_field", comp_id, &point_id, 1, 1, "1",
458 YAC_TIME_UNIT_SECOND, &tgt_field_id);
459 }
460
461 int interp_stack;
462 yac_cget_interp_stack_config(&interp_stack);
463 yac_cadd_interp_stack_config_nnn(interp_stack, YAC_NNN_AVG, 1, 0.0, 0.0);
464
465 // Regression: on rank 1, grid_a is absent. The bug produced two distinct
466 // yac_basic_grid_empty_new() pointers for the same name, crashing in
467 // yac_dist_grid_pair_new() with "both grids use that same name but point to
468 // different basic grids".
470 "comp_a", "grid_a", "src_field",
471 "comp_a", "grid_a", "tgt_field",
473 interp_stack, 0, 0);
474 yac_cfree_interp_stack_config(interp_stack);
475
476 yac_cenddef();
477
478 if (!is_comp_b) {
479 double src_data[NUM_CELLS], tgt_data[NUM_CELLS];
480 for (int i = 0; i < NUM_CELLS; ++i) {
481 src_data[i] = (double)i;
482 tgt_data[i] = -1.0;
483 }
484 int send_info, recv_info, ierr;
486 src_field_id, tgt_field_id, 1,
487 src_data, tgt_data, &send_info, &recv_info, &ierr);
488 if (send_info != YAC_ACTION_COUPLING)
489 PUT_ERR("wrong send_info in self-coupling exchange");
490 if (recv_info != YAC_ACTION_COUPLING)
491 PUT_ERR("wrong recv_info in self-coupling exchange");
492 if (ierr)
493 PUT_ERR("error in self-coupling exchange");
494 for (int i = 0; i < NUM_CELLS; ++i)
495 if (tgt_data[i] != src_data[i])
496 PUT_ERR("wrong value in self-coupled tgt_field");
497 }
498
499 yac_ccleanup();
500}
@ COLLECTION_SIZE
#define YAC_RAD
#define FRAC_MASK_VALUE
int point_id
#define TEST_EXIT_CODE
Definition tests.h:15
#define PUT_ERR(string)
Definition tests.h:10
const char * fieldName[]
Definition toy_common.h:5
int info
int grid_id
int ierror
int * field_id
int comp_id
int const YAC_REDUCTION_TIME_AVERAGE
Definition yac.c:59
void yac_cenddef(void)
Definition yac.c:5006
void yac_cenable_field_frac_mask(const char *comp_name, const char *grid_name, const char *field_name, double frac_mask_fallback_value)
Definition yac.c:1969
void yac_cget_(int const field_id, int const collection_size, double *recv_field, int *info, int *ierr)
Definition yac.c:3280
void yac_cdef_datetime(const char *start_datetime, const char *end_datetime)
Definition of job start and end datetime for the default YAC instance.
Definition yac.c:1028
int const YAC_LOCATION_CELL
Definition yac.c:38
int const YAC_REDUCTION_TIME_MINIMUM
Definition yac.c:60
void yac_cput_frac_(int const field_id, int const collection_size, double *send_field, double *send_frac_mask, int *info, int *ierr)
Definition yac.c:3621
void yac_ccleanup()
Clean-up default YAC instance (see Restarting YAC)
Definition yac.c:977
int const YAC_REDUCTION_TIME_MAXIMUM
Definition yac.c:61
void yac_cinit(void)
Definition yac.c:737
void yac_cfinalize()
Finalises YAC.
Definition yac.c:1007
int const YAC_LOCATION_CORNER
Definition yac.c:39
void yac_cexchange_(int const send_field_id, int const recv_field_id, int const collection_size, double *send_field, double *recv_field, int *send_info, int *recv_info, int *ierr)
Definition yac.c:4284
int const YAC_REDUCTION_TIME_ACCUMULATE
Definition yac.c:58
void yac_cdef_grid_reg2d(const char *grid_name, int nbr_vertices[2], int cyclic[2], double *x_vertices, double *y_vertices, int *grid_id)
Definition yac.c:5511
void yac_cdef_points_reg2d(int const grid_id, int const *nbr_points, int const located, double const *x_points, double const *y_points, int *point_id)
Definition yac.c:1540
int const YAC_TIME_UNIT_SECOND
Definition yac.c:64
int const YAC_ACTION_NONE
no data exchanges
Definition yac.c:47
void yac_cdef_calendar(int calendar)
Definition yac.c:1036
int const YAC_ACTION_REDUCTION
data reduction, but data exchange
Definition yac.c:48
int const YAC_PROLEPTIC_GREGORIAN
Definition yac.c:73
int const YAC_ACTION_COUPLING
data exchange
Definition yac.c:49
void yac_cfree_interp_stack_config(int interp_stack_config_id)
Definition yac.c:5834
int const YAC_NNN_AVG
Definition yac.c:84
void yac_cadd_interp_stack_config_nnn(int interp_stack_config_id, int type, size_t n, double max_search_distance, double scale)
Definition yac.c:5890
void yac_cget_interp_stack_config(int *interp_stack_config_id)
Definition yac.c:5818
void yac_cdef_comp(char const *comp_name, int *comp_id)
Definition yac.c:1345
void yac_cdef_field(char const *name, int const comp_id, int const *point_ids, int const num_pointsets, int collection_size, const char *timestep, int time_unit, int *field_id)
Definition yac.c:1929
void yac_cdef_couple(char const *src_comp_name, char const *src_grid_name, char const *src_field_name, char const *tgt_comp_name, char const *tgt_grid_name, char const *tgt_field_name, char const *coupling_timestep, int time_unit, int time_reduction, int interp_stack_config_id, int src_lag, int tgt_lag)
Definition yac.c:2475
int const YAC_REDUCTION_TIME_NONE
Definition yac.c:57