YAC 3.13.0
Yet Another Coupler
Loading...
Searching...
No Matches
toy_reg2d_atm.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#undef VERBOSE
6
7#define T106
8
9/* supported is
10 *
11 * T21, T31, T42, T63, T85, T106, T127, T159, T255
12 */
13
14#include <mpi.h>
15#include <stdlib.h>
16#include <stdio.h>
17#include <unistd.h>
18#include <math.h>
19#include "yac.h"
20#include "yac_core.h"
21#include "yac_utils.h"
22
23/* ------------------------------------------------- */
24
25/* For simplicity we define the same 4 fields that are in the
26 * coupling configuration */
27
28#include "toy_common.h"
29
30#define STR_USAGE "Usage: %s -c configFilename -x num_cells_x -y num_cells_y\n"
31#define YAC_ASSERT_ARGS(exp, msg) \
32 { \
33 if(!((exp))) { \
34 fprintf(stderr, "ERROR: %s\n" STR_USAGE, msg, argv[0]); \
35 exit(EXIT_FAILURE); \
36 } \
37 }
38
39static void parse_arguments(
40 int argc, char ** argv,
41 char const ** configFilename, size_t * num_cells_x, size_t * num_cells_y);
42
43int main (int argc, char *argv[]) {
44
45 // Initialisation of MPI
46
47 MPI_Init (0, NULL);
48
49 /* The initialisation phase includes the reading of the
50 * coupling configuration */
51#ifdef VERBOSE
52 printf (". main: calling yac_cinit\n");
53#endif
54
55 double tic, toc, time;
56 double time_min, time_max, time_ave;
57 double time_min_acc, time_max_acc, time_ave_acc;
58
59 MPI_Barrier(MPI_COMM_WORLD);
60
61 tic=MPI_Wtime();
62
63#if defined T21
64
65#define NUM_CELLS_X (64)
66#define NUM_CELLS_Y (32)
67
68#elif defined T31
69
70#define NUM_CELLS_X (96)
71#define NUM_CELLS_Y (48)
72
73#elif defined T42
74
75#define NUM_CELLS_X (128)
76#define NUM_CELLS_Y (64)
77
78#elif defined T63
79
80#define NUM_CELLS_X (192)
81#define NUM_CELLS_Y (96)
82
83#elif defined T85
84
85#define NUM_CELLS_X (256)
86#define NUM_CELLS_Y (128)
87
88#elif defined T106
89
90#define NUM_CELLS_X (320)
91#define NUM_CELLS_Y (160)
92
93#elif defined T127
94
95#define NUM_CELLS_X (384)
96#define NUM_CELLS_Y (192)
97
98#elif defined T159
99
100#define NUM_CELLS_X (480)
101#define NUM_CELLS_Y (240)
102
103#elif defined T255
104
105#define NUM_CELLS_X (768)
106#define NUM_CELLS_Y (384)
107
108#else
109
110#define NUM_CELLS_X (256)
111#define NUM_CELLS_Y (128)
112
113#endif
114
115 char const * configFilename = "toy_atm_ocn.yaml"; // default configuration file
116 size_t num_cells_x = NUM_CELLS_X; // default horizontal resolution
117 size_t num_cells_y = NUM_CELLS_Y; // default vertical resolution
118 parse_arguments(argc, argv, &configFilename, &num_cells_x, &num_cells_y);
119 yac_cinit ();
120 yac_cread_config_yaml(configFilename);
121
122 /* The usual component definition, here for two sequential components on the same process */
123
124#ifdef VERBOSE
125 printf (". main: calling yac_cdef_comp\n");
126#endif
127
128 int comp_id;
129 char * comp_name = "ATMOS";
130
131 yac_cdef_comp ( comp_name, &comp_id );
132#ifdef VERBOSE
133 printf ( ". main: defined %s with local comp ID %i \n", "toy-reg2d-atmosphere", comp_id );
134#endif
135
136 MPI_Comm local_comm;
137
138 yac_cget_comp_comm(comp_id, &local_comm);
139
140 int rank, size;
141
142 MPI_Comm_rank(local_comm,&rank);
143 MPI_Comm_size(local_comm,&size);
144
145
146 int point_id;
147
148 int field_ids[4];
149 int grid_id;
150
151 /* Grid definition for toy-reg2d-atmosphere */
152
153 int total_nbr_cells[2];
154 int num_procs[2];
155 total_nbr_cells[0] = (int)num_cells_x;
156 total_nbr_cells[1] = (int)num_cells_y;
157
158 yac_generate_reg2d_decomp(total_nbr_cells, size, num_procs);
159
160 int block_pos[2];
161 int block_size[2];
162 int nbr_cells[2];
163 int cyclic[2];
164 int nbr_points[2];
165 block_pos[0] = rank%num_procs[0];
166 block_pos[1] = rank/num_procs[0];
167 block_size[0] = (NUM_CELLS_X + num_procs[0] - 1)/num_procs[0];
168 block_size[1] = (NUM_CELLS_Y + num_procs[1] - 1)/num_procs[1];
169 nbr_cells[0] = MIN(block_size[0], NUM_CELLS_X - block_size[0] * block_pos[0]);
170 nbr_cells[1] = MIN(block_size[1], NUM_CELLS_Y - block_size[1] * block_pos[1]);
171 cyclic[0] = num_procs[0] == 1;
172 cyclic[1] = 0;
173
174 // halo
175 if (!cyclic[0])
176 nbr_cells[0] += 2;
177 if (num_procs[1] > 1) {
178
179 nbr_cells[1] += 1;
180
181 if ((block_pos[1] != 0) && (block_pos[1] != num_procs[1]-1))
182 nbr_cells[1] += 1;
183 }
184
185 if (cyclic[0])
186 nbr_points[0] = nbr_cells[0];
187 else
188 nbr_points[0] = nbr_cells[0]+1;
189
190 nbr_points[1] = nbr_cells[1] + 1;
191
192 // the grid
193
194 double * global_x_vertices = malloc((NUM_CELLS_X+3) * sizeof(*global_x_vertices));
195 double * global_y_vertices = malloc((NUM_CELLS_Y+1) * sizeof(*global_y_vertices));;
196 double * x_vertices;
197 double * y_vertices;
198 double * x_points;
199 double * y_points;
200
201 for (int i = 0; i < NUM_CELLS_X+3; ++i)
202 global_x_vertices[i] =
203 0.0 + (2.0 * M_PI / ((double)NUM_CELLS_X)) * (double)(i-1);
204 for (int i = 0; i < NUM_CELLS_Y; ++i)
205 global_y_vertices[i] =
206 -M_PI_2 + (M_PI / ((double)NUM_CELLS_Y)) * (double)i;
207 global_y_vertices[NUM_CELLS_Y] = M_PI_2;
208
209 if (cyclic[0])
210 x_vertices = x_points = global_x_vertices + 1;
211 else
212 x_vertices = x_points = global_x_vertices + block_size[0] * block_pos[0];
213
214 y_vertices = y_points = global_y_vertices + block_size[1] * block_pos[1];
215
216 if (block_pos[1] != 0) {
217 y_vertices -= 1;
218 y_points -= 1;
219 }
220
221 yac_cdef_grid_reg2d ( "grid1", nbr_points, cyclic, x_vertices, y_vertices, &grid_id);
222
223 int * global_global_cell_id = malloc(NUM_CELLS_Y * (NUM_CELLS_X + 2) * sizeof(*global_global_cell_id));
224 int * global_global_cell_id_rank = malloc(NUM_CELLS_Y * (NUM_CELLS_X + 2) * sizeof(*global_global_cell_id_rank));
225 int * global_global_corner_id = malloc((NUM_CELLS_Y + 1) * (NUM_CELLS_X + 3) * sizeof(*global_global_corner_id));
226
227 for (int j = 0, id = 0; j < NUM_CELLS_Y; ++j) {
228
229 for (int i = 1; i <= NUM_CELLS_X; ++i) {
230
231 global_global_cell_id[j * (NUM_CELLS_X + 2) + i] = id;
232 global_global_cell_id_rank[j * (NUM_CELLS_X + 2) + i] =
233 (id%NUM_CELLS_X) / block_size[0] + ((id/NUM_CELLS_X) / block_size[1]) * num_procs[0];
234 id++;
235 }
236 }
237 for (int i = 0; i < NUM_CELLS_Y; ++i) {
238 global_global_cell_id[i * (NUM_CELLS_X + 2) + 0] = global_global_cell_id[i * (NUM_CELLS_X + 2) + NUM_CELLS_X];
239 global_global_cell_id_rank[i * (NUM_CELLS_X + 2) + 0] = global_global_cell_id_rank[i * (NUM_CELLS_X + 2) + NUM_CELLS_X];
240 global_global_cell_id[i * (NUM_CELLS_X + 2) + NUM_CELLS_X+1] = global_global_cell_id[i * (NUM_CELLS_X + 2) + 1];
241 global_global_cell_id_rank[i * (NUM_CELLS_X + 2) + NUM_CELLS_X+1] = global_global_cell_id_rank[i * (NUM_CELLS_X + 2) + 1];
242 }
243
244 if (num_procs[0] == 1) {
245
246 for (int j = 0, id = 0; j < NUM_CELLS_Y + 1; ++j) {
247
248 for (int i = 1; i <= NUM_CELLS_X; ++i)
249 global_global_corner_id[j * (NUM_CELLS_X + 3) + i] = id++;
250 }
251
252 } else {
253
254 for (int j = 0, id = 0; j < NUM_CELLS_Y + 1; ++j) {
255
256 for (int i = 1; i <= NUM_CELLS_X + 1; ++i)
257 global_global_corner_id[j * (NUM_CELLS_X + 3) + i] = id++;
258 }
259
260 for (int i = 0; i < NUM_CELLS_Y + 1; ++i) {
261
262 global_global_corner_id[i * (NUM_CELLS_X + 3) + 0] = global_global_corner_id[i * (NUM_CELLS_X + 3) + NUM_CELLS_X + 1];
263 global_global_corner_id[i * (NUM_CELLS_X + 3) + NUM_CELLS_X + 2] = global_global_corner_id[i * (NUM_CELLS_X + 3) + 1];
264 }
265 }
266
267 for (int i = 0; i < (NUM_CELLS_Y)*(NUM_CELLS_X + 2); ++i)
268 if(global_global_cell_id_rank[i] == rank) global_global_cell_id_rank[i] = -1;
269
270 int * local_global_cell_id = malloc(nbr_cells[1] * nbr_cells[0] * sizeof(*local_global_cell_id));
271 int * local_global_cell_id_rank = malloc(nbr_cells[1] * nbr_cells[0] * sizeof(*local_global_cell_id_rank));
272 int * local_global_corner_id = malloc(nbr_points[1] * nbr_points[0] * sizeof(*local_global_corner_id));
273 int * local_global_corner_id_rank = malloc(nbr_points[1] * nbr_points[0] * sizeof(*local_global_corner_id_rank));
274
275 int offset_x, offset_y;
276
277 if (cyclic[0])
278 offset_x = 1;
279 else
280 offset_x = block_size[0] * block_pos[0];
281 offset_y = block_size[1] * block_pos[1] - (block_pos[1] != 0);
282
283 for (int j = 0; j < nbr_cells[1]; ++j) {
284 for (int i = 0; i < nbr_cells[0]; ++i) {
285 local_global_cell_id[j * nbr_cells[0] + i] = global_global_cell_id[(j+offset_y) * (NUM_CELLS_X + 2) + i+offset_x];
286 local_global_cell_id_rank[j * nbr_cells[0] + i] = global_global_cell_id_rank[(j+offset_y) * (NUM_CELLS_X + 2) +i+offset_x];
287 }
288 }
289
290 free(global_global_cell_id_rank);
291 free(global_global_cell_id);
292
293 for (int j = 0; j < nbr_points[1]; ++j)
294 for (int i = 0; i < nbr_points[0]; ++i)
295 local_global_corner_id[j * nbr_points[0] + i] = global_global_corner_id[(j+offset_y) * (NUM_CELLS_X + 3) +i+offset_x];
296
297 free(global_global_corner_id);
298
299 for (int j = 0; j < nbr_points[1]; ++j)
300 for (int i = 0; i < nbr_points[0]; ++i)
301 local_global_corner_id_rank[j * nbr_points[0] + i] = -1;
302
303 for (int i = 0; i < nbr_points[0]-2; ++i) {
304 local_global_corner_id_rank[0 * nbr_points[0] + i] = local_global_cell_id_rank[0 * nbr_cells[0] + i];
305 local_global_corner_id_rank[(nbr_points[1]-1) * nbr_points[0] + i] = local_global_cell_id_rank[(nbr_cells[1]-1) * nbr_cells[0] + i];
306 }
307 for (int j = 0; j < nbr_points[1]-2; ++j) {
308 local_global_corner_id_rank[j * nbr_points[0] + 0] = local_global_cell_id_rank[j * nbr_cells[0] + 0];
309 local_global_corner_id_rank[j * nbr_points[0] + nbr_points[0]-1] = local_global_cell_id_rank[j * nbr_cells[0] + nbr_cells[0]-1];
310 }
311 local_global_corner_id_rank[(nbr_points[1]-1)* nbr_points[0] + nbr_points[0]-1] = local_global_cell_id_rank[(nbr_cells[1]-1) * nbr_cells[0] + nbr_cells[0]-1];
312 local_global_corner_id_rank[(nbr_points[1]-2)* nbr_points[0] + 0] = local_global_cell_id_rank[(nbr_cells[1]-2) * nbr_cells[0] + 0];
313 local_global_corner_id_rank[0 * nbr_points[0] + nbr_points[0]-2] = local_global_cell_id_rank[0 * nbr_cells[0] + nbr_cells[0]-2];
314 local_global_corner_id_rank[(nbr_points[1]-2) * nbr_points[0] + nbr_points[0]-1] = local_global_cell_id_rank[(nbr_cells[1]-2) * nbr_cells[0] + nbr_cells[0]-1];
315 local_global_corner_id_rank[(nbr_points[1]-1) * nbr_points[0] + nbr_points[0]-2] = local_global_cell_id_rank[(nbr_cells[1]-1) * nbr_cells[0] + nbr_cells[0]-2];
316
317 int * cell_core_mask = malloc(nbr_cells[1] * nbr_cells[0] * sizeof(*local_global_corner_id_rank));
318 int * corner_core_mask = malloc(nbr_points[1] * nbr_points[0] * sizeof(*local_global_corner_id_rank));
319
320 for (int i = 0; i < nbr_cells[1]*nbr_cells[0]; ++i)
321 cell_core_mask[i] = local_global_cell_id_rank[i] == -1;
322 for (int i = 0; i < nbr_points[1]*nbr_points[0]; ++i)
323 corner_core_mask[i] = local_global_corner_id_rank[i] == -1;
324
325 yac_cset_global_index(local_global_cell_id, YAC_LOCATION_CELL, grid_id);
327 yac_cset_global_index(local_global_corner_id, YAC_LOCATION_CORNER, grid_id);
329
331 grid_id, nbr_cells, YAC_LOCATION_CELL, x_points, y_points, &point_id);
332
333 /* Field definition for toy-reg2d-atmosphere */
334
335 for ( int i = 0; i < num_fields; i++ )
337 fieldName[i], comp_id, &point_id, 1, 1, "2", YAC_TIME_UNIT_SECOND,
338 &field_ids[i]);
339
340 toc=MPI_Wtime();
341 time = toc-tic;
342
343 MPI_Reduce ( &time, &time_min, 1, MPI_DOUBLE, MPI_MIN, 0, local_comm);
344 MPI_Reduce ( &time, &time_max, 1, MPI_DOUBLE, MPI_MAX, 0, local_comm);
345 MPI_Reduce ( &time, &time_ave, 1, MPI_DOUBLE, MPI_SUM, 0, local_comm);
346
347 time_ave /= (double) size;
348
349 if ( rank == 0 )
350 printf ("toy-reg2d-atmosphere: Time for initialisation %f %f %f \n", time_min, time_max, time_ave );
351
352 /* Search. */
353
354#ifdef VERBOSE
355 printf (". main: calling yac_cenddef\n");
356#endif
357
358 MPI_Barrier(MPI_COMM_WORLD);
359
360 tic=MPI_Wtime();
361
362 yac_cenddef ( );
363
364 toc=MPI_Wtime();
365 time = toc-tic;
366
367 MPI_Reduce ( &time, &time_min, 1, MPI_DOUBLE, MPI_MIN, 0, local_comm);
368 MPI_Reduce ( &time, &time_max, 1, MPI_DOUBLE, MPI_MAX, 0, local_comm);
369 MPI_Reduce ( &time, &time_ave, 1, MPI_DOUBLE, MPI_SUM, 0, local_comm);
370
371 time_ave /= (double) size;
372
373 if ( rank == 0 )
374 printf ("toy-icon-atmosphere: Time for search %f %f %f \n", time_min, time_max, time_ave );
375
376 unsigned num_cells = nbr_cells[0] * nbr_cells[1];
377
378 double * cell_data_field = malloc(num_cells * sizeof(*cell_data_field));
379 double * cell_out_conserv_data = malloc(num_cells * sizeof(*cell_out_conserv_data));
380 double * cell_out_hcsbb_data = malloc(num_cells * sizeof(*cell_out_hcsbb_data));
381 double * cell_in_conserv_data = malloc(num_cells * sizeof(*cell_in_conserv_data));
382 double * cell_in_hcsbb_data = malloc(num_cells * sizeof(*cell_in_hcsbb_data));
383 double * cell_in_conserv_err_abs = malloc(num_cells * sizeof(*cell_in_conserv_err_abs));
384 double * cell_in_hcsbb_err_abs = malloc(num_cells * sizeof(*cell_in_hcsbb_err_abs));
385 double * cell_in_conserv_err_rel = malloc(num_cells * sizeof(*cell_in_conserv_err_rel));
386 double * cell_in_hcsbb_err_rel = malloc(num_cells * sizeof(*cell_in_hcsbb_err_rel));
387
388 int err;
389 int info;
390
391 unsigned * cell_corners = malloc(num_cells * 4 * sizeof(*cell_corners));
392 unsigned * num_points_per_cell = malloc(num_cells * sizeof(*num_points_per_cell));
393
394 for (unsigned i = 0; i < num_cells; ++i) {
395
396 {
397 unsigned x_index, y_index;
398
399 y_index = i / nbr_cells[0];
400 x_index = i - y_index * nbr_cells[0];
401
402 if (!cyclic[0]) {
403
404 cell_corners[i*4+0] = y_index * (nbr_cells[0] + 1) + x_index;
405 cell_corners[i*4+1] = y_index * (nbr_cells[0] + 1) + x_index + 1;
406 cell_corners[i*4+2] = (y_index + 1) * (nbr_cells[0] + 1) + x_index + 1;
407 cell_corners[i*4+3] = (y_index + 1) * (nbr_cells[0] + 1) + x_index;
408
409 } else {
410
411 cell_corners[i*4+0] = y_index * nbr_cells[0] + x_index;
412 if (x_index + 1 != (unsigned)(nbr_cells[0])) {
413 cell_corners[i*4+1] = y_index * nbr_cells[0] + x_index + 1;
414 cell_corners[i*4+2] = (y_index + 1) * nbr_cells[0] + x_index + 1;
415 } else {
416 cell_corners[i*4+1] = y_index * nbr_cells[0];
417 cell_corners[i*4+2] = (y_index + 1) * nbr_cells[0];
418 }
419 cell_corners[i*4+3] = (y_index + 1) * nbr_cells[0] + x_index;
420 }
421 }
422 num_points_per_cell[i] = 4;
423 }
424
425 for (unsigned i = 0; i < num_cells; ++i) {
426
427 double curr_x, curr_y;
428
429 curr_x = (x_vertices[cell_corners[i*4+0]%nbr_points[0]] +
430 x_vertices[cell_corners[i*4+1]%nbr_points[0]]) * 0.5;
431 curr_y = (y_vertices[cell_corners[i*4+1]/nbr_points[0]] +
432 y_vertices[cell_corners[i*4+2]/nbr_points[0]]) * 0.5;
433
434 cell_data_field[i] = yac_test_func(curr_x, curr_y);
435 }
436
437 for (unsigned i = 0; i < num_cells; ++i) {
438
439 cell_out_conserv_data[i] = cell_data_field[i];
440 cell_out_hcsbb_data[i] = cell_data_field[i];
441
442 cell_in_conserv_data[i] = -999;
443 cell_in_hcsbb_data[i] = -999;
444 }
445
446 double * point_data = malloc(nbr_points[0] * nbr_points[1] * 3 * sizeof(*point_data));
447 for (int i = 0; i < nbr_points[1]; ++i) {
448 for (int j = 0; j < nbr_points[0]; ++j) {
449
450 LLtoXYZ(x_vertices[j], y_vertices[i], &point_data[3*(i * nbr_points[0] + j)]);
451 }
452 }
453
454 time_min_acc = 0.0;
455 time_max_acc = 0.0;
456 time_ave_acc = 0.0;
457
458 for (int step = 0; step < num_steps; ++step) {
459
460 MPI_Barrier(MPI_COMM_WORLD);
461
462 tic=MPI_Wtime();
463
464 {
465 double *point_set_data[1];
466 double **collection_data[1] = {point_set_data};
467
468 point_set_data[0] = cell_out_conserv_data;
469 yac_cput(field_ids[0], 1, collection_data, &info, &err);
470 YAC_ASSERT(info, "check coupling period")
471
472 point_set_data[0] = cell_out_hcsbb_data;
473 yac_cput(field_ids[1], 1, collection_data, &info, &err);
474 YAC_ASSERT(info, "check coupling period")
475 }
476
477 {
478 double *collection_data[1];
479
480 collection_data[0] = cell_in_conserv_data;
481 yac_cget(field_ids[2], 1, collection_data, &info, &err);
482 YAC_ASSERT(info, "check coupling period")
483
484 collection_data[0] = cell_in_hcsbb_data;
485 yac_cget(field_ids[3], 1, collection_data, &info, &err);
486 YAC_ASSERT(info, "check coupling period")
487 }
488
489 toc=MPI_Wtime();
490 time = toc-tic;
491
492 MPI_Reduce ( &time, &time_min, 1, MPI_DOUBLE, MPI_MIN, 0, local_comm);
493 MPI_Reduce ( &time, &time_max, 1, MPI_DOUBLE, MPI_MAX, 0, local_comm);
494 MPI_Reduce ( &time, &time_ave, 1, MPI_DOUBLE, MPI_SUM, 0, local_comm);
495
496 time_ave /= (double) size;
497
498 if ( rank == 0 ) {
499 time_min_acc += time_min;
500 time_max_acc += time_max;
501 time_ave_acc += time_ave;
502 }
503
504 //----------------------------------------------------------
505 // compute error
506 //----------------------------------------------------------
507
508 for ( unsigned i = 0; i < num_cells; ++i ) {
509
510 cell_in_conserv_err_abs[i] = fabs(cell_in_conserv_data[i] - cell_data_field[i]);
511 cell_in_hcsbb_err_abs[i] = fabs(cell_in_hcsbb_data[i] - cell_data_field[i]);
512 cell_in_conserv_err_rel[i] = fabs(cell_in_conserv_err_abs[i] / cell_data_field[i]);
513 cell_in_hcsbb_err_rel[i] = fabs(cell_in_hcsbb_err_abs[i] / cell_data_field[i]);
514 }
515
516#if 0
517
518 double f1err_1 = 0.0;
519 double f1err_100 = 0.0;
520 int if1err_1 = 0;
521 int if1err_100 = 0;
522
523 double f2err_1 = 0.0;
524 double f2err_100 = 0.0;
525 int if2err_1 = 0;
526 int if2err_100 = 0;
527
528 for ( int i = 0; i < num_cells; ++i ) {
529
530 if ( cell_in_conserv_err_rel[i] < 100.0 ) {
531 f1err_1 += cell_in_conserv_err_rel[i];
532 if1err_1++;
533 } else {
534 f1err_100 += cell_in_conserv_err_rel[i];
535 if1err_100++;
536 }
537
538 if ( cell_in_hcsbb_err_rel[i] < 100.0 ) {
539 f2err_1 += cell_in_hcsbb_err_rel[i];
540 if2err_1++;
541 } else {
542 f2err_100 += cell_in_hcsbb_err_rel[i];
543 if2err_100++;
544 }
545 }
546
547 if (if1err_1 > 0.0) f1err_1 = f1err_1 / (double)if1err_1;
548 if (if1err_1 > 0.0) f2err_1 = f2err_1 / (double)if2err_1;
549
550 if (if1err_100 > 0.0) f1err_100 = f1err_100 / (double)if1err_100;
551 if (if2err_100 > 0.0) f2err_100 = f2err_100 / (double)if2err_100;
552
553 printf ("avg. rel err. for 1st field %f\n", f1err_1);
554 printf ("avg. rel err. for 2nd field %f\n", f2err_1);
555
556 printf ("extreme avg. rel err. for 1st field %f\n", f1err_100);
557 printf ("extreme avg. rel err. for 2nd field %f\n", f2err_100);
558#endif
559
560 //----------------------------------------------------------
561 // write field to vtk output file
562 //----------------------------------------------------------
563
564 char vtk_filename[36];
565
566 sprintf(vtk_filename, "toy-reg2d-atmosphere_out_%d_%d.vtk", rank, step);
567
568 YAC_VTK_FILE *vtk_file =
569 yac_vtk_open(vtk_filename, "toy-reg2d-atmosphere_out");
570
572 vtk_file, point_data, nbr_points[0]*nbr_points[1]);
573
575 vtk_file, cell_corners, num_points_per_cell, num_cells);
576
578 vtk_file, corner_core_mask, nbr_points[0]*nbr_points[1],
579 "corner_core_mask");
581 vtk_file, local_global_corner_id, nbr_points[0]*nbr_points[1],
582 "global_corner_id");
584 vtk_file, cell_core_mask, num_cells, "cell_core_mask");
586 vtk_file, local_global_cell_id, num_cells, "global_cell_id");
587
589 vtk_file, cell_in_conserv_data, num_cells, "cell_in_conserv");
591 vtk_file, cell_in_hcsbb_data, num_cells, "cell_in_hcsbb");
593 vtk_file, cell_out_conserv_data, num_cells, "cell_out_conserv");
595 vtk_file, cell_out_hcsbb_data, num_cells, "cell_out_hcsbb");
597 vtk_file, cell_in_conserv_err_abs, num_cells, "cell_in_conserv_err_abs");
599 vtk_file, cell_in_hcsbb_err_abs, num_cells, "cell_in_hcsbb_err_abs");
601 vtk_file, cell_in_conserv_err_rel, num_cells, "cell_in_conserv_err_rel");
603 vtk_file, cell_in_hcsbb_err_rel, num_cells, "cell_in_hcsbb_err_rel");
604
605 yac_vtk_close(vtk_file);
606
607 for ( unsigned i = 0; i < num_cells; ++i ) {
608 cell_out_conserv_data[i] = cell_in_conserv_data[i];
609 cell_out_hcsbb_data[i] = cell_in_hcsbb_data[i];
610 cell_in_conserv_data[i] = -10;
611 cell_in_hcsbb_data[i] = -10;
612 }
613 }
614
615 if ( rank == 0 ) {
616 time_min_acc /= (double) num_steps;
617 time_max_acc /= (double) num_steps;
618 time_ave_acc /= (double) num_steps;
619 printf ("toy-icon-atmosphere: Time for ping-pong %f %f %f \n", time_min_acc, time_max_acc, time_ave_acc );
620 }
621
623
624 MPI_Finalize();
625
626 free(cell_in_hcsbb_err_rel);
627 free(cell_in_conserv_err_rel);
628 free(cell_in_hcsbb_err_abs);
629 free(cell_in_conserv_err_abs);
630 free(cell_in_conserv_data);
631 free(cell_in_hcsbb_data);
632 free(cell_out_conserv_data);
633 free(cell_out_hcsbb_data);
634 free(cell_data_field);
635
636 free(num_points_per_cell);
637 free(cell_corners);
638 free(point_data);
639 free(corner_core_mask);
640 free(local_global_corner_id_rank);
641 free(local_global_corner_id);
642 free(cell_core_mask);
643 free(local_global_cell_id_rank);
644 free(local_global_cell_id);
645 free(global_x_vertices);
646 free(global_y_vertices);
647
648 return EXIT_SUCCESS;
649}
650
651static void parse_arguments(
652 int argc, char ** argv,
653 char const ** configFilename, size_t * num_cells_x, size_t * num_cells_y) {
654
655 int opt;
656 while ((opt = getopt(argc, argv, "c:x:y:")) != -1) {
657 YAC_ASSERT((opt == 'c') || (opt == 'x') || (opt == 'y'), "invalid command argument")
658 switch (opt) {
659 default:
660 case 'c':
661 *configFilename = optarg;
662 break;
663 case 'x':
664 *num_cells_x = atoi(optarg);
665 YAC_ASSERT_ARGS(*num_cells_x > 0, "invalid horizontal resolution");
666 break;
667 case 'y':
668 *num_cells_y = atoi(optarg);
669 YAC_ASSERT_ARGS(*num_cells_y > 0, "invalid vertical resolution");
670 break;
671 }
672 }
673}
unsigned num_fields
#define YAC_ASSERT(exp, msg)
void yac_generate_reg2d_decomp(int num_points[2], int total_num_procs, int *num_procs)
const char * fieldName[]
int point_id
double yac_test_func(double lon, double lat)
size_t num_cells[2]
unsigned cyclic[2]
#define MIN(a, b)
Definition toy_common.h:29
const int num_steps
Definition toy_common.h:10
int info
int * cell_core_mask
int grid_id
int comp_id
#define YAC_ASSERT_ARGS(exp, msg)
#define NUM_CELLS_X
#define NUM_CELLS_Y
static void parse_arguments(int argc, char **argv, char const **configFilename, size_t *num_cells_x, size_t *num_cells_y)
static void LLtoXYZ(double lon, double lat, double p_out[])
Definition toy_scrip.c:587
void yac_vtk_write_cell_scalars_double(YAC_VTK_FILE *vtk_file, double *scalars, unsigned num_cells, char *name)
Definition vtk_output.c:264
void yac_vtk_write_cell_scalars_int(YAC_VTK_FILE *vtk_file, int *scalars, unsigned num_cells, char *name)
Definition vtk_output.c:250
void yac_vtk_write_point_data(YAC_VTK_FILE *vtk_file, double *point_data, unsigned num_points)
Definition vtk_output.c:69
void yac_vtk_write_point_scalars_int(YAC_VTK_FILE *vtk_file, int *scalars, unsigned num_points, char *name)
Definition vtk_output.c:278
void yac_vtk_close(YAC_VTK_FILE *vtk_file)
Definition vtk_output.c:403
void yac_vtk_write_cell_data(YAC_VTK_FILE *vtk_file, unsigned *cell_corners, unsigned *num_points_per_cell, unsigned num_cells)
Definition vtk_output.c:88
YAC_VTK_FILE * yac_vtk_open(const char *filename, const char *title)
Definition vtk_output.c:45
void yac_cenddef(void)
Definition yac.c:4400
void yac_cset_global_index(int const *global_index, int location, int grid_id)
Definition yac.c:5049
int const YAC_LOCATION_CELL
Definition yac.c:34
void yac_cinit(void)
Definition yac.c:510
void yac_cfinalize()
Finalises YAC.
Definition yac.c:740
int const YAC_LOCATION_CORNER
Definition yac.c:35
void yac_cget_comp_comm(int comp_id, MPI_Comm *comp_comm)
Definition yac.c:856
void yac_cget(int const field_id, int collection_size, double **recv_field, int *info, int *ierr)
Definition yac.c:2759
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:4831
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:1103
int const YAC_TIME_UNIT_SECOND
Definition yac.c:59
void yac_cput(int const field_id, int const collection_size, double ***const send_field, int *info, int *ierr)
Definition yac.c:3721
void yac_cread_config_yaml(const char *yaml_filename)
Definition yac.c:595
void yac_cset_core_mask(int const *is_core, int location, int grid_id)
Definition yac.c:5065
void yac_cdef_comp(char const *comp_name, int *comp_id)
Definition yac.c:1013
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:1396