30#define STR_USAGE "Usage: %s -c configFilename -x num_cells_x -y num_cells_y\n"
31#define YAC_ASSERT_ARGS(exp, msg) \
34 fprintf(stderr, "ERROR: %s\n" STR_USAGE, msg, argv[0]); \
40 int argc,
char ** argv,
41 char const ** configFilename,
size_t * num_cells_x,
size_t * num_cells_y);
43int main (
int argc,
char *argv[]) {
52 printf (
". main: calling yac_cinit\n");
55 double tic, toc, time;
56 double time_min, time_max, time_ave;
57 double time_min_acc, time_max_acc, time_ave_acc;
59 MPI_Barrier(MPI_COMM_WORLD);
65#define NUM_CELLS_X (64)
66#define NUM_CELLS_Y (32)
70#define NUM_CELLS_X (96)
71#define NUM_CELLS_Y (48)
75#define NUM_CELLS_X (128)
76#define NUM_CELLS_Y (64)
80#define NUM_CELLS_X (192)
81#define NUM_CELLS_Y (96)
85#define NUM_CELLS_X (256)
86#define NUM_CELLS_Y (128)
90#define NUM_CELLS_X (320)
91#define NUM_CELLS_Y (160)
95#define NUM_CELLS_X (384)
96#define NUM_CELLS_Y (192)
100#define NUM_CELLS_X (480)
101#define NUM_CELLS_Y (240)
105#define NUM_CELLS_X (768)
106#define NUM_CELLS_Y (384)
110#define NUM_CELLS_X (256)
111#define NUM_CELLS_Y (128)
115 char const * configFilename =
"toy_atm_ocn.yaml";
118 parse_arguments(argc, argv, &configFilename, &num_cells_x, &num_cells_y);
125 printf (
". main: calling yac_cdef_comp\n");
129 char * comp_name =
"ATMOS";
133 printf (
". main: defined %s with local comp ID %i \n",
"toy-reg2d-atmosphere",
comp_id );
142 MPI_Comm_rank(local_comm,&rank);
143 MPI_Comm_size(local_comm,&size);
153 int total_nbr_cells[2];
155 total_nbr_cells[0] = (int)num_cells_x;
156 total_nbr_cells[1] = (int)num_cells_y;
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;
177 if (num_procs[1] > 1) {
181 if ((block_pos[1] != 0) && (block_pos[1] != num_procs[1]-1))
186 nbr_points[0] = nbr_cells[0];
188 nbr_points[0] = nbr_cells[0]+1;
190 nbr_points[1] = nbr_cells[1] + 1;
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));;
202 global_x_vertices[i] =
203 0.0 + (2.0 * M_PI / ((
double)
NUM_CELLS_X)) * (double)(i-1);
205 global_y_vertices[i] =
210 x_vertices = x_points = global_x_vertices + 1;
212 x_vertices = x_points = global_x_vertices + block_size[0] * block_pos[0];
214 y_vertices = y_points = global_y_vertices + block_size[1] * block_pos[1];
216 if (block_pos[1] != 0) {
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));
231 global_global_cell_id[j * (
NUM_CELLS_X + 2) + i] =
id;
232 global_global_cell_id_rank[j * (
NUM_CELLS_X + 2) + i] =
244 if (num_procs[0] == 1) {
246 for (
int j = 0,
id = 0; j <
NUM_CELLS_Y + 1; ++j) {
249 global_global_corner_id[j * (
NUM_CELLS_X + 3) +
i] =
id++;
254 for (
int j = 0,
id = 0; j <
NUM_CELLS_Y + 1; ++j) {
257 global_global_corner_id[j * (
NUM_CELLS_X + 3) +
i] =
id++;
268 if(global_global_cell_id_rank[i] == rank) global_global_cell_id_rank[
i] = -1;
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));
275 int offset_x, offset_y;
280 offset_x = block_size[0] * block_pos[0];
281 offset_y = block_size[1] * block_pos[1] - (block_pos[1] != 0);
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];
290 free(global_global_cell_id_rank);
291 free(global_global_cell_id);
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];
297 free(global_global_corner_id);
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;
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];
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];
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];
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));
320 for (
int i = 0;
i < nbr_cells[1]*nbr_cells[0]; ++
i)
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;
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);
347 time_ave /= (double) size;
350 printf (
"toy-reg2d-atmosphere: Time for initialisation %f %f %f \n", time_min, time_max, time_ave );
355 printf (
". main: calling yac_cenddef\n");
358 MPI_Barrier(MPI_COMM_WORLD);
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);
371 time_ave /= (double) size;
374 printf (
"toy-icon-atmosphere: Time for search %f %f %f \n", time_min, time_max, time_ave );
376 unsigned num_cells = nbr_cells[0] * nbr_cells[1];
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));
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));
397 unsigned x_index, y_index;
399 y_index =
i / nbr_cells[0];
400 x_index =
i - y_index * nbr_cells[0];
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;
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;
416 cell_corners[
i*4+1] = y_index * nbr_cells[0];
417 cell_corners[
i*4+2] = (y_index + 1) * nbr_cells[0];
419 cell_corners[
i*4+3] = (y_index + 1) * nbr_cells[0] + x_index;
422 num_points_per_cell[
i] = 4;
427 double curr_x, curr_y;
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;
439 cell_out_conserv_data[
i] = cell_data_field[
i];
440 cell_out_hcsbb_data[
i] = cell_data_field[
i];
442 cell_in_conserv_data[
i] = -999;
443 cell_in_hcsbb_data[
i] = -999;
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) {
450 LLtoXYZ(x_vertices[j], y_vertices[i], &point_data[3*(i * nbr_points[0] + j)]);
458 for (
int step = 0; step <
num_steps; ++step) {
460 MPI_Barrier(MPI_COMM_WORLD);
465 double *point_set_data[1];
466 double **collection_data[1] = {point_set_data};
468 point_set_data[0] = cell_out_conserv_data;
469 yac_cput(field_ids[0], 1, collection_data, &info, &err);
472 point_set_data[0] = cell_out_hcsbb_data;
473 yac_cput(field_ids[1], 1, collection_data, &info, &err);
478 double *collection_data[1];
480 collection_data[0] = cell_in_conserv_data;
481 yac_cget(field_ids[2], 1, collection_data, &info, &err);
484 collection_data[0] = cell_in_hcsbb_data;
485 yac_cget(field_ids[3], 1, collection_data, &info, &err);
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);
496 time_ave /= (
double) size;
499 time_min_acc += time_min;
500 time_max_acc += time_max;
501 time_ave_acc += time_ave;
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]);
518 double f1err_1 = 0.0;
519 double f1err_100 = 0.0;
523 double f2err_1 = 0.0;
524 double f2err_100 = 0.0;
530 if ( cell_in_conserv_err_rel[i] < 100.0 ) {
531 f1err_1 += cell_in_conserv_err_rel[
i];
534 f1err_100 += cell_in_conserv_err_rel[
i];
538 if ( cell_in_hcsbb_err_rel[i] < 100.0 ) {
539 f2err_1 += cell_in_hcsbb_err_rel[
i];
542 f2err_100 += cell_in_hcsbb_err_rel[
i];
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;
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;
553 printf (
"avg. rel err. for 1st field %f\n", f1err_1);
554 printf (
"avg. rel err. for 2nd field %f\n", f2err_1);
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);
564 char vtk_filename[36];
566 sprintf(vtk_filename,
"toy-reg2d-atmosphere_out_%d_%d.vtk", rank, step);
572 vtk_file, point_data, nbr_points[0]*nbr_points[1]);
575 vtk_file, cell_corners, num_points_per_cell,
num_cells);
578 vtk_file, corner_core_mask, nbr_points[0]*nbr_points[1],
581 vtk_file, local_global_corner_id, nbr_points[0]*nbr_points[1],
586 vtk_file, local_global_cell_id,
num_cells,
"global_cell_id");
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");
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;
619 printf (
"toy-icon-atmosphere: Time for ping-pong %f %f %f \n", time_min_acc, time_max_acc, time_ave_acc );
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);
636 free(num_points_per_cell);
639 free(corner_core_mask);
640 free(local_global_corner_id_rank);
641 free(local_global_corner_id);
643 free(local_global_cell_id_rank);
644 free(local_global_cell_id);
645 free(global_x_vertices);
646 free(global_y_vertices);
652 int argc,
char ** argv,
653 char const ** configFilename,
size_t * num_cells_x,
size_t * num_cells_y) {
656 while ((opt = getopt(argc, argv,
"c:x:y:")) != -1) {
657 YAC_ASSERT((opt ==
'c') || (opt ==
'x') || (opt ==
'y'),
"invalid command argument")
661 *configFilename = optarg;
664 *num_cells_x = atoi(optarg);
668 *num_cells_y = atoi(optarg);
#define YAC_ASSERT(exp, msg)
void yac_generate_reg2d_decomp(int num_points[2], int total_num_procs, int *num_procs)
double yac_test_func(double lon, double lat)
#define YAC_ASSERT_ARGS(exp, msg)
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[])
void yac_vtk_write_cell_scalars_double(YAC_VTK_FILE *vtk_file, double *scalars, unsigned num_cells, char *name)
void yac_vtk_write_cell_scalars_int(YAC_VTK_FILE *vtk_file, int *scalars, unsigned num_cells, char *name)
void yac_vtk_write_point_data(YAC_VTK_FILE *vtk_file, double *point_data, unsigned num_points)
void yac_vtk_write_point_scalars_int(YAC_VTK_FILE *vtk_file, int *scalars, unsigned num_points, char *name)
void yac_vtk_close(YAC_VTK_FILE *vtk_file)
void yac_vtk_write_cell_data(YAC_VTK_FILE *vtk_file, unsigned *cell_corners, unsigned *num_points_per_cell, unsigned num_cells)
YAC_VTK_FILE * yac_vtk_open(const char *filename, const char *title)
void yac_cset_global_index(int const *global_index, int location, int grid_id)
int const YAC_LOCATION_CELL
void yac_cfinalize()
Finalises YAC.
int const YAC_LOCATION_CORNER
void yac_cget_comp_comm(int comp_id, MPI_Comm *comp_comm)
void yac_cget(int const field_id, int collection_size, double **recv_field, int *info, int *ierr)
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)
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)
int const YAC_TIME_UNIT_SECOND
void yac_cput(int const field_id, int const collection_size, double ***const send_field, int *info, int *ierr)
void yac_cread_config_yaml(const char *yaml_filename)
void yac_cset_core_mask(int const *is_core, int location, int grid_id)
void yac_cdef_comp(char const *comp_name, int *comp_id)
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)