YAC 3.8.0
Yet Another Coupler
Loading...
Searching...
No Matches
read_icon_grid.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#ifdef HAVE_CONFIG_H
6#include "config.h"
7#endif
8
9#include <stdlib.h>
10#include <stdio.h>
11#include <string.h>
12#include <assert.h>
13#include <math.h>
14
15#include "read_icon_grid.h"
16#include "utils_common.h"
17#include "io_utils.h"
18#include "geometry.h"
19#include "read_grid.h"
20
21#ifdef YAC_NETCDF_ENABLED
22#include <netcdf.h>
23
24static int * get_icon_cell_mask( int ncid, size_t nbr_cells );
25static void get_icon_connect(
26 int ncid, size_t nbr_cells, int ** vertex_of_cell, size_t * nv);
27
29 const char * filename, int * nbr_vertices, int * nbr_cells,
30 int ** num_vertices_per_cell, int ** cell_to_vertex,
31 double ** x_vertices, double ** y_vertices,
32 double ** x_cells, double ** y_cells, int ** global_cell_id,
33 int ** cell_mask, int ** cell_core_mask,
34 int ** global_corner_id, int ** corner_core_mask, int rank, int size) {
35
36 // read the global grid
38 filename, nbr_vertices, nbr_cells, num_vertices_per_cell, cell_to_vertex,
39 x_vertices, y_vertices, x_cells, y_cells, cell_mask);
40
41 // determine local range
42 int local_start =
43 ((unsigned long)(*nbr_cells) * (unsigned long)rank) /
44 (unsigned long)size;
45 int num_local_cells =
46 ((unsigned long)(*nbr_cells) * ((unsigned long)rank+1)) /
47 (unsigned long)size - (unsigned long)local_start;
48
49 // mask for required vertices and cells
50 int * required_vertices = xcalloc(*nbr_vertices, sizeof(*required_vertices));
51 int * required_cells = xcalloc(*nbr_cells, sizeof(*required_cells));
52
53 int offset = 0;
54 for (int i = 0; i < local_start; ++i)
55 offset += (*num_vertices_per_cell)[i];
56
57 // mark all local cells and their vertices as required
58 for (int i = local_start; i < local_start + num_local_cells; ++i) {
59
60 required_cells[i] = 2;
61 for (int j = 0; j < (*num_vertices_per_cell)[i]; ++j)
62 required_vertices[(*cell_to_vertex)[offset+j]] = 2;
63
64 offset += (*num_vertices_per_cell)[i];
65 }
66
67 // mark all halo cells as required
68 offset = 0;
69 for (int i = 0; i < *nbr_cells; ++i) {
70
71 if (!required_cells[i]) {
72
73 for (int j = 0; j < (*num_vertices_per_cell)[i]; ++j) {
74
75 if (required_vertices[(*cell_to_vertex)[offset+j]]) {
76
77 required_cells[i] = 1;
78 break;
79 }
80 }
81
82 }
83
84 offset += (*num_vertices_per_cell)[i];
85 }
86
87 // mark all halo vertices as required
88 offset = 0;
89 for (int i = 0; i < *nbr_cells; ++i) {
90
91 if (required_cells[i] == 1) {
92
93 for (int j = 0; j < (*num_vertices_per_cell)[i]; ++j) {
94
95 if (!required_vertices[(*cell_to_vertex)[offset+j]]) {
96
97 required_vertices[(*cell_to_vertex)[offset+j]] = 1;
98 }
99 }
100 }
101
102 offset += (*num_vertices_per_cell)[i];
103 }
104
105 // count the number of cells and vertices
106 int part_num_vertices = 0;
107 int part_num_cells = 0;
108 for (int i = 0; i < *nbr_vertices; ++i)
109 if (required_vertices[i])
110 part_num_vertices++;
111 for (int i = 0; i < *nbr_cells; ++i)
112 if(required_cells[i])
113 part_num_cells++;
114
115 *global_cell_id = xmalloc(part_num_cells * sizeof(**global_cell_id));
116 *cell_core_mask = xmalloc(part_num_cells * sizeof(**cell_core_mask));
117 *global_corner_id = xmalloc(part_num_vertices * sizeof(**global_corner_id));
118 *corner_core_mask = xmalloc(part_num_vertices * sizeof(**corner_core_mask));
119
120 // generate final vertex data
121 part_num_vertices = 0;
122 int * global_to_local_vertex =
123 xmalloc(*nbr_vertices * sizeof(*global_to_local_vertex));
124 for (int i = 0; i < *nbr_vertices; ++i) {
125
126 if (required_vertices[i]) {
127
128 (*global_corner_id)[part_num_vertices] = i;
129 (*corner_core_mask)[part_num_vertices] = required_vertices[i] == 2;
130 (*x_vertices)[part_num_vertices] = (*x_vertices)[i];
131 (*y_vertices)[part_num_vertices] = (*y_vertices)[i];
132 global_to_local_vertex[i] = part_num_vertices;
133 part_num_vertices++;
134 }
135 }
136
137 *x_vertices = xrealloc(*x_vertices, part_num_vertices * sizeof(**x_vertices));
138 *y_vertices = xrealloc(*y_vertices, part_num_vertices * sizeof(**y_vertices));
139 *nbr_vertices = part_num_vertices;
140 free(required_vertices);
141
142 // generate final cell data
143 int num_cell_vertex_dependencies = 0;
144 part_num_cells = 0;
145 offset = 0;
146 for (int i = 0; i < *nbr_cells; ++i) {
147
148 if (required_cells[i]) {
149
150 (*global_cell_id)[part_num_cells] = i;
151 (*cell_core_mask)[part_num_cells] = required_cells[i] == 2;
152 (*x_cells)[part_num_cells] = (*x_cells)[i];
153 (*y_cells)[part_num_cells] = (*y_cells)[i];
154 (*cell_mask)[part_num_cells] = (*cell_mask)[i];
155
156 for (int j = 0; j < (*num_vertices_per_cell)[i]; ++j)
157 (*cell_to_vertex)[num_cell_vertex_dependencies++] =
158 global_to_local_vertex[(*cell_to_vertex)[offset+j]];
159
160 (*num_vertices_per_cell)[part_num_cells] = (*num_vertices_per_cell)[i];
161
162 part_num_cells++;
163 }
164
165 offset += (*num_vertices_per_cell)[i];
166 }
167
168 *x_cells = xrealloc(*x_cells, part_num_cells * sizeof(**x_cells));
169 *y_cells = xrealloc(*y_cells, part_num_cells * sizeof(**y_cells));
170 *cell_mask = xrealloc(*cell_mask, part_num_cells * sizeof(**cell_mask));
171
172 *num_vertices_per_cell = xrealloc(*num_vertices_per_cell, part_num_cells *
173 sizeof(**num_vertices_per_cell));
174 *cell_to_vertex = xrealloc(*cell_to_vertex, num_cell_vertex_dependencies *
175 sizeof(**cell_to_vertex));
176 *nbr_cells = part_num_cells;
177 free(required_cells);
178 free(global_to_local_vertex);
179}
180
181void yac_read_icon_grid_information(const char * filename, int * nbr_vertices,
182 int * nbr_cells, int ** num_vertices_per_cell,
183 int ** cell_to_vertex, double ** x_vertices,
184 double ** y_vertices, double ** x_cells,
185 double ** y_cells, int ** cell_mask) {
186
187 /* Open file */
188 int ncid;
189 yac_nc_open(filename, NC_NOWRITE, &ncid);
190
191 /* Get vertex longitudes and latitudes of cells */
192 size_t nbr_vertices_;
194 ncid, "vlon", "vlat", x_vertices, y_vertices, &nbr_vertices_);
195 *nbr_vertices = (int)nbr_vertices_;
196
197 /* Get cell center longitudes and latitudes of cells */
198 size_t nbr_cells_;
199 yac_read_coords(ncid, "clon", "clat", x_cells, y_cells, &nbr_cells_);
200 *nbr_cells = (int)nbr_cells_;
201
202 /* Get mask of cells */
203 *cell_mask = get_icon_cell_mask ( ncid, *nbr_cells );
204
205 /* Get relations between vertices and cells */
206 int * vertex_of_cell;
207 size_t nv;
208 get_icon_connect(ncid, *nbr_cells, &vertex_of_cell, &nv);
209
210 /* Close file */
211 YAC_HANDLE_ERROR(nc_close(ncid));
212
213 //-------------------------------------------------------------------------//
214
215 *num_vertices_per_cell = xmalloc(*nbr_cells * sizeof(**num_vertices_per_cell));
216 for (int i = 0; i < *nbr_cells; (*num_vertices_per_cell)[i++] = (int)nv);
217
218 *cell_to_vertex = xmalloc(*nbr_cells * nv * sizeof(**cell_to_vertex));
219
220 // Unfortunately the data is only available in Fortran order
221 for (int i = 0; i < *nbr_cells; ++i) {
222
223 for (size_t j = 0; j < nv; ++j)
224 (*cell_to_vertex)[nv * i + j] =
225 vertex_of_cell[i + j * (*nbr_cells)] - 1;
226 }
227
228 free(vertex_of_cell);
229}
230
231// taken from scales-ppm library
232// https://www.dkrz.de/redmine/projects/scales-ppm
233static inline int
234partition_idx_from_element_idx(unsigned element_idx, unsigned num_elements,
235 int num_partitions) {
236
237 return (int)((((unsigned long)element_idx) * ((unsigned long)num_partitions) +
238 (unsigned long)num_partitions - 1) /
239 ((unsigned long)num_elements));
240}
241
242static void convert_to_rad(int ncid, int varid, double * array, size_t count) {
243
244 if (yac_check_coord_units(ncid, varid))
245 for (size_t i = 0; i < count; ++i) array[i] *= YAC_RAD;
246}
247
249 const char * filename, MPI_Comm comm, int * nbr_vertices, int * nbr_cells,
250 int ** num_vertices_per_cell, int ** cell_to_vertex, int ** global_cell_ids,
251 int ** cell_owner, int ** global_vertex_ids, int ** vertex_owner,
252 double ** x_vertices, double ** y_vertices,
253 double ** x_cells, double ** y_cells, int ** cell_msk) {
254
256 ((x_cells == NULL) && (y_cells == NULL)) ||
257 ((x_cells != NULL) && (y_cells != NULL)),
258 "ERROR(yac_read_icon_grid_information_parallel): "
259 "arguments x_cells and y_cells have to be either both NULL or "
260 "both different from NULL");
261
262 const int with_cell_coords = x_cells != NULL;
263 const int with_cell_mask = cell_msk != NULL;
264
265 int comm_rank, comm_size;
266
267 MPI_Comm_rank(comm, &comm_rank);
268 MPI_Comm_size(comm, &comm_size);
269
270 int local_is_io, * io_ranks, num_io_ranks;
271 yac_get_io_ranks(comm, &local_is_io, &io_ranks, &num_io_ranks);
272
273 size_t ncells, nvertices, nv, ne;
274
275 size_t read_local_start_cell = 0;
276 size_t read_num_local_cells = 0;
277 size_t read_local_start_vertex = 0;
278 size_t read_num_local_vertices = 0;
279 double * read_cell_lon = NULL;
280 double * read_cell_lat = NULL;
281 int * read_cell_mask = NULL;
282 double * read_vertex_lon = NULL;
283 double * read_vertex_lat = NULL;
284 int * read_dist_vertex_of_cell = NULL;
285 int * read_dist_cells_of_vertex = NULL;
286
287 if (local_is_io) {
288
289 unsigned long io_proc_idx = ULONG_MAX;
290 for (int i = 0; (i < num_io_ranks) && (io_proc_idx == ULONG_MAX); ++i)
291 if (io_ranks[i] == comm_rank)
292 io_proc_idx = (unsigned long)i;
293
294 // open file
295 int ncid;
296 yac_nc_open(filename, NC_NOWRITE, &ncid);
297
298 // get number of cells and vertices
299 int dim_id;
300 yac_nc_inq_dimid(ncid, "cell", &dim_id);
301 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dim_id, &ncells));
302 yac_nc_inq_dimid(ncid, "vertex", &dim_id);
303 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dim_id, &nvertices));
304 yac_nc_inq_dimid(ncid, "nv", &dim_id);
305 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dim_id, &nv));
306 yac_nc_inq_dimid(ncid, "ne", &dim_id);
307 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dim_id, &ne));
308
309 // determine local range for cell and vertex data
310 read_local_start_cell =
311 ((unsigned long)ncells * io_proc_idx) / (unsigned long)num_io_ranks;
312 read_num_local_cells =
313 ((unsigned long)ncells * (io_proc_idx+1)) / (unsigned long)num_io_ranks -
314 (unsigned long)read_local_start_cell;
315 read_local_start_vertex =
316 ((unsigned long)nvertices * io_proc_idx) / (unsigned long)num_io_ranks;
317 read_num_local_vertices =
318 ((unsigned long)nvertices * (io_proc_idx+1)) / (unsigned long)num_io_ranks -
319 (unsigned long)read_local_start_vertex;
320
321 // read basic grid data (each process its individual part)
322 read_cell_lon =
323 with_cell_coords?
324 xmalloc(read_num_local_cells * sizeof(*read_cell_lon)):NULL;
325 read_cell_lat =
326 with_cell_coords?
327 xmalloc(read_num_local_cells * sizeof(*read_cell_lat)):NULL;
328 read_cell_mask =
329 with_cell_mask?
330 xmalloc(read_num_local_cells * sizeof(*read_cell_mask)):NULL;
331 read_vertex_lon = xmalloc(read_num_local_vertices * sizeof(*read_vertex_lon));
332 read_vertex_lat = xmalloc(read_num_local_vertices * sizeof(*read_vertex_lat));
333 read_dist_vertex_of_cell =
334 xmalloc(read_num_local_cells * nv * sizeof(*read_dist_vertex_of_cell));
335 read_dist_cells_of_vertex =
336 xmalloc(read_num_local_vertices * ne * sizeof(*read_dist_cells_of_vertex));
337 int varid;
338 if (with_cell_coords) {
339 yac_nc_inq_varid(ncid, "clon", &varid);
340 YAC_HANDLE_ERROR(nc_get_vara_double(ncid, varid, &read_local_start_cell,
341 &read_num_local_cells, read_cell_lon));
342 convert_to_rad(ncid, varid, read_cell_lon, read_num_local_cells);
343 yac_nc_inq_varid(ncid, "clat", &varid);
344 YAC_HANDLE_ERROR(nc_get_vara_double(ncid, varid, &read_local_start_cell,
345 &read_num_local_cells, read_cell_lat));
346 convert_to_rad(ncid, varid, read_cell_lat, read_num_local_cells);
347 }
348 if (with_cell_mask) {
349 yac_nc_inq_varid(ncid, "cell_sea_land_mask", &varid);
350 YAC_HANDLE_ERROR(nc_get_vara_int(ncid, varid, &read_local_start_cell,
351 &read_num_local_cells, read_cell_mask));
352 }
353 yac_nc_inq_varid(ncid, "vlon", &varid);
354 YAC_HANDLE_ERROR(nc_get_vara_double(ncid, varid, &read_local_start_vertex,
355 &read_num_local_vertices, read_vertex_lon));
356 convert_to_rad(ncid, varid, read_vertex_lon, read_num_local_vertices);
357 yac_nc_inq_varid(ncid, "vlat", &varid);
358 YAC_HANDLE_ERROR(nc_get_vara_double(ncid, varid, &read_local_start_vertex,
359 &read_num_local_vertices, read_vertex_lat));
360 convert_to_rad(ncid, varid, read_vertex_lat, read_num_local_vertices);
361 {
362 int * buffer = xmalloc(MAX(nv*read_num_local_cells,
363 ne*read_num_local_vertices) * sizeof(*buffer));
364 {
365 size_t tmp_start[2] = {0, read_local_start_cell};
366 size_t tmp_count[2] = {nv, read_num_local_cells};
367 yac_nc_inq_varid(ncid, "vertex_of_cell", &varid);
368 YAC_HANDLE_ERROR(nc_get_vara_int(ncid, varid, tmp_start, tmp_count, buffer));
369 for (size_t i = 0; i < read_num_local_cells; ++i)
370 for (size_t j = 0; j < nv; ++j)
371 read_dist_vertex_of_cell[i * nv + j] =
372 buffer[i + j * read_num_local_cells];
373 }
374 {
375 size_t tmp_start[2] = {0, read_local_start_vertex};
376 size_t tmp_count[2] = {ne, read_num_local_vertices};
377 yac_nc_inq_varid(ncid, "cells_of_vertex", &varid);
378 YAC_HANDLE_ERROR(nc_get_vara_int(ncid, varid, tmp_start, tmp_count, buffer));
379 for (size_t i = 0; i < read_num_local_vertices; ++i)
380 for (size_t j = 0; j < ne; ++j)
381 read_dist_cells_of_vertex[i * ne + j] =
382 buffer[i + j * read_num_local_vertices];
383 }
384 free(buffer);
385
386 // adjust for c indices
387 for (size_t i = 0; i < read_num_local_cells * nv; ++i)
388 if (read_dist_vertex_of_cell[i] > 0) read_dist_vertex_of_cell[i]--;
389 for (size_t i = 0; i < read_num_local_vertices * ne; ++i)
390 read_dist_cells_of_vertex[i]--;
391 }
392
393 YAC_HANDLE_ERROR(nc_close(ncid));
394 } else {
395 read_cell_lon =
396 with_cell_coords?xmalloc(1 * sizeof(*read_cell_lon)):NULL;
397 read_cell_lat =
398 with_cell_coords?xmalloc(1 * sizeof(*read_cell_lat)):NULL;
399 read_cell_mask =
400 with_cell_mask?xmalloc(1 * sizeof(*read_cell_mask)):NULL;
401 read_vertex_lon = xmalloc(1 * sizeof(*read_vertex_lon));
402 read_vertex_lat = xmalloc(1 * sizeof(*read_vertex_lat));
403 read_dist_vertex_of_cell = xmalloc(1 * sizeof(*read_dist_vertex_of_cell));
404 read_dist_cells_of_vertex = xmalloc(1 * sizeof(*read_dist_cells_of_vertex));
405 }
406
407 free(io_ranks);
408
409 {
410 int tmp;
411 if (comm_rank == 0) tmp = (int)ncells;
412 MPI_Bcast(&tmp, 1, MPI_INT, 0, comm);
413 ncells = (size_t)tmp;
414 if (comm_rank == 0) tmp = (int)nvertices;
415 MPI_Bcast(&tmp, 1, MPI_INT, 0, comm);
416 nvertices = (size_t)tmp;
417 if (comm_rank == 0) tmp = (int)nv;
418 MPI_Bcast(&tmp, 1, MPI_INT, 0, comm);
419 nv = (size_t)tmp;
420 if (comm_rank == 0) tmp = (int)ne;
421 MPI_Bcast(&tmp, 1, MPI_INT, 0, comm);
422 ne = (size_t)tmp;
423 }
424
425 // determine local range for cell and vertex data
426 size_t local_start_cell =
427 ((unsigned long)ncells * (unsigned long)comm_rank) /
428 (unsigned long)comm_size;
429 size_t num_local_cells =
430 ((unsigned long)ncells * ((unsigned long)comm_rank+1)) /
431 (unsigned long)comm_size - (unsigned long)local_start_cell;
432 size_t local_start_vertex =
433 ((unsigned long)nvertices * (unsigned long)comm_rank) /
434 (unsigned long)comm_size;
435 size_t num_local_vertices =
436 ((unsigned long)nvertices * ((unsigned long)comm_rank+1)) /
437 (unsigned long)comm_size - (unsigned long)local_start_vertex;
438
439 // redistribute basic cell data (from io decomposition)
440 double * cell_lon =
441 with_cell_coords?xmalloc(num_local_cells * sizeof(*cell_lon)):NULL;
442 double * cell_lat =
443 with_cell_coords?xmalloc(num_local_cells * sizeof(*cell_lat)):NULL;
444 int * cell_mask =
445 with_cell_mask?xmalloc(num_local_cells * sizeof(*cell_mask)):NULL;
446 int * dist_vertex_of_cell =
447 xmalloc(num_local_cells * nv * sizeof(*dist_vertex_of_cell));
448 {
449 int * send_count = xcalloc(comm_size, sizeof(*send_count));
450 int * recv_count = xcalloc(comm_size, sizeof(*recv_count));
451
452 for (size_t i = 0; i < read_num_local_cells; ++i)
453 send_count[
455 read_local_start_cell + i, ncells, comm_size)]++;
456
457 MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
458
459 int * send_displ = xmalloc(comm_size * sizeof(*send_displ));
460 int * recv_displ = xmalloc(comm_size * sizeof(*recv_displ));
461 int send_accum = 0, recv_accum = 0;
462 for (int i = 0; i < comm_size; ++i) {
463 send_displ[i] = send_accum;
464 recv_displ[i] = recv_accum;
465 send_accum += send_count[i];
466 recv_accum += recv_count[i];
467 }
468
469 if (with_cell_coords) {
470 MPI_Alltoallv(read_cell_lon, send_count, send_displ, MPI_DOUBLE,
471 cell_lon, recv_count, recv_displ, MPI_DOUBLE, comm);
472 MPI_Alltoallv(read_cell_lat, send_count, send_displ, MPI_DOUBLE,
473 cell_lat, recv_count, recv_displ, MPI_DOUBLE, comm);
474 }
475 if (with_cell_mask) {
476 MPI_Alltoallv(read_cell_mask, send_count, send_displ, MPI_INT,
477 cell_mask, recv_count, recv_displ, MPI_INT, comm);
478 }
479
480 for (int i = 0; i < comm_size; ++i) {
481 send_count[i] *= nv;
482 send_displ[i] *= nv;
483 recv_count[i] *= nv;
484 recv_displ[i] *= nv;
485 }
486
487 MPI_Alltoallv(read_dist_vertex_of_cell, send_count, send_displ, MPI_INT,
488 dist_vertex_of_cell, recv_count, recv_displ, MPI_INT, comm);
489
490 free(recv_displ);
491 free(send_displ);
492 free(recv_count);
493 free(send_count);
494 free(read_cell_lon);
495 free(read_cell_lat);
496 free(read_cell_mask);
497 free(read_dist_vertex_of_cell);
498 }
499
500 // redistribute basic vertex data (from io decomposition)
501 double * vertex_lon = xmalloc(num_local_vertices * sizeof(*vertex_lon));
502 double * vertex_lat = xmalloc(num_local_vertices * sizeof(*vertex_lat));
503 int * dist_cells_of_vertex =
504 xmalloc(num_local_vertices * ne * sizeof(*dist_cells_of_vertex));
505 {
506 int * send_count = xcalloc(comm_size, sizeof(*send_count));
507 int * recv_count = xcalloc(comm_size, sizeof(*recv_count));
508
509 for (size_t i = 0; i < read_num_local_vertices; ++i)
510 send_count[
512 read_local_start_vertex + i, nvertices, comm_size)]++;
513
514 MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
515
516 int * send_displ = xmalloc(comm_size * sizeof(*send_displ));
517 int * recv_displ = xmalloc(comm_size * sizeof(*recv_displ));
518 int send_accum = 0, recv_accum = 0;
519 for (int i = 0; i < comm_size; ++i) {
520 send_displ[i] = send_accum;
521 recv_displ[i] = recv_accum;
522 send_accum += send_count[i];
523 recv_accum += recv_count[i];
524 }
525
526 MPI_Alltoallv(read_vertex_lon, send_count, send_displ, MPI_DOUBLE,
527 vertex_lon, recv_count, recv_displ, MPI_DOUBLE, comm);
528 MPI_Alltoallv(read_vertex_lat, send_count, send_displ, MPI_DOUBLE,
529 vertex_lat, recv_count, recv_displ, MPI_DOUBLE, comm);
530
531 for (int i = 0; i < comm_size; ++i) {
532 send_count[i] *= ne;
533 send_displ[i] *= ne;
534 recv_count[i] *= ne;
535 recv_displ[i] *= ne;
536 }
537
538 MPI_Alltoallv(read_dist_cells_of_vertex, send_count, send_displ, MPI_INT,
539 dist_cells_of_vertex, recv_count, recv_displ, MPI_INT, comm);
540
541 free(recv_displ);
542 free(send_displ);
543 free(recv_count);
544 free(send_count);
545 free(read_vertex_lon);
546 free(read_vertex_lat);
547 free(read_dist_cells_of_vertex);
548 }
549
550 // determine required vertices for core cells
551 size_t num_core_vertices = num_local_cells * nv;
552 int * core_vertices = xmalloc(num_core_vertices * sizeof(*core_vertices));
553 {
554 memcpy(core_vertices, dist_vertex_of_cell,
555 num_core_vertices * sizeof(*core_vertices));
556 yac_quicksort_index(core_vertices, num_core_vertices, NULL);
557 yac_remove_duplicates_int(core_vertices, &num_core_vertices);
558 core_vertices =
559 xrealloc(core_vertices, num_core_vertices * sizeof(*core_vertices));
560 }
561
562 // get cells_of_vertex for core vertices and compute dist_vertex_owner
563 int * cells_of_vertex_core =
564 xmalloc(num_core_vertices * ne * sizeof(*cells_of_vertex_core));
565 int * dist_vertex_owner =
566 xmalloc(num_local_vertices * sizeof(*dist_vertex_owner));
567 {
568 int * send_count = xcalloc(comm_size, sizeof(*send_count));
569 int * recv_count = xcalloc(comm_size, sizeof(*recv_count));
570
571 for (size_t i = 0; i < num_core_vertices; ++i)
572 send_count[((unsigned long)(core_vertices[i]) *
573 (unsigned long)comm_size + (unsigned long)comm_size - 1) /
574 (unsigned long)nvertices]++;
575
576 MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
577
578 int * send_displ = xmalloc(comm_size * sizeof(*send_displ));
579 int * recv_displ = xmalloc(comm_size * sizeof(*recv_displ));
580 int send_accum = 0, recv_accum = 0;
581 for (int i = 0; i < comm_size; ++i) {
582 send_displ[i] = send_accum;
583 recv_displ[i] = recv_accum;
584 send_accum += send_count[i];
585 recv_accum += recv_count[i];
586 }
587
588 int num_core_vertices_remote = 0;
589 for (int i = 0; i < comm_size; ++i)
590 num_core_vertices_remote += recv_count[i];
591
592 int * remote_vertex_buffer =
593 xmalloc(num_core_vertices_remote * sizeof(*remote_vertex_buffer));
594
595 MPI_Alltoallv(core_vertices, send_count, send_displ, MPI_INT,
596 remote_vertex_buffer, recv_count, recv_displ, MPI_INT, comm);
597
598 for (size_t i = 0; i < num_local_vertices; ++i) dist_vertex_owner[i] = -1;
599
600 for (int i = 0, j = 0; i < comm_size; ++i)
601 for (int k = 0; k < recv_count[i]; ++k, ++j)
602 dist_vertex_owner[remote_vertex_buffer[j] - local_start_vertex] = i;
603
604 int * send_cell_of_vertex =
605 xmalloc(num_core_vertices_remote * ne * sizeof(*send_cell_of_vertex));
606
607 for (int i = 0, l = 0, m = 0; i < comm_size; ++i) {
608 for (int j = 0; j < recv_count[i]; ++j, ++l) {
609 int idx = remote_vertex_buffer[l] - local_start_vertex;
610 for (size_t k = 0; k < ne; ++k, ++m)
611 send_cell_of_vertex[m] = dist_cells_of_vertex[idx * ne + k];
612 }
613 }
614 free(dist_cells_of_vertex);
615
616 for (int i = 0; i < comm_size; ++i) {
617 send_count[i] *= ne;
618 send_displ[i] *= ne;
619 recv_count[i] *= ne;
620 recv_displ[i] *= ne;
621 }
622
623 MPI_Alltoallv(send_cell_of_vertex, recv_count, recv_displ, MPI_INT,
624 cells_of_vertex_core, send_count, send_displ, MPI_INT, comm);
625
626 free(send_cell_of_vertex);
627 free(remote_vertex_buffer);
628 free(recv_displ);
629 free(send_displ);
630 free(recv_count);
631 free(send_count);
632 }
633
634 // determine halo cells
635 int * halo_cells;
636 size_t num_halo_cells = 0;
637 {
638 int * tmp_cells = cells_of_vertex_core;
639 size_t num_tmp_cells = num_core_vertices * ne;
640
641 yac_quicksort_index(tmp_cells, num_tmp_cells, NULL);
642 yac_remove_duplicates_int(tmp_cells, &num_tmp_cells);
643
644 if ((num_tmp_cells > 0) && (tmp_cells[0] == -1)) {
645 num_tmp_cells--;
646 tmp_cells++;
647 }
648
649 halo_cells = xmalloc((num_tmp_cells - num_local_cells) * sizeof(*halo_cells));
650
651 size_t i = 0;
652 for (; i < num_tmp_cells && tmp_cells[i] < (int)local_start_cell; ++i)
653 halo_cells[num_halo_cells++] = tmp_cells[i];
654 i += num_local_cells;
655 for (; i < num_tmp_cells; ++i) halo_cells[num_halo_cells++] = tmp_cells[i];
656
657 assert(num_halo_cells == num_tmp_cells - num_local_cells);
658
659 free(cells_of_vertex_core);
660 }
661
662 // determine all vertices and get coordinates of halo cells
663 size_t num_all_local_vertices = num_halo_cells * nv + num_core_vertices;
664 int * all_cell_to_vertex;
665 int * all_local_vertices = xrealloc(core_vertices, num_all_local_vertices *
666 sizeof(*all_local_vertices));
667 if (with_cell_coords) {
668 cell_lon = xrealloc(cell_lon, (num_local_cells + num_halo_cells) *
669 sizeof(*cell_lon));
670 cell_lat = xrealloc(cell_lat, (num_local_cells + num_halo_cells) *
671 sizeof(*cell_lat));
672 }
673 if (with_cell_mask) {
674 cell_mask = xrealloc(cell_mask, (num_local_cells + num_halo_cells) *
675 sizeof(*cell_mask));
676 }
677 {
678 int * send_count = xcalloc(comm_size, sizeof(*send_count));
679 int * recv_count = xcalloc(comm_size, sizeof(*recv_count));
680
681 for (size_t i = 0; i < num_halo_cells; ++i)
682 send_count[((unsigned long)(halo_cells[i]) *
683 (unsigned long)comm_size + (unsigned long)comm_size - 1) /
684 (unsigned long)ncells]++;
685
686 MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
687
688 int * send_displ = xmalloc(comm_size * sizeof(*send_displ));
689 int * recv_displ = xmalloc(comm_size * sizeof(*recv_displ));
690 int send_accum = 0, recv_accum = 0;
691 for (int i = 0; i < comm_size; ++i) {
692 send_displ[i] = send_accum;
693 recv_displ[i] = recv_accum;
694 send_accum += send_count[i];
695 recv_accum += recv_count[i];
696 }
697
698 int num_halo_cells_remote = 0;
699 for (int i = 0; i < comm_size; ++i)
700 num_halo_cells_remote += recv_count[i];
701
702 int * remote_halo_cell_buffer =
703 xmalloc(num_halo_cells_remote * sizeof(*remote_halo_cell_buffer));
704
705 MPI_Alltoallv(halo_cells, send_count, send_displ, MPI_INT,
706 remote_halo_cell_buffer, recv_count, recv_displ, MPI_INT,
707 comm);
708
709 int * send_halo_cell_vertices =
710 xmalloc(num_halo_cells_remote * nv * sizeof(*send_halo_cell_vertices));
711 double * send_cell_lon =
712 with_cell_coords?
713 xmalloc(num_halo_cells_remote * sizeof(*send_cell_lon)):NULL;
714 double * send_cell_lat =
715 with_cell_coords?
716 xmalloc(num_halo_cells_remote * sizeof(*send_cell_lat)):NULL;
717 int * send_cell_mask =
718 with_cell_mask?
719 xmalloc(num_halo_cells_remote * sizeof(*send_cell_mask)):NULL;
720
721 for (int i = 0, l = 0, m = 0; i < comm_size; ++i) {
722 for (int j = 0; j < recv_count[i]; ++j, ++l) {
723 int idx = remote_halo_cell_buffer[l] - local_start_cell;
724 for (size_t k = 0; k < nv; ++k, ++m)
725 send_halo_cell_vertices[m] = dist_vertex_of_cell[idx * nv + k];
726 }
727 }
728 if (with_cell_coords) {
729 for (int i = 0, l = 0; i < comm_size; ++i) {
730 for (int j = 0; j < recv_count[i]; ++j, ++l) {
731 int idx = remote_halo_cell_buffer[l] - local_start_cell;
732 send_cell_lon[l] = cell_lon[idx];
733 send_cell_lat[l] = cell_lat[idx];
734 }
735 }
736 }
737 if (with_cell_mask) {
738 for (int i = 0, l = 0; i < comm_size; ++i) {
739 for (int j = 0; j < recv_count[i]; ++j, ++l) {
740 int idx = remote_halo_cell_buffer[l] - local_start_cell;
741 send_cell_mask[l] = cell_mask[idx];
742 }
743 }
744 }
745
746 if (with_cell_coords) {
747 MPI_Alltoallv(send_cell_lon, recv_count, recv_displ, MPI_DOUBLE,
748 cell_lon + num_local_cells, send_count, send_displ,
749 MPI_DOUBLE, comm);
750 MPI_Alltoallv(send_cell_lat, recv_count, recv_displ, MPI_DOUBLE,
751 cell_lat + num_local_cells, send_count, send_displ,
752 MPI_DOUBLE, comm);
753 }
754 if (with_cell_mask) {
755 MPI_Alltoallv(send_cell_mask, recv_count, recv_displ, MPI_INT,
756 cell_mask + num_local_cells, send_count, send_displ,
757 MPI_INT, comm);
758 }
759
760 for (int i = 0; i < comm_size; ++i) {
761 send_count[i] *= nv;
762 send_displ[i] *= nv;
763 recv_count[i] *= nv;
764 recv_displ[i] *= nv;
765 }
766
767 all_cell_to_vertex =
768 xrealloc(dist_vertex_of_cell, (num_local_cells + num_halo_cells) * nv *
769 sizeof(*all_cell_to_vertex));
770
771 MPI_Alltoallv(send_halo_cell_vertices, recv_count, recv_displ, MPI_INT,
772 all_cell_to_vertex + num_local_cells * nv, send_count,
773 send_displ, MPI_INT, comm);
774
775 memcpy(all_local_vertices + num_core_vertices,
776 all_cell_to_vertex + num_local_cells * nv,
777 num_halo_cells * nv * sizeof(*all_local_vertices));
778
779 free(send_cell_mask);
780 free(send_cell_lat);
781 free(send_cell_lon);
782 free(send_halo_cell_vertices);
783 free(remote_halo_cell_buffer);
784 free(recv_displ);
785 free(send_displ);
786 free(recv_count);
787 free(send_count);
788
790 all_local_vertices, num_all_local_vertices, NULL);
791 yac_remove_duplicates_int(all_local_vertices, &num_all_local_vertices);
792 }
793
794 // determine owner and coordinates for all vertices
795 double * all_vertex_lon =
796 xmalloc(num_all_local_vertices * sizeof(*all_vertex_lon));
797 double * all_vertex_lat =
798 xmalloc(num_all_local_vertices * sizeof(*all_vertex_lat));
799 int * all_local_vertices_owner =
800 xmalloc(num_all_local_vertices * sizeof(*all_local_vertices_owner));
801 {
802 int * send_count = xcalloc(comm_size, sizeof(*send_count));
803 int * recv_count = xcalloc(comm_size, sizeof(*recv_count));
804
805 for (size_t i = 0; i < num_all_local_vertices; ++i)
806 send_count[((unsigned long)(all_local_vertices[i]) *
807 (unsigned long)comm_size + (unsigned long)comm_size - 1) /
808 (unsigned long)nvertices]++;
809
810 MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
811
812 int * send_displ = xmalloc(comm_size * sizeof(*send_displ));
813 int * recv_displ = xmalloc(comm_size * sizeof(*recv_displ));
814 int send_accum = 0, recv_accum = 0;
815 for (int i = 0; i < comm_size; ++i) {
816 send_displ[i] = send_accum;
817 recv_displ[i] = recv_accum;
818 send_accum += send_count[i];
819 recv_accum += recv_count[i];
820 }
821
822 int num_all_local_vertices_remote = 0;
823 for (int i = 0; i < comm_size; ++i)
824 num_all_local_vertices_remote += recv_count[i];
825
826 int * remote_vertex_buffer =
827 xmalloc(num_all_local_vertices_remote * sizeof(*remote_vertex_buffer));
828
829 MPI_Alltoallv(all_local_vertices, send_count, send_displ, MPI_INT,
830 remote_vertex_buffer, recv_count, recv_displ, MPI_INT, comm);
831
832 int * send_vertex_owner = remote_vertex_buffer;
833 double * send_vertex_lon =
834 xmalloc(num_all_local_vertices_remote * sizeof(*send_vertex_lon));
835 double * send_vertex_lat =
836 xmalloc(num_all_local_vertices_remote * sizeof(*send_vertex_lat));
837
838 for (int i = 0, l = 0; i < comm_size; ++i) {
839 for (int j = 0; j < recv_count[i]; ++j, ++l) {
840 int idx = remote_vertex_buffer[l] - local_start_vertex;
841 send_vertex_owner[l] = dist_vertex_owner[idx];
842 send_vertex_lon[l] = vertex_lon[idx];
843 send_vertex_lat[l] = vertex_lat[idx];
844 }
845 }
846
847 free(dist_vertex_owner);
848 free(vertex_lon);
849 free(vertex_lat);
850
851 MPI_Alltoallv(send_vertex_owner, recv_count, recv_displ, MPI_INT,
852 all_local_vertices_owner, send_count, send_displ, MPI_INT,
853 comm);
854 MPI_Alltoallv(send_vertex_lon, recv_count, recv_displ, MPI_DOUBLE,
855 all_vertex_lon, send_count, send_displ, MPI_DOUBLE, comm);
856 MPI_Alltoallv(send_vertex_lat, recv_count, recv_displ, MPI_DOUBLE,
857 all_vertex_lat, send_count, send_displ, MPI_DOUBLE, comm);
858
859 free(send_vertex_lat);
860 free(send_vertex_lon);
861 free(send_vertex_owner);
862 free(recv_displ);
863 free(send_displ);
864 free(recv_count);
865 free(send_count);
866 }
867
868 // convert global ids within all_cell_to_vertex into local ids
869 {
870 size_t vertex_count = nv * (num_local_cells + num_halo_cells);
871 int * permutation = xmalloc(vertex_count * sizeof(*permutation));
872 for (size_t i = 0; i < vertex_count; ++i) permutation[i] = (int)i;
873
874 yac_quicksort_index(all_cell_to_vertex, vertex_count, permutation);
875
876 for (size_t i = 0, j = 0; i < vertex_count; ++i) {
877 while (all_local_vertices[j] != all_cell_to_vertex[i]) ++j;
878 all_cell_to_vertex[i] = (int)j;
879 }
880
881 yac_quicksort_index(permutation, vertex_count, all_cell_to_vertex);
882 free(permutation);
883 }
884
885 *nbr_vertices = num_all_local_vertices;
886 *nbr_cells = num_local_cells + num_halo_cells;
887
888 *num_vertices_per_cell =
889 xmalloc((num_local_cells + num_halo_cells) * sizeof(**num_vertices_per_cell));
890 for (size_t i = 0; i < num_local_cells + num_halo_cells; ++i)
891 (*num_vertices_per_cell)[i] = nv;
892
893 *cell_to_vertex = all_cell_to_vertex;
894
895 *global_cell_ids =
896 xmalloc((num_local_cells + num_halo_cells) * sizeof(**global_cell_ids));
897 for (size_t i = 0; i < num_local_cells; ++i)
898 (*global_cell_ids)[i] = local_start_cell + i;
899 memcpy(*global_cell_ids + num_local_cells, halo_cells,
900 num_halo_cells * sizeof(*halo_cells));
901 *global_vertex_ids = xrealloc(all_local_vertices, num_all_local_vertices *
902 sizeof(*all_local_vertices));
903
904 *cell_owner =
905 xmalloc((num_local_cells + num_halo_cells) * sizeof(**cell_owner));
906 for (size_t i = 0; i < num_local_cells; ++i)
907 (*cell_owner)[i] = -1;
908 for (size_t i = 0; i < num_halo_cells; ++i)
909 (*cell_owner)[num_local_cells + i] =
910 ((unsigned long)(halo_cells[i]) * (unsigned long)comm_size +
911 (unsigned long)comm_size - 1) / (unsigned long)ncells;
912 free(halo_cells);
913
914 for (size_t i = 0; i < num_all_local_vertices; ++i)
915 if (all_local_vertices_owner[i] == comm_rank)
916 all_local_vertices_owner[i] = -1;
917 *vertex_owner = all_local_vertices_owner;
918
919 *x_vertices = all_vertex_lon;
920 *y_vertices = all_vertex_lat;
921 if (with_cell_coords) {
922 *x_cells = cell_lon;
923 *y_cells = cell_lat;
924 }
925 if (with_cell_mask) *cell_msk = cell_mask;
926}
927
929 const char * filename, MPI_Fint comm, int * nbr_vertices, int * nbr_cells,
930 int ** num_vertices_per_cell, int ** cell_to_vertex, int ** global_cell_ids,
931 int ** cell_owner, int ** global_vertex_ids, int ** vertex_owner,
932 double ** x_vertices, double ** y_vertices,
933 double ** x_cells, double ** y_cells, int ** cell_msk) {
934
936 filename, MPI_Comm_f2c(comm), nbr_vertices, nbr_cells,
937 num_vertices_per_cell, cell_to_vertex, global_cell_ids,
938 cell_owner, global_vertex_ids, vertex_owner,
939 x_vertices, y_vertices, x_cells, y_cells, cell_msk);
940}
941
942static int * generate_simple_core_mask(size_t N) {
943 int * mask = xmalloc(N * sizeof(*mask));
944 for (size_t i = 0; i < N; ++i) mask[i] = 1;
945 return mask;
946}
947
948static size_t * generate_offsets(size_t N, int * counts) {
949
950 size_t * offsets = xmalloc(N * sizeof(*offsets));
951 for (size_t i = 0, accu = 0; i < N; ++i) {
952 offsets[i] = accu;
953 accu += (size_t)(counts[i]);
954 }
955 return offsets;
956}
957
959 const char * filename, MPI_Comm comm,
960 double ** x_vertices, double ** y_vertices,
961 yac_int ** cell_ids, yac_int ** vertex_ids, yac_int ** edge_ids,
962 size_t * num_cells, size_t * num_vertices, size_t * num_edges,
963 int ** num_vertices_per_cell, int ** num_cells_per_vertex,
964 size_t ** cell_to_vertex, size_t ** cell_to_edge, size_t ** vertex_to_cell,
965 size_t ** edge_to_vertex, enum yac_edge_type ** edge_type,
966 double ** x_cells, double ** y_cells, int ** cell_msk) {
967
969 ((x_cells == NULL) && (y_cells == NULL)) ||
970 ((x_cells != NULL) && (y_cells != NULL)),
971 "ERROR(yac_read_icon_grid_information_parallel_2): "
972 "arguments x_cells and y_cells have to be either both NULL or "
973 "both different from NULL");
974
975 const int with_cell_coords = x_cells != NULL;
976 const int with_cell_mask = cell_msk != NULL;
977
978 int comm_rank, comm_size;
979
980 MPI_Comm_rank(comm, &comm_rank);
981 MPI_Comm_size(comm, &comm_size);
982
983 int local_is_io, * io_ranks, num_io_ranks;
984 yac_get_io_ranks(comm, &local_is_io, &io_ranks, &num_io_ranks);
985
986 size_t ncells, nvertices, nedges, nv, ne;
987
988 size_t read_local_start_cell = 0;
989 size_t read_num_local_cells = 0;
990 size_t read_local_start_vertex = 0;
991 size_t read_num_local_vertices = 0;
992 size_t read_local_start_edge = 0;
993 size_t read_num_local_edges = 0;
994 double * read_cell_lon = NULL;
995 double * read_cell_lat = NULL;
996 int * read_cell_mask = NULL;
997 double * read_vertex_lon = NULL;
998 double * read_vertex_lat = NULL;
999 int * read_dist_vertex_of_cell = NULL;
1000 int * read_dist_edge_of_cell = NULL;
1001 int * read_dist_edge_vertices = NULL;
1002
1003 if (local_is_io) {
1004
1005 unsigned long io_proc_idx = ULONG_MAX;
1006 for (int i = 0; (i < num_io_ranks) && (io_proc_idx == ULONG_MAX); ++i)
1007 if (io_ranks[i] == comm_rank)
1008 io_proc_idx = (unsigned long)i;
1009
1010 // open file
1011 int ncid;
1012 yac_nc_open(filename, NC_NOWRITE, &ncid);
1013
1014 // get number of cells and vertices
1015 int dim_id;
1016 yac_nc_inq_dimid(ncid, "cell", &dim_id);
1017 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dim_id, &ncells));
1018 yac_nc_inq_dimid(ncid, "vertex", &dim_id);
1019 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dim_id, &nvertices));
1020 yac_nc_inq_dimid(ncid, "edge", &dim_id);
1021 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dim_id, &nedges));
1022 yac_nc_inq_dimid(ncid, "nv", &dim_id);
1023 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dim_id, &nv));
1024 yac_nc_inq_dimid(ncid, "ne", &dim_id);
1025 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dim_id, &ne));
1026
1027 // determine local range for cell and vertex data
1028 read_local_start_cell =
1029 ((unsigned long)ncells * io_proc_idx) / (unsigned long)num_io_ranks;
1030 read_num_local_cells =
1031 ((unsigned long)ncells * (io_proc_idx+1)) / (unsigned long)num_io_ranks -
1032 (unsigned long)read_local_start_cell;
1033 read_local_start_vertex =
1034 ((unsigned long)nvertices * io_proc_idx) / (unsigned long)num_io_ranks;
1035 read_num_local_vertices =
1036 ((unsigned long)nvertices * (io_proc_idx+1)) / (unsigned long)num_io_ranks -
1037 (unsigned long)read_local_start_vertex;
1038 read_local_start_edge =
1039 ((unsigned long)nedges * io_proc_idx) / (unsigned long)num_io_ranks;
1040 read_num_local_edges =
1041 ((unsigned long)nedges * (io_proc_idx+1)) / (unsigned long)num_io_ranks -
1042 (unsigned long)read_local_start_edge;
1043
1044 // read basic grid data (each process its individual part)
1045 read_vertex_lon = xmalloc(read_num_local_vertices * sizeof(*read_vertex_lon));
1046 read_vertex_lat = xmalloc(read_num_local_vertices * sizeof(*read_vertex_lat));
1047 read_cell_lon =
1048 with_cell_coords?
1049 xmalloc(read_num_local_cells * sizeof(*read_cell_lon)):NULL;
1050 read_cell_lat =
1051 with_cell_coords?
1052 xmalloc(read_num_local_cells * sizeof(*read_cell_lat)):NULL;
1053 read_cell_mask =
1054 with_cell_mask?
1055 xmalloc(read_num_local_cells * sizeof(*read_cell_mask)):NULL;
1056 read_dist_vertex_of_cell =
1057 xmalloc(read_num_local_cells * nv * sizeof(*read_dist_vertex_of_cell));
1058 read_dist_edge_of_cell =
1059 xmalloc(read_num_local_cells * nv * sizeof(*read_dist_edge_of_cell));
1060 read_dist_edge_vertices =
1061 xmalloc(read_num_local_edges * 2 * sizeof(*read_dist_edge_vertices));
1062 int varid;
1063 yac_nc_inq_varid(ncid, "vlon", &varid);
1064 YAC_HANDLE_ERROR(nc_get_vara_double(ncid, varid, &read_local_start_vertex,
1065 &read_num_local_vertices, read_vertex_lon));
1066 convert_to_rad(ncid, varid, read_vertex_lon, read_num_local_vertices);
1067 yac_nc_inq_varid(ncid, "vlat", &varid);
1068 YAC_HANDLE_ERROR(nc_get_vara_double(ncid, varid, &read_local_start_vertex,
1069 &read_num_local_vertices, read_vertex_lat));
1070 convert_to_rad(ncid, varid, read_vertex_lat, read_num_local_vertices);
1071 if (with_cell_coords) {
1072 yac_nc_inq_varid(ncid, "clon", &varid);
1074 nc_get_vara_double(
1075 ncid, varid, &read_local_start_cell,
1076 &read_num_local_cells, read_cell_lon));
1077 convert_to_rad(ncid, varid, read_cell_lon, read_num_local_cells);
1078 yac_nc_inq_varid(ncid, "clat", &varid);
1080 nc_get_vara_double(
1081 ncid, varid, &read_local_start_cell,
1082 &read_num_local_cells, read_cell_lat));
1083 convert_to_rad(ncid, varid, read_cell_lat, read_num_local_cells);
1084 }
1085 if (read_cell_mask) {
1086 yac_nc_inq_varid(ncid, "cell_sea_land_mask", &varid);
1088 nc_get_vara_int(
1089 ncid, varid, &read_local_start_cell,
1090 &read_num_local_cells, read_cell_mask));
1091 }
1092 {
1093 int * buffer =
1094 xmalloc(
1095 MAX(nv*read_num_local_cells, 2*read_num_local_edges) *
1096 sizeof(*buffer));
1097 {
1098 size_t tmp_start[2] = {0, read_local_start_cell};
1099 size_t tmp_count[2] = {nv, read_num_local_cells};
1100 yac_nc_inq_varid(ncid, "vertex_of_cell", &varid);
1101 YAC_HANDLE_ERROR(nc_get_vara_int(ncid, varid, tmp_start, tmp_count, buffer));
1102 for (size_t i = 0; i < read_num_local_cells; ++i)
1103 for (size_t j = 0; j < nv; ++j)
1104 read_dist_vertex_of_cell[i * nv + j] =
1105 buffer[i + j * read_num_local_cells];
1106 yac_nc_inq_varid(ncid, "edge_of_cell", &varid);
1107 YAC_HANDLE_ERROR(nc_get_vara_int(ncid, varid, tmp_start, tmp_count, buffer));
1108 for (size_t i = 0; i < read_num_local_cells; ++i)
1109 for (size_t j = 0; j < nv; ++j)
1110 read_dist_edge_of_cell[i * nv + j] =
1111 buffer[i + j * read_num_local_cells];
1112 }
1113 {
1114 size_t tmp_start[2] = {0, read_local_start_edge};
1115 size_t tmp_count[2] = {2, read_num_local_edges};
1116 yac_nc_inq_varid(ncid, "edge_vertices", &varid);
1117 YAC_HANDLE_ERROR(nc_get_vara_int(ncid, varid, tmp_start, tmp_count, buffer));
1118 for (size_t i = 0; i < read_num_local_edges; ++i)
1119 for (int j = 0; j < 2; ++j)
1120 read_dist_edge_vertices[i * 2 + j] =
1121 buffer[i + j * read_num_local_edges];
1122 }
1123 free(buffer);
1124
1125 // adjust for c indices
1126 for (size_t i = 0; i < read_num_local_cells * nv; ++i)
1127 if (read_dist_vertex_of_cell[i] > 0) read_dist_vertex_of_cell[i]--;
1128 for (size_t i = 0; i < read_num_local_cells * nv; ++i)
1129 if (read_dist_edge_of_cell[i] > 0) read_dist_edge_of_cell[i]--;
1130 for (size_t i = 0; i < read_num_local_edges * 2; ++i)
1131 if (read_dist_edge_vertices[i] > 0) read_dist_edge_vertices[i]--;
1132 }
1133
1134 YAC_HANDLE_ERROR(nc_close(ncid));
1135 } else {
1136 read_cell_lon = with_cell_coords?xmalloc(1 * sizeof(*read_cell_lon)):NULL;
1137 read_cell_lat = with_cell_coords?xmalloc(1 * sizeof(*read_cell_lat)):NULL;
1138 read_cell_mask = with_cell_mask?xmalloc(1 * sizeof(*read_cell_mask)):NULL;
1139 read_vertex_lon = xmalloc(1 * sizeof(*read_vertex_lon));
1140 read_vertex_lat = xmalloc(1 * sizeof(*read_vertex_lat));
1141 read_dist_vertex_of_cell = xmalloc(1 * sizeof(*read_dist_vertex_of_cell));
1142 read_dist_edge_of_cell = xmalloc(1 * sizeof(*read_dist_edge_of_cell));
1143 read_dist_edge_vertices = xmalloc(1 * sizeof(*read_dist_edge_vertices));
1144 }
1145
1146 free(io_ranks);
1147
1148 {
1149 int tmp;
1150 if (comm_rank == 0) tmp = (int)ncells;
1151 MPI_Bcast(&tmp, 1, MPI_INT, 0, comm);
1152 ncells = (size_t)tmp;
1153 if (comm_rank == 0) tmp = (int)nvertices;
1154 MPI_Bcast(&tmp, 1, MPI_INT, 0, comm);
1155 nvertices = (size_t)tmp;
1156 if (comm_rank == 0) tmp = (int)nedges;
1157 MPI_Bcast(&tmp, 1, MPI_INT, 0, comm);
1158 nedges = (size_t)tmp;
1159 if (comm_rank == 0) tmp = (int)nv;
1160 MPI_Bcast(&tmp, 1, MPI_INT, 0, comm);
1161 nv = (size_t)tmp;
1162 if (comm_rank == 0) tmp = (int)ne;
1163 MPI_Bcast(&tmp, 1, MPI_INT, 0, comm);
1164 ne = (size_t)tmp;
1165 }
1166
1167 // determine local range for cell and vertex data
1168 size_t local_start_cell =
1169 ((unsigned long)ncells * (unsigned long)comm_rank) /
1170 (unsigned long)comm_size;
1171 size_t num_local_cells =
1172 ((unsigned long)ncells * ((unsigned long)comm_rank+1)) /
1173 (unsigned long)comm_size - (unsigned long)local_start_cell;
1174 size_t local_start_vertex =
1175 ((unsigned long)nvertices * (unsigned long)comm_rank) /
1176 (unsigned long)comm_size;
1177 size_t num_local_vertices =
1178 ((unsigned long)nvertices * ((unsigned long)comm_rank+1)) /
1179 (unsigned long)comm_size - (unsigned long)local_start_vertex;
1180 size_t local_start_edge =
1181 ((unsigned long)nedges * (unsigned long)comm_rank) /
1182 (unsigned long)comm_size;
1183 size_t num_local_edges =
1184 ((unsigned long)nedges * ((unsigned long)comm_rank+1)) /
1185 (unsigned long)comm_size - (unsigned long)local_start_edge;
1186
1187 // redistribute basic cell data (from io decomposition)
1188 if (with_cell_coords) {
1189 *x_cells = xmalloc(num_local_cells * sizeof(**x_cells));
1190 *y_cells = xmalloc(num_local_cells * sizeof(**y_cells));
1191 }
1192 if (with_cell_mask) {
1193 *cell_msk = xmalloc(num_local_cells * sizeof(**cell_msk));
1194 }
1195 int * dist_vertex_of_cell =
1196 xmalloc(num_local_cells * nv * sizeof(*dist_vertex_of_cell));
1197 int * dist_edge_of_cell =
1198 xmalloc(num_local_cells * nv * sizeof(*dist_edge_of_cell));
1199 {
1200 int * send_count = xcalloc(comm_size, sizeof(*send_count));
1201 int * recv_count = xcalloc(comm_size, sizeof(*recv_count));
1202
1203 for (unsigned i = 0; i < read_num_local_cells; ++i)
1204 send_count[
1206 read_local_start_cell + i, ncells, comm_size)]++;
1207
1208 MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
1209
1210 int * send_displ = xmalloc(comm_size * sizeof(*send_displ));
1211 int * recv_displ = xmalloc(comm_size * sizeof(*recv_displ));
1212 int send_accum = 0, recv_accum = 0;
1213 for (int i = 0; i < comm_size; ++i) {
1214 send_displ[i] = send_accum;
1215 recv_displ[i] = recv_accum;
1216 send_accum += send_count[i];
1217 recv_accum += recv_count[i];
1218 }
1219
1220 if (with_cell_coords) {
1221 MPI_Alltoallv(read_cell_lon, send_count, send_displ, MPI_DOUBLE,
1222 *x_cells, recv_count, recv_displ, MPI_DOUBLE, comm);
1223 MPI_Alltoallv(read_cell_lat, send_count, send_displ, MPI_DOUBLE,
1224 *y_cells, recv_count, recv_displ, MPI_DOUBLE, comm);
1225 }
1226 if (with_cell_mask) {
1227 MPI_Alltoallv(read_cell_mask, send_count, send_displ, MPI_INT,
1228 *cell_msk, recv_count, recv_displ, MPI_INT, comm);
1229 }
1230
1231 for (int i = 0; i < comm_size; ++i) {
1232 send_count[i] *= nv;
1233 send_displ[i] *= nv;
1234 recv_count[i] *= nv;
1235 recv_displ[i] *= nv;
1236 }
1237
1238 MPI_Alltoallv(read_dist_vertex_of_cell, send_count, send_displ, MPI_INT,
1239 dist_vertex_of_cell, recv_count, recv_displ, MPI_INT, comm);
1240 MPI_Alltoallv(read_dist_edge_of_cell, send_count, send_displ, MPI_INT,
1241 dist_edge_of_cell, recv_count, recv_displ, MPI_INT, comm);
1242
1243 free(recv_displ);
1244 free(send_displ);
1245 free(recv_count);
1246 free(send_count);
1247 free(read_dist_vertex_of_cell);
1248 free(read_dist_edge_of_cell);
1249 free(read_cell_mask);
1250 free(read_cell_lat);
1251 free(read_cell_lon);
1252 }
1253
1254 // redistribute basic vertex data (from io decomposition)
1255 double * vertex_lon = xmalloc(num_local_vertices * sizeof(*vertex_lon));
1256 double * vertex_lat = xmalloc(num_local_vertices * sizeof(*vertex_lat));
1257 {
1258 int * send_count = xcalloc(comm_size, sizeof(*send_count));
1259 int * recv_count = xcalloc(comm_size, sizeof(*recv_count));
1260
1261 for (unsigned i = 0; i < read_num_local_vertices; ++i)
1262 send_count[
1264 read_local_start_vertex + i, nvertices, comm_size)]++;
1265
1266 MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
1267
1268 int * send_displ = xmalloc(comm_size * sizeof(*send_displ));
1269 int * recv_displ = xmalloc(comm_size * sizeof(*recv_displ));
1270 int send_accum = 0, recv_accum = 0;
1271 for (int i = 0; i < comm_size; ++i) {
1272 send_displ[i] = send_accum;
1273 recv_displ[i] = recv_accum;
1274 send_accum += send_count[i];
1275 recv_accum += recv_count[i];
1276 }
1277
1278 MPI_Alltoallv(read_vertex_lon, send_count, send_displ, MPI_DOUBLE,
1279 vertex_lon, recv_count, recv_displ, MPI_DOUBLE, comm);
1280 MPI_Alltoallv(read_vertex_lat, send_count, send_displ, MPI_DOUBLE,
1281 vertex_lat, recv_count, recv_displ, MPI_DOUBLE, comm);
1282
1283 free(recv_displ);
1284 free(send_displ);
1285 free(recv_count);
1286 free(send_count);
1287 free(read_vertex_lon);
1288 free(read_vertex_lat);
1289 }
1290
1291 // redistribute basic edge data (from io decomposition)
1292 int * dist_edge_vertices =
1293 xmalloc(num_local_edges * 2 * sizeof(*dist_edge_vertices));
1294 {
1295 int * send_count = xcalloc(comm_size, sizeof(*send_count));
1296 int * recv_count = xcalloc(comm_size, sizeof(*recv_count));
1297
1298 for (unsigned i = 0; i < read_num_local_edges; ++i)
1299 send_count[
1301 read_local_start_edge + i, nedges, comm_size)] += 2;
1302
1303 MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
1304
1305 int * send_displ = xmalloc(comm_size * sizeof(*send_displ));
1306 int * recv_displ = xmalloc(comm_size * sizeof(*recv_displ));
1307 int send_accum = 0, recv_accum = 0;
1308 for (int i = 0; i < comm_size; ++i) {
1309 send_displ[i] = send_accum;
1310 recv_displ[i] = recv_accum;
1311 send_accum += send_count[i];
1312 recv_accum += recv_count[i];
1313 }
1314
1315 MPI_Alltoallv(read_dist_edge_vertices, send_count, send_displ, MPI_INT,
1316 dist_edge_vertices, recv_count, recv_displ, MPI_INT, comm);
1317
1318 free(recv_displ);
1319 free(send_displ);
1320 free(recv_count);
1321 free(send_count);
1322 free(read_dist_edge_vertices);
1323 }
1324
1325 // determine required vertices for core cells
1326 // in additional compute num_cells_per_vertex, vertex_to_cell,
1327 // and cell_to_vertex
1328 size_t num_core_vertices;
1329 {
1330 size_t N = num_local_cells * nv;
1331 *vertex_ids = xmalloc(N * sizeof(**vertex_ids));
1332 *num_cells_per_vertex = xmalloc(N * sizeof(**num_cells_per_vertex));
1333 *vertex_to_cell = xmalloc(N * sizeof(**vertex_to_cell));
1334 *cell_to_vertex = xmalloc(N * sizeof(**cell_to_vertex));
1335 for (size_t i = 0; i < N; ++i)
1336 (*vertex_ids)[i] = (yac_int)dist_vertex_of_cell[i];
1337 size_t * permutation = *vertex_to_cell;
1338 for (size_t i = 0; i < N; ++i) permutation[i] = i;
1339 yac_quicksort_index_yac_int_size_t(*vertex_ids, N, permutation);
1340 // remove duplicated core vertices and count number of cells per vertex
1341 yac_int prev_vertex_id = XT_INT_MAX;
1342 num_core_vertices = 0;
1343 for (size_t i = 0; i < N; ++i) {
1344 yac_int curr_vertex_id = (*vertex_ids)[i];
1345 if (prev_vertex_id == curr_vertex_id) {
1346 (*num_cells_per_vertex)[num_core_vertices-1]++;
1347 } else {
1348 (*num_cells_per_vertex)[num_core_vertices] = 1;
1349 (*vertex_ids)[num_core_vertices] = (prev_vertex_id = curr_vertex_id);
1350 ++num_core_vertices;
1351 }
1352 (*cell_to_vertex)[permutation[i]] = num_core_vertices-1;
1353 permutation[i] /= nv;
1354 }
1355 *vertex_ids =
1356 xrealloc(*vertex_ids, num_core_vertices * sizeof(**vertex_ids));
1357 *num_cells_per_vertex =
1358 xrealloc(*num_cells_per_vertex,
1359 num_core_vertices * sizeof(**num_cells_per_vertex));
1360 free(dist_vertex_of_cell);
1361 }
1362
1363 // determine required edges for core cells
1364 // in additional compute edge_to_cell and cell_to_edge
1365 size_t num_core_edges;
1366 {
1367 size_t N = num_local_cells * nv;
1368 *edge_ids = xmalloc(N * sizeof(**edge_ids));
1369 size_t * permutation = xmalloc(N * sizeof(*permutation));
1370 *cell_to_edge = xmalloc(N * sizeof(**cell_to_edge));
1371 for (size_t i = 0; i < N; ++i)
1372 (*edge_ids)[i] = (yac_int)dist_edge_of_cell[i];
1373 for (size_t i = 0; i < N; ++i) permutation[i] = i;
1374 yac_quicksort_index_yac_int_size_t(*edge_ids, N, permutation);
1375 // remove duplicated core edges and count number of cells per edge
1376 yac_int prev_edge_id = XT_INT_MAX;
1377 num_core_edges = 0;
1378 for (size_t i = 0; i < N; ++i) {
1379 yac_int curr_edge_id = (*edge_ids)[i];
1380 if (prev_edge_id != curr_edge_id)
1381 (*edge_ids)[num_core_edges++] = (prev_edge_id = curr_edge_id);
1382 (*cell_to_edge)[permutation[i]] = num_core_edges-1;
1383 permutation[i] /= nv;
1384 }
1385 free(permutation);
1386 *edge_ids =
1387 xrealloc(*edge_ids, num_core_edges * sizeof(**edge_ids));
1388 free(dist_edge_of_cell);
1389 }
1390
1391 // generate vertex coordinate data
1392 {
1393 *x_vertices = xmalloc(num_core_vertices * sizeof(**x_vertices));
1394 *y_vertices = xmalloc(num_core_vertices * sizeof(**y_vertices));
1395 int * send_count = xcalloc(comm_size, sizeof(*send_count));
1396 int * recv_count = xcalloc(comm_size, sizeof(*recv_count));
1397
1398 for (size_t i = 0; i < num_core_vertices; ++i)
1399 send_count[((unsigned long)((*vertex_ids)[i]) *
1400 (unsigned long)comm_size + (unsigned long)comm_size - 1) /
1401 (unsigned long)nvertices]++;
1402
1403 MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
1404
1405 int * send_displ = xmalloc(comm_size * sizeof(*send_displ));
1406 int * recv_displ = xmalloc(comm_size * sizeof(*recv_displ));
1407 int send_accum = 0, recv_accum = 0;
1408 for (int i = 0; i < comm_size; ++i) {
1409 send_displ[i] = send_accum;
1410 recv_displ[i] = recv_accum;
1411 send_accum += send_count[i];
1412 recv_accum += recv_count[i];
1413 }
1414
1415 int num_all_local_vertices_remote = 0;
1416 for (int i = 0; i < comm_size; ++i)
1417 num_all_local_vertices_remote += recv_count[i];
1418
1419 yac_int * remote_vertex_buffer =
1420 xmalloc(num_all_local_vertices_remote * sizeof(*remote_vertex_buffer));
1421
1422 MPI_Alltoallv(
1423 *vertex_ids, send_count, send_displ, yac_int_dt,
1424 remote_vertex_buffer, recv_count, recv_displ, yac_int_dt, comm);
1425
1426 double * send_vertex_lon =
1427 xmalloc(num_all_local_vertices_remote * sizeof(*send_vertex_lon));
1428 double * send_vertex_lat =
1429 xmalloc(num_all_local_vertices_remote * sizeof(*send_vertex_lat));
1430
1431 for (int i = 0, l = 0; i < comm_size; ++i) {
1432 for (int j = 0; j < recv_count[i]; ++j, ++l) {
1433 size_t idx = (size_t)(remote_vertex_buffer[l]) - local_start_vertex;
1434 send_vertex_lon[l] = vertex_lon[idx];
1435 send_vertex_lat[l] = vertex_lat[idx];
1436 }
1437 }
1438
1439 free(remote_vertex_buffer);
1440 free(vertex_lon);
1441 free(vertex_lat);
1442
1443 MPI_Alltoallv(send_vertex_lon, recv_count, recv_displ, MPI_DOUBLE,
1444 *x_vertices, send_count, send_displ, MPI_DOUBLE, comm);
1445 MPI_Alltoallv(send_vertex_lat, recv_count, recv_displ, MPI_DOUBLE,
1446 *y_vertices, send_count, send_displ, MPI_DOUBLE, comm);
1447
1448 free(send_vertex_lat);
1449 free(send_vertex_lon);
1450 free(recv_displ);
1451 free(send_displ);
1452 free(recv_count);
1453 free(send_count);
1454 }
1455
1456 // generate edge vertex data
1457 *edge_to_vertex =
1458 xmalloc(2 * num_core_edges * sizeof(**edge_to_vertex));
1459 {
1460 int * local_edge_to_vertex =
1461 xmalloc(2 * num_core_edges * sizeof(*local_edge_to_vertex));
1462 int * send_count = xcalloc(comm_size, sizeof(*send_count));
1463 int * recv_count = xcalloc(comm_size, sizeof(*recv_count));
1464
1465 for (size_t i = 0; i < num_core_edges; ++i)
1466 send_count[((unsigned long)((*edge_ids)[i]) *
1467 (unsigned long)comm_size + (unsigned long)comm_size - 1) /
1468 (unsigned long)nedges]++;
1469
1470 MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
1471
1472 int * send_displ = xmalloc(comm_size * sizeof(*send_displ));
1473 int * recv_displ = xmalloc(comm_size * sizeof(*recv_displ));
1474 int send_accum = 0, recv_accum = 0;
1475 for (int i = 0; i < comm_size; ++i) {
1476 send_displ[i] = send_accum;
1477 recv_displ[i] = recv_accum;
1478 send_accum += send_count[i];
1479 recv_accum += recv_count[i];
1480 }
1481
1482 int num_all_local_edges_remote = 0;
1483 for (int i = 0; i < comm_size; ++i)
1484 num_all_local_edges_remote += recv_count[i];
1485
1486 yac_int * remote_edge_buffer =
1487 xmalloc(num_all_local_edges_remote * sizeof(*remote_edge_buffer));
1488
1489 MPI_Alltoallv(
1490 *edge_ids, send_count, send_displ, yac_int_dt,
1491 remote_edge_buffer, recv_count, recv_displ, yac_int_dt, comm);
1492
1493 int * send_edge_vertices =
1494 xmalloc(2 * num_all_local_edges_remote * sizeof(*send_edge_vertices));
1495
1496 for (int i = 0, l = 0; i < comm_size; ++i) {
1497 for (int j = 0; j < recv_count[i]; ++j, ++l) {
1498 size_t idx = (size_t)(remote_edge_buffer[l]) - local_start_edge;
1499 send_edge_vertices[2*l+0] = dist_edge_vertices[2*idx+0];
1500 send_edge_vertices[2*l+1] = dist_edge_vertices[2*idx+1];
1501
1502 }
1503 send_count[i] *= 2;
1504 recv_count[i] *= 2;
1505 send_displ[i] *= 2;
1506 recv_displ[i] *= 2;
1507 }
1508
1509 free(remote_edge_buffer);
1510 free(dist_edge_vertices);
1511
1512 MPI_Alltoallv(send_edge_vertices, recv_count, recv_displ, MPI_INT,
1513 local_edge_to_vertex, send_count, send_displ, MPI_INT, comm);
1514
1515 size_t * permutation = xmalloc(2 * num_core_edges * sizeof(*permutation));
1516 for (size_t i = 0; i < 2 * num_core_edges; ++i) permutation[i] = i;
1517
1519 local_edge_to_vertex, 2 * num_core_edges, permutation);
1520
1521 for (size_t i = 0, j = 0; i < 2 * num_core_edges; ++i) {
1522 yac_int curr_vertex = (yac_int)(local_edge_to_vertex[i]);
1523 while ((j < nvertices) && ((*vertex_ids)[j] < curr_vertex)) ++j;
1524 YAC_ASSERT(
1525 (j < nvertices) && ((*vertex_ids)[j] == curr_vertex),
1526 "ERROR(yac_read_icon_basic_grid_data_parallel): vertex id missmatch")
1527 (*edge_to_vertex)[permutation[i]] = j;
1528 }
1529
1530 free(permutation);
1531 free(send_edge_vertices);
1532 free(recv_displ);
1533 free(send_displ);
1534 free(recv_count);
1535 free(send_count);
1536 free(local_edge_to_vertex);
1537 }
1538
1539 // generate cell ids for local partition
1540 *cell_ids = xmalloc(num_local_cells * sizeof(**cell_ids));
1541 for (size_t i = 0; i < num_local_cells; ++i)
1542 (*cell_ids)[i] = (yac_int)(local_start_cell + i);
1543
1544 // generate num_vertices_per_cell
1545 *num_vertices_per_cell =
1546 xmalloc(num_local_cells * sizeof(**num_vertices_per_cell));
1547 for (size_t i = 0; i < num_local_cells; ++i) (*num_vertices_per_cell)[i] = nv;
1548
1549 // generate edge_type (for icon grids this is always YAC_GREAT_CIRCLE_EDGE)
1550 *edge_type = xmalloc(num_core_edges * sizeof(**edge_type));
1551 for (size_t i = 0; i < num_core_edges; ++i)
1552 (*edge_type)[i] = YAC_GREAT_CIRCLE_EDGE;
1553
1554 *num_cells = num_local_cells;
1555 *num_vertices = num_core_vertices;
1556 *num_edges = num_core_edges;
1557}
1558
1560 const char * filename, MPI_Comm comm) {
1561
1562 double * x_vertices;
1563 double * y_vertices;
1564 yac_int * cell_ids;
1566 yac_int * edge_ids;
1567 size_t num_cells;
1568 size_t num_vertices;
1569 size_t num_edges;
1572 size_t * cell_to_vertex;
1573 size_t * cell_to_edge;
1574 size_t * vertex_to_cell;
1575 size_t * edge_to_vertex;
1576 enum yac_edge_type * edge_type;
1577
1579 filename, comm, &x_vertices, &y_vertices, &cell_ids, &vertex_ids, &edge_ids,
1582 &edge_to_vertex, &edge_type, NULL, NULL, NULL);
1583
1586 for (size_t i = 0; i < num_vertices; ++i)
1587 LLtoXYZ(x_vertices[i], y_vertices[i], vertex_coordinates[i]);
1588 free(x_vertices);
1589 free(y_vertices);
1590
1592 grid_data.vertex_coordinates = vertex_coordinates;
1593 grid_data.cell_ids = cell_ids;
1594 grid_data.vertex_ids = vertex_ids;
1595 grid_data.edge_ids = edge_ids;
1596 grid_data.num_cells = num_cells;
1597 grid_data.num_vertices = num_vertices;
1598 grid_data.num_edges = num_edges;
1603 grid_data.num_cells_per_vertex = num_cells_per_vertex;
1605 grid_data.cell_to_vertex_offsets = generate_offsets(num_cells, num_vertices_per_cell);
1606 grid_data.cell_to_edge = cell_to_edge;
1607 grid_data.cell_to_edge_offsets = grid_data.cell_to_vertex_offsets;
1608 grid_data.vertex_to_cell = vertex_to_cell;
1610 grid_data.edge_to_vertex = (yac_size_t_2_pointer)&(edge_to_vertex[0]);
1611 grid_data.edge_type = edge_type;
1612 grid_data.num_total_cells = num_cells;
1613 grid_data.num_total_vertices = num_vertices;
1614 grid_data.num_total_edges = num_edges;
1615
1616 return grid_data;
1617}
1618
1620 char const * filename, char const * gridname, MPI_Comm comm,
1621 struct yac_basic_grid ** basic_grid, size_t * cell_coordinate_idx,
1622 int ** cell_mask) {
1623
1624 double * x_vertices;
1625 double * y_vertices;
1626 yac_int * cell_ids;
1628 yac_int * edge_ids;
1629 size_t num_cells;
1630 size_t num_vertices;
1631 size_t num_edges;
1634 size_t * cell_to_vertex;
1635 size_t * cell_to_edge;
1636 size_t * vertex_to_cell;
1637 size_t * edge_to_vertex;
1638 enum yac_edge_type * edge_type;
1639 double * x_cells;
1640 double * y_cells;
1641
1643 filename, comm, &x_vertices, &y_vertices, &cell_ids, &vertex_ids, &edge_ids,
1646 &edge_to_vertex, &edge_type, &x_cells, &y_cells, cell_mask);
1647
1650 for (size_t i = 0; i < num_vertices; ++i)
1651 LLtoXYZ(x_vertices[i], y_vertices[i], vertex_coordinates[i]);
1652 free(x_vertices);
1653 free(y_vertices);
1654
1655 struct yac_basic_grid_data basic_grid_data;
1656 basic_grid_data.vertex_coordinates = vertex_coordinates;
1657 basic_grid_data.cell_ids = cell_ids;
1658 basic_grid_data.vertex_ids = vertex_ids;
1659 basic_grid_data.edge_ids = edge_ids;
1660 basic_grid_data.num_cells = num_cells;
1661 basic_grid_data.num_vertices = num_vertices;
1662 basic_grid_data.num_edges = num_edges;
1668 basic_grid_data.cell_to_vertex = cell_to_vertex;
1670 basic_grid_data.cell_to_edge = cell_to_edge;
1671 basic_grid_data.cell_to_edge_offsets = basic_grid_data.cell_to_vertex_offsets;
1672 basic_grid_data.vertex_to_cell = vertex_to_cell;
1674 basic_grid_data.edge_to_vertex = (yac_size_t_2_pointer)&(edge_to_vertex[0]);
1675 basic_grid_data.edge_type = edge_type;
1676 basic_grid_data.num_total_cells = num_cells;
1677 basic_grid_data.num_total_vertices = num_vertices;
1678 basic_grid_data.num_total_edges = num_edges;
1679
1680 *basic_grid = yac_basic_grid_new(gridname, basic_grid_data);
1681
1682 yac_coordinate_pointer cell_coordinates =
1683 xmalloc(num_cells * sizeof(*cell_coordinates));
1684 for (size_t i = 0; i < num_cells; ++i)
1685 LLtoXYZ(x_cells[i], y_cells[i], cell_coordinates[i]);
1686 *cell_coordinate_idx =
1688 *basic_grid, YAC_LOC_CELL, cell_coordinates);
1689 free(x_cells);
1690 free(y_cells);
1691}
1692
1694 char const * filename, char const * gridname, MPI_Comm comm) {
1695
1696 return
1698 gridname,
1700}
1701
1703 char const * filename, char const * gridname, MPI_Fint comm) {
1704
1705 return
1707 filename, gridname, MPI_Comm_f2c(comm));
1708}
1709
1711 char const * filename) {
1712
1713 int nbr_vertices;
1714 int nbr_cells;
1715 int * num_vertices_per_cell = NULL;
1716 int * cell_to_vertex = NULL;
1717 int * cell_mask = NULL;
1718
1719 double * x_vertices = NULL;
1720 double * y_vertices = NULL;
1721 double * x_cells = NULL;
1722 double * y_cells = NULL;
1723
1724 yac_read_icon_grid_information(filename, &nbr_vertices, &nbr_cells,
1726 &x_vertices, &y_vertices,
1727 &x_cells, &y_cells,
1728 &cell_mask);
1729
1730 free(x_cells);
1731 free(y_cells);
1732 free(cell_mask);
1733
1734 struct yac_basic_grid_data grid =
1736 (size_t)nbr_vertices, (int)nbr_cells, num_vertices_per_cell,
1737 x_vertices, y_vertices, cell_to_vertex);
1739 free(x_vertices);
1740 free(y_vertices);
1741 free(cell_to_vertex);
1742
1743 return grid;
1744}
1745
1747 char const * filename, char const * gridname) {
1748
1749 return
1751 gridname, yac_read_icon_basic_grid_data(filename));
1752}
1753
1754/* ---------------------------------------------------------------- */
1755
1756static int * get_icon_cell_mask ( int ncid, size_t nbr_cells ) {
1757
1758 // get variable id
1759 int mask_id;
1760 yac_nc_inq_varid (ncid, "cell_sea_land_mask", &mask_id);
1761
1762 // check number of dimension (has to be 1)
1763 int ndims;
1764 YAC_HANDLE_ERROR(nc_inq_varndims(ncid, mask_id, &ndims));
1765 YAC_ASSERT(
1766 ndims == 1,
1767 "ERROR(get_icon_cell_mask): mask array has more than one dimension")
1768
1769 // get id of dimension
1770 int dimid;
1771 YAC_HANDLE_ERROR(nc_inq_vardimid(ncid, mask_id, &dimid));
1772
1773 // check size of mask (has to be equal to nbr_cells)
1774 size_t dimlen;
1775 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dimid, &dimlen));
1776 YAC_ASSERT(
1777 dimlen == nbr_cells,
1778 "ERROR(get_icon_cell_mask): invalid size of mask array")
1779
1780 // check mask type (has to be NC_INT)
1781 nc_type mask_type;
1782 YAC_HANDLE_ERROR(nc_inq_vartype(ncid, mask_id, &mask_type));
1783 YAC_ASSERT(
1784 mask_type == NC_INT, "ERROR(get_icon_cell_mask): invalid mask type")
1785
1786 // get and return mask
1787 int * cell_mask = xmalloc(nbr_cells * sizeof(*cell_mask));
1788 YAC_HANDLE_ERROR(nc_get_var_int (ncid, mask_id, cell_mask));
1789 return cell_mask;
1790}
1791
1792/* ---------------------------------------------------------------- */
1793
1795 int ncid, size_t nbr_cells, int ** vertex_of_cell, size_t * nv) {
1796
1797 // get variable id
1798 int conn_id;
1799 yac_nc_inq_varid(ncid, "vertex_of_cell", &conn_id);
1800
1801 // check number of dimension (has to be 1)
1802 int ndims;
1803 YAC_HANDLE_ERROR(nc_inq_varndims(ncid, conn_id, &ndims));
1804 YAC_ASSERT(
1805 ndims == 2,
1806 "ERROR(get_icon_connect): "
1807 "connectivity array has invalid number of dimensions")
1808
1809 // get ids of dimensions
1810 int dimids[2];
1811 YAC_HANDLE_ERROR(nc_inq_vardimid(ncid, conn_id, dimids));
1812
1813 // check size of dimensions
1814 size_t dimlen;
1815 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dimids[0], nv));
1816 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dimids[1], &dimlen));
1817 YAC_ASSERT(
1818 dimlen == nbr_cells,
1819 "ERROR(get_icon_connect): invalid size of second dimension of "
1820 "connectivity array (has to be nbr_cells)")
1821
1822 // get and return connectivity array
1823 *vertex_of_cell = xmalloc(*nv * nbr_cells * sizeof(**vertex_of_cell));
1824 YAC_HANDLE_ERROR(nc_get_var_int (ncid, conn_id, *vertex_of_cell));
1825}
1826
1828 int ** global_cell_id,
1829 int ** global_cell_id_rank,
1830 int ** num_vertices_per_cell,
1831 int ** global_corner_id,
1832 int ** global_corner_id_rank,
1833 int ** cell_to_vertex,
1834 double ** x_cells,
1835 double ** y_cells,
1836 double ** x_vertices,
1837 double ** y_vertices) {
1838
1839 free (*cell_mask);
1840 free (*global_cell_id);
1841 free (*global_cell_id_rank);
1842 free (*num_vertices_per_cell);
1843 free (*global_corner_id);
1844 free (*global_corner_id_rank);
1845 free (*cell_to_vertex);
1846 free (*x_cells);
1847 free (*y_cells);
1848 free (*x_vertices);
1849 free (*y_vertices);
1850}
1851
1852#else
1853
1855 char const * filename) {
1856
1857 UNUSED(filename);
1858 die(
1859 "ERROR(yac_basic_grid_data): "
1860 "YAC is built without the NetCDF support");
1861
1862 return
1864 (size_t[]){0,0}, (int[]){0,0}, NULL, NULL);
1865}
1866
1868 char const * filename, char const * gridname) {
1869
1870 UNUSED(filename);
1871 UNUSED(gridname);
1872 die(
1873 "ERROR(yac_read_icon_basic_grid): "
1874 "YAC is built without the NetCDF support");
1875
1876 return NULL;
1877}
1878
1879void yac_read_icon_grid_information(const char * filename, int * nbr_vertices,
1880 int * nbr_cells, int ** num_vertices_per_cell,
1881 int ** cell_to_vertex, double ** x_vertices,
1882 double ** y_vertices, double ** x_cells,
1883 double ** y_cells, int ** cell_mask) {
1884
1885 UNUSED(filename);
1886 UNUSED(nbr_vertices);
1887 UNUSED(nbr_cells);
1888 UNUSED(num_vertices_per_cell);
1890 UNUSED(x_vertices);
1891 UNUSED(y_vertices);
1892 UNUSED(x_cells);
1893 UNUSED(y_cells);
1895 die(
1896 "ERROR(yac_read_icon_grid_information): "
1897 "YAC is built without the NetCDF support");
1898}
1899
1901 const char * filename, int * num_vertices, int * num_cells,
1902 int ** num_vertices_per_cell, int ** cell_to_vertex,
1903 double ** x_vertices, double ** y_vertices,
1904 double ** x_cells, double ** y_cells, int ** global_cell_id,
1905 int ** cell_mask, int ** cell_core_mask,
1906 int ** global_corner_id, int ** corner_core_mask, int rank, int size) {
1907
1908 UNUSED(filename);
1909 UNUSED(num_vertices);
1911 UNUSED(num_vertices_per_cell);
1913 UNUSED(x_vertices);
1914 UNUSED(y_vertices);
1915 UNUSED(x_cells);
1916 UNUSED(y_cells);
1917 UNUSED(global_cell_id);
1920 UNUSED(global_corner_id);
1921 UNUSED(corner_core_mask);
1922 UNUSED(rank);
1923 UNUSED(size);
1924 die(
1925 "ERROR(yac_read_part_icon_grid_information): "
1926 "YAC is built without the NetCDF support");
1927}
1928
1930 const char * filename, MPI_Comm comm, int * nbr_vertices, int * nbr_cells,
1931 int ** num_vertices_per_cell, int ** cell_to_vertex, int ** global_cell_ids,
1932 int ** cell_owner, int ** global_vertex_ids, int ** vertex_owner,
1933 double ** x_vertices, double ** y_vertices,
1934 double ** x_cells, double ** y_cells, int ** cell_msk) {
1935
1936 UNUSED(filename);
1937 UNUSED(comm);
1938 UNUSED(nbr_vertices);
1939 UNUSED(nbr_cells);
1940 UNUSED(num_vertices_per_cell);
1942 UNUSED(global_cell_ids);
1943 UNUSED(cell_owner);
1944 UNUSED(global_vertex_ids);
1945 UNUSED(vertex_owner);
1946 UNUSED(x_vertices);
1947 UNUSED(y_vertices);
1948 UNUSED(x_cells);
1949 UNUSED(y_cells);
1950 UNUSED(cell_msk);
1951 die(
1952 "ERROR(yac_read_icon_grid_information_parallel): "
1953 "YAC is built without the NetCDF support");
1954}
1955
1957 const char * filename, MPI_Comm comm) {
1958
1959 UNUSED(filename);
1960 UNUSED(comm);
1961 die(
1962 "ERROR(yac_read_icon_basic_grid_data_parallel): "
1963 "YAC is built without the NetCDF support");
1964
1965 return
1967 (size_t[]){0,0}, (int[]){0,0}, NULL, NULL);
1968}
1969
1971 char const * filename, char const * gridname, MPI_Comm comm) {
1972
1973 UNUSED(filename);
1974 UNUSED(gridname);
1975 UNUSED(comm);
1976 die(
1977 "ERROR(yac_read_icon_basic_grid_parallel): "
1978 "YAC is built without the NetCDF support");
1979
1980 return NULL;
1981}
1982
1984 int ** global_cell_id,
1985 int ** cell_core_mask,
1986 int ** num_vertices_per_cell,
1987 int ** global_corner_id,
1988 int ** corner_core_mask,
1989 int ** cell_to_vertex,
1990 double ** x_cells,
1991 double ** y_cells,
1992 double ** x_vertices,
1993 double ** y_vertices) {
1994
1996 UNUSED(global_cell_id);
1998 UNUSED(num_vertices_per_cell);
1999 UNUSED(global_corner_id);
2000 UNUSED(corner_core_mask);
2002 UNUSED(x_vertices);
2003 UNUSED(y_vertices);
2004 UNUSED(x_cells);
2005 UNUSED(y_cells);
2006 die(
2007 "ERROR(yac_delete_icon_grid_data): "
2008 "YAC is built without the NetCDF support");
2009}
2010
2011#endif // YAC_NETCDF_ENABLED
#define YAC_ASSERT(exp, msg)
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_data yac_generate_basic_grid_data_reg_2d(size_t nbr_vertices[2], int cyclic[2], double *lon_vertices, double *lat_vertices)
Definition grid_reg2d.c:65
struct yac_basic_grid_data yac_generate_basic_grid_data_unstruct(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell, double *x_vertices, double *y_vertices, int *cell_to_vertex)
#define UNUSED(x)
Definition core.h:73
#define YAC_RAD
yac_edge_type
Definition grid_cell.h:12
@ YAC_GREAT_CIRCLE_EDGE
great circle
Definition grid_cell.h:13
void yac_get_io_ranks(MPI_Comm comm, int *local_is_io_, int **io_ranks_, int *num_io_ranks_)
Definition io_utils.c:309
void yac_nc_inq_varid(int ncid, char const *name, int *varidp)
Definition io_utils.c:411
void yac_nc_open(const char *path, int omode, int *ncidp)
Definition io_utils.c:350
void yac_nc_inq_dimid(int ncid, char const *name, int *dimidp)
Definition io_utils.c:385
@ YAC_LOC_CELL
Definition location.h:14
#define xrealloc(ptr, size)
Definition ppm_xfuncs.h:67
#define xcalloc(nmemb, size)
Definition ppm_xfuncs.h:64
#define xmalloc(size)
Definition ppm_xfuncs.h:66
int yac_check_coord_units(int ncid, int varid)
Definition read_grid.c:45
void yac_read_coords(int ncid, char const *lon_name, char const *lat_name, double **lon, double **lat, size_t *len)
Definition read_grid.c:94
struct yac_basic_grid * yac_read_icon_basic_grid_parallel_f2c(char const *filename, char const *gridname, MPI_Fint comm)
static void convert_to_rad(int ncid, int varid, double *array, size_t count)
void yac_delete_icon_grid_data(int **cell_mask, int **global_cell_id, int **global_cell_id_rank, int **num_vertices_per_cell, int **global_corner_id, int **global_corner_id_rank, int **cell_to_vertex, double **x_cells, double **y_cells, double **x_vertices, double **y_vertices)
struct yac_basic_grid * yac_read_icon_basic_grid(char const *filename, char const *gridname)
void yac_read_part_icon_grid_information(const char *filename, int *nbr_vertices, int *nbr_cells, int **num_vertices_per_cell, int **cell_to_vertex, double **x_vertices, double **y_vertices, double **x_cells, double **y_cells, int **global_cell_id, int **cell_mask, int **cell_core_mask, int **global_corner_id, int **corner_core_mask, int rank, int size)
struct yac_basic_grid_data yac_read_icon_basic_grid_data_parallel(const char *filename, MPI_Comm comm)
static int * get_icon_cell_mask(int ncid, size_t nbr_cells)
void yac_read_icon_grid_information_parallel(const char *filename, MPI_Comm comm, int *nbr_vertices, int *nbr_cells, int **num_vertices_per_cell, int **cell_to_vertex, int **global_cell_ids, int **cell_owner, int **global_vertex_ids, int **vertex_owner, double **x_vertices, double **y_vertices, double **x_cells, double **y_cells, int **cell_msk)
void yac_read_icon_basic_grid_parallel_2(char const *filename, char const *gridname, MPI_Comm comm, struct yac_basic_grid **basic_grid, size_t *cell_coordinate_idx, int **cell_mask)
void yac_read_icon_grid_information_parallel_f2c(const char *filename, MPI_Fint comm, int *nbr_vertices, int *nbr_cells, int **num_vertices_per_cell, int **cell_to_vertex, int **global_cell_ids, int **cell_owner, int **global_vertex_ids, int **vertex_owner, double **x_vertices, double **y_vertices, double **x_cells, double **y_cells, int **cell_msk)
static int partition_idx_from_element_idx(unsigned element_idx, unsigned num_elements, int num_partitions)
void yac_read_icon_grid_information_parallel_2(const char *filename, MPI_Comm comm, double **x_vertices, double **y_vertices, yac_int **cell_ids, yac_int **vertex_ids, yac_int **edge_ids, size_t *num_cells, size_t *num_vertices, size_t *num_edges, int **num_vertices_per_cell, int **num_cells_per_vertex, size_t **cell_to_vertex, size_t **cell_to_edge, size_t **vertex_to_cell, size_t **edge_to_vertex, enum yac_edge_type **edge_type, double **x_cells, double **y_cells, int **cell_msk)
static size_t * generate_offsets(size_t N, int *counts)
static int * generate_simple_core_mask(size_t N)
void yac_read_icon_grid_information(const char *filename, int *nbr_vertices, int *nbr_cells, int **num_vertices_per_cell, int **cell_to_vertex, double **x_vertices, double **y_vertices, double **x_cells, double **y_cells, int **cell_mask)
static void get_icon_connect(int ncid, size_t nbr_cells, int **vertex_of_cell, size_t *nv)
struct yac_basic_grid * yac_read_icon_basic_grid_parallel(char const *filename, char const *gridname, MPI_Comm comm)
struct yac_basic_grid_data yac_read_icon_basic_grid_data(char const *filename)
yac_coordinate_pointer vertex_coordinates
size_t * vertex_to_cell_offsets
yac_size_t_2_pointer edge_to_vertex
enum yac_edge_type * edge_type
size_t * cell_to_vertex_offsets
int * cell_to_vertex
int * cell_mask
double * cell_lat
double * cell_lon
#define N
size_t num_cells[2]
int vertex_of_cell[3][16]
static int mask[16]
int * cell_core_mask
double * buffer
#define YAC_HANDLE_ERROR(exp)
Definition toy_output.c:13
static void LLtoXYZ(double lon, double lat, double p_out[])
Definition toy_scrip.c:587
static void yac_remove_duplicates_int(int *array, size_t *n)
#define MAX(a, b)
void yac_quicksort_index_yac_int_size_t(yac_int *a, size_t n, size_t *idx)
void yac_quicksort_index_int_size_t(int *a, size_t n, size_t *idx)
void yac_quicksort_index(int *a, size_t n, int *idx)
#define die(msg)
Definition yac_assert.h:12
Xt_int yac_int
Definition yac_types.h:15
size_t(* yac_size_t_2_pointer)[2]
Definition yac_types.h:23
#define yac_int_dt
Definition yac_types.h:16
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19