YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
read_scrip_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 <stdint.h>
12#include <string.h>
13
14#include <mpi.h>
15
16#include "read_scrip_grid.h"
17#include "geometry.h"
18#include "utils_common.h"
19#include "io_utils.h"
20#include "yac_mpi_internal.h"
21
22#ifdef YAC_NETCDF_ENABLED
23
24#include <netcdf.h>
25
26struct point_with_index {
27
28 int32_t lon, lat;
29 int rank;
31 size_t i;
32 double dlon, dlat;
33};
34
40
47
49 double ** vertex_lon, double ** vertex_lat,
50 size_t * nbr_vertices, int * old_to_new_id);
52 double ** vertex_lon, double ** vertex_lat,
53 size_t * nbr_vertices, int * old_to_new_id,
54 yac_int ** vertex_ids, MPI_Comm comm);
55static void detect_duplicated_cells(
56 int * cell_to_vertex, int * num_vertices_per_cell, int * cell_mask,
57 size_t nbr_cells, size_t ** duplicated_cell_idx,
58 size_t ** orig_cell_idx, size_t * nbr_duplicated_cells);
60 int * cell_to_vertex, int * num_vertices_per_cell, int * cell_mask,
61 yac_int * cell_ids, size_t num_cells, yac_int * vertex_ids,
62 MPI_Comm comm, size_t ** duplicated_cell_idx,
63 yac_int ** orig_cell_ids, size_t * num_duplicated_cells_);
64
65static inline int compare_int(
66 const void * a,const void * b) {
67 return *(int const *)a - *(int const*)b;
68}
69
70static inline int compare_yac_int(
71 const void * a, const void * b) {
72 return (*(yac_int const *)a > *(yac_int const*)b) -
73 (*(yac_int const *)a < *(yac_int const*)b);
74}
75
77 const void * a_, const void * b_) {
78
79 struct cell_vertex_ids_with_index const * a =
80 (struct cell_vertex_ids_with_index const *)a_;
81 struct cell_vertex_ids_with_index const * b =
82 (struct cell_vertex_ids_with_index const *)b_;
83
84 size_t num_vertices_a = a->num_vertices;
85 size_t num_vertices_b = b->num_vertices;
86 int ret = (num_vertices_a > num_vertices_b) -
87 (num_vertices_a < num_vertices_b);
88
89 for (size_t i = 0; (i < num_vertices_a) && !ret; ++i)
90 ret = (a->vertex_ids[i] > b->vertex_ids[i]) -
91 (a->vertex_ids[i] < b->vertex_ids[i]);
92
93 return ret;
94}
95
97 const void * a_, const void * b_) {
98
100 if (ret) return ret;
101
102 struct cell_vertex_ids_with_index const * a =
103 (struct cell_vertex_ids_with_index const *)a_;
104 struct cell_vertex_ids_with_index const * b =
105 (struct cell_vertex_ids_with_index const *)b_;
106
107 return (a->cell_id > b->cell_id) - (a->cell_id < b->cell_id);
108}
109
111 const void * a_, const void * b_) {
112
113 struct point_with_index * a = (struct point_with_index *)a_;
114 struct point_with_index * b = (struct point_with_index *)b_;
115
116 int ret = (a->lon > b->lon) - (a->lon < b->lon);
117 if (ret) return ret;
118 ret = (a->lat > b->lat) - (a->lat < b->lat);
119 if (ret) return ret;
120 return (a->id > b->id) - (a->id < b->id);
121}
122
124 const void * a_, const void * b_) {
125
126 struct point_with_index * a = (struct point_with_index *)a_;
127 struct point_with_index * b = (struct point_with_index *)b_;
128
129 int ret = (a->lon > b->lon) - (a->lon < b->lon);
130 if (ret) return ret;
131 return (a->lat > b->lat) - (a->lat < b->lat);
132}
133
135 const void * a_, const void * b_) {
136
137 struct point_with_index * a = (struct point_with_index *)a_;
138 struct point_with_index * b = (struct point_with_index *)b_;
139
140 return (a->rank > b->rank) - (a->rank < b->rank);
141}
142
144 const void * a_, const void * b_) {
145
146 struct point_with_index * a = (struct point_with_index *)a_;
147 struct point_with_index * b = (struct point_with_index *)b_;
148
149 return (a->i > b->i) - (a->i < b->i);
150}
151
153 const char * filename, const char * grid_name,
154 size_t cell_mask_size, int * cell_mask,
155 size_t * num_vertices_, size_t * num_cells_, int ** num_vertices_per_cell_,
156 int ** cell_to_vertex_, double ** x_vertices, double ** y_vertices,
157 double ** x_cells, double ** y_cells,
158 size_t ** duplicated_cell_idx_, size_t ** orig_cell_idx_,
159 size_t * nbr_duplicated_cells_) {
160
161 size_t grid_name_len = strlen(grid_name) + 1;
162 char cla_var_name[4 + grid_name_len];
163 char clo_var_name[4 + grid_name_len];
164 char lat_var_name[4 + grid_name_len];
165 char lon_var_name[4 + grid_name_len];
166
167 snprintf(cla_var_name, 4 + grid_name_len, "%s.cla", grid_name);
168 snprintf(clo_var_name, 4 + grid_name_len, "%s.clo", grid_name);
169 snprintf(lat_var_name, 4 + grid_name_len, "%s.lat", grid_name);
170 snprintf(lon_var_name, 4 + grid_name_len, "%s.lon", grid_name);
171
172 int ncid;
173 yac_nc_open(filename, NC_NOWRITE, &ncid);
174
175 // get variable ids
176 int cla_var_id;
177 int clo_var_id;
178 int lat_var_id;
179 int lon_var_id;
180 yac_nc_inq_varid(ncid, cla_var_name, &cla_var_id);
181 yac_nc_inq_varid(ncid, clo_var_name, &clo_var_id);
182 yac_nc_inq_varid(ncid, lat_var_name, &lat_var_id);
183 yac_nc_inq_varid(ncid, lon_var_name, &lon_var_id);
184
185 // get dimension ids
186 int ndims;
187 int dim_ids[NC_MAX_VAR_DIMS];
189 nc_inq_var(ncid, cla_var_id, NULL, NULL, &ndims, dim_ids, NULL));
191 ndims == 3,
192 "ERROR(yac_read_scrip_basic_grid_information): "
193 "invalid number of dimensions for variable \"%s\" (has to be 3 not %d)",
194 cla_var_name, ndims);
195 int crn_dim_id = dim_ids[0];
196 int x_dim_id = dim_ids[2];
197 int y_dim_id = dim_ids[1];
198
199 // get dimension length
200 size_t crn_dim_len;
201 size_t x_dim_len;
202 size_t y_dim_len;
203 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, crn_dim_id, &crn_dim_len));
204 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, x_dim_id, &x_dim_len));
205 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, y_dim_id, &y_dim_len));
206
207 size_t num_cells = x_dim_len * y_dim_len;
208 size_t num_vertices = num_cells * crn_dim_len;
209
211 num_cells == cell_mask_size,
212 "ERROR(yac_read_scrip_basic_grid_information): "
213 "cell mask size is inconsistent with number of grid cells")
214
215 // allocate variables
216 double * cla = xmalloc(num_vertices * sizeof(*cla));
217 double * clo = xmalloc(num_vertices * sizeof(*clo));
218
219 //read variables
220 YAC_HANDLE_ERROR(nc_get_var_double(ncid, cla_var_id, cla));
221 YAC_HANDLE_ERROR(nc_get_var_double(ncid, clo_var_id, clo));
222
223 if ((x_cells != NULL) && (y_cells != NULL)) {
224 double * lat = xmalloc(num_cells * sizeof(*lat));
225 double * lon = xmalloc(num_cells * sizeof(*lon));
226 YAC_HANDLE_ERROR(nc_get_var_double(ncid, lat_var_id, lat));
227 YAC_HANDLE_ERROR(nc_get_var_double(ncid, lon_var_id, lon));
228 *x_cells = lon;
229 *y_cells = lat;
230 }
231
232 YAC_HANDLE_ERROR(nc_close(ncid));
233
234 size_t * reorder_idx = xmalloc(num_vertices * sizeof(*reorder_idx));
235 for (size_t y = 0, l = 0; y < y_dim_len; ++y)
236 for (size_t x = 0; x < x_dim_len; ++x)
237 for (size_t n = 0; n < crn_dim_len; ++n, ++l)
238 reorder_idx[x + y * x_dim_len + n * x_dim_len * y_dim_len] = l;
239
240 // remove duplicated vertices
241 int * cell_to_vertex = xmalloc(num_vertices * sizeof(*cell_to_vertex));
242 remove_duplicated_vertices(&clo, &cla, &num_vertices, cell_to_vertex);
243
244 *x_vertices = clo;
245 *y_vertices = cla;
246
247 // we have to reorder cell_to_vertex
249 reorder_idx, num_cells * crn_dim_len, cell_to_vertex);
250 free(reorder_idx);
251
252 // determine number of vertices per cell and compact cell_to_vertex
253 int * num_vertices_per_cell =
254 xmalloc(num_cells * sizeof(*num_vertices_per_cell));
255 size_t total_num_cell_vertices = 0;
256 int * to_vertices = cell_to_vertex;
257 int * from_vertices = cell_to_vertex;
258 for (size_t i = 0; i < num_cells; ++i) {
259 size_t curr_num_vertices = 0;
260 if (cell_mask[i]) {
261 int prev_vertex = from_vertices[crn_dim_len-1];
262 for (size_t j = 0; j < crn_dim_len; ++j, ++from_vertices) {
263 int curr_vertex = *from_vertices;
264 if (prev_vertex != curr_vertex) {
265 prev_vertex = curr_vertex;
266 if (to_vertices != from_vertices) *to_vertices = curr_vertex;
267 ++curr_num_vertices;
268 ++to_vertices;
269 }
270 }
271 } else {
272 from_vertices += crn_dim_len;
273 }
274 num_vertices_per_cell[i] = (int)curr_num_vertices;
275 total_num_cell_vertices += curr_num_vertices;
276 }
277
278 if (total_num_cell_vertices != num_cells * crn_dim_len)
279 cell_to_vertex =
280 xrealloc(
281 cell_to_vertex, total_num_cell_vertices * sizeof(*cell_to_vertex));
282
283 // make a copy of cell_to_vertex and sort the vertex indices of each cell,
284 // this is required by the detection of duplicated cells
285 int * cell_to_vertex_copy =
286 xmalloc(total_num_cell_vertices * sizeof(*cell_to_vertex_copy));
287 memcpy(
288 cell_to_vertex_copy, cell_to_vertex,
289 total_num_cell_vertices * sizeof(*cell_to_vertex_copy));
290 for (size_t i = 0, offset = 0; i < num_cells; ++i) {
291 qsort(
292 cell_to_vertex_copy + offset, (size_t)(num_vertices_per_cell[i]),
293 sizeof(*cell_to_vertex_copy), compare_int);
294 offset += num_vertices_per_cell[i];
295 }
296
297 // detect duplicated cells
298 size_t * duplicated_cell_idx;
299 size_t * orig_cell_idx;
300 size_t nbr_duplicated_cells;
302 cell_to_vertex_copy, num_vertices_per_cell, cell_mask, num_cells,
303 &duplicated_cell_idx, &orig_cell_idx, &nbr_duplicated_cells);
304 free(cell_to_vertex_copy);
305
306 // mask out duplicated cells
307 for (size_t i = 0; i < nbr_duplicated_cells; ++i)
308 cell_mask[duplicated_cell_idx[i]] = 0;
309
310 if ((duplicated_cell_idx_ != NULL) && (orig_cell_idx_ != NULL) &&
311 (nbr_duplicated_cells_ != NULL)) {
312 *duplicated_cell_idx_ = duplicated_cell_idx;
313 *orig_cell_idx_ = orig_cell_idx;
314 *nbr_duplicated_cells_ = nbr_duplicated_cells;
315 } else {
316 free(duplicated_cell_idx);
317 free(orig_cell_idx);
318 }
319
320 *num_vertices_ = num_vertices;
321 *num_cells_ = num_cells;
322 *num_vertices_per_cell_ = num_vertices_per_cell;
323 *cell_to_vertex_ = cell_to_vertex;
324}
325
327 const char * filename, const char * grid_name,
328 size_t cell_mask_size, int * cell_mask,
329 size_t * num_vertices_, size_t * num_cells_, int ** num_vertices_per_cell_,
330 int ** cell_to_vertex_,
331 double ** x_vertices_, double ** y_vertices_, yac_int ** vertex_ids_,
332 double ** x_cells_, double ** y_cells_, yac_int ** cell_ids_,
333 size_t ** duplicated_cell_idx_, yac_int ** orig_cell_ids_,
334 size_t * nbr_duplicated_cells_,
335 int io_rank_idx, int num_io_ranks, MPI_Comm comm) {
336
337 size_t num_cells = 0;
338 size_t num_vertices = 0;
339 size_t * reorder_idx = NULL;
340 size_t max_num_vertices_per_cell = 0;
341
342 double * x_vertices = NULL;
343 double * y_vertices = NULL;
344 yac_int * vertex_ids = NULL;
345 double * x_cells = NULL;
346 double * y_cells = NULL;
347 yac_int * cell_ids = NULL;
348
349 if ((io_rank_idx >= 0) && (io_rank_idx < num_io_ranks)) {
350
351 size_t grid_name_len = strlen(grid_name) + 1;
352 char cla_var_name[4 + grid_name_len];
353 char clo_var_name[4 + grid_name_len];
354 char lat_var_name[4 + grid_name_len];
355 char lon_var_name[4 + grid_name_len];
356
357 snprintf(cla_var_name, 4 + grid_name_len, "%s.cla", grid_name);
358 snprintf(clo_var_name, 4 + grid_name_len, "%s.clo", grid_name);
359 snprintf(lat_var_name, 4 + grid_name_len, "%s.lat", grid_name);
360 snprintf(lon_var_name, 4 + grid_name_len, "%s.lon", grid_name);
361
362 int ncid;
363 yac_nc_open(filename, NC_NOWRITE, &ncid);
364
365 // get variable ids
366 int cla_var_id, clo_var_id, lat_var_id, lon_var_id;
367 yac_nc_inq_varid(ncid, cla_var_name, &cla_var_id);
368 yac_nc_inq_varid(ncid, clo_var_name, &clo_var_id);
369 yac_nc_inq_varid(ncid, lat_var_name, &lat_var_id);
370 yac_nc_inq_varid(ncid, lon_var_name, &lon_var_id);
371
372 // get dimension ids
373 int ndims;
374 int dim_ids[NC_MAX_VAR_DIMS];
376 nc_inq_var(ncid, cla_var_id, NULL, NULL, &ndims, dim_ids, NULL));
378 ndims == 3,
379 "ERROR(yac_read_part_scrip_basic_grid_information): "
380 "invalid number of dimensions for variable \"%s\" (has to be 3 not %d)",
381 cla_var_name, ndims);
382 int crn_dim_id = dim_ids[0];
383 int x_dim_id = dim_ids[2];
384 int y_dim_id = dim_ids[1];
385
386 // get dimension length
387 size_t crn_dim_len, x_dim_len, y_dim_len;
388 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, crn_dim_id, &crn_dim_len));
389 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, x_dim_id, &x_dim_len));
390 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, y_dim_id, &y_dim_len));
391
392 // decompose x dimension among io ranks
393 size_t x_start_local =
394 (size_t)
395 (((unsigned long)x_dim_len * (unsigned long)io_rank_idx) /
396 (unsigned long)num_io_ranks);
397 size_t x_count_local =
398 (size_t)
399 (((unsigned long)x_dim_len * (io_rank_idx + 1)) /
400 (unsigned long)num_io_ranks) - x_start_local;
401
402 num_cells = x_count_local * y_dim_len;
403 num_vertices = num_cells * crn_dim_len;
404 max_num_vertices_per_cell = crn_dim_len;
405
407 num_cells == cell_mask_size,
408 "ERROR(yac_read_part_scrip_basic_grid_information): "
409 "cell mask size is inconsistent with number of grid cells")
410
411 // allocate variables
412 x_vertices = xmalloc(num_vertices * sizeof(*x_vertices));
413 y_vertices = xmalloc(num_vertices * sizeof(*y_vertices));
414
415 //read variables
416 size_t start[3] = {0, 0, x_start_local};
417 size_t count[3] = {crn_dim_len, y_dim_len, x_count_local};
419 nc_get_vara_double(ncid, clo_var_id, start, count, x_vertices));
421 nc_get_vara_double(ncid, cla_var_id, start, count, y_vertices));
422
423 if ((x_cells_ != NULL) && (y_cells_ != NULL)) {
424 x_cells = xmalloc(num_cells * sizeof(*x_cells));
425 y_cells = xmalloc(num_cells * sizeof(*y_cells));
426 size_t start[2] = {0, x_start_local};
427 size_t count[2] = {y_dim_len, x_count_local};
429 nc_get_vara_double(ncid, lon_var_id, start, count, x_cells));
431 nc_get_vara_double(ncid, lat_var_id, start, count, y_cells));
432 }
433
434 YAC_HANDLE_ERROR(nc_close(ncid));
435
436 reorder_idx = xmalloc(num_vertices * sizeof(*reorder_idx));
437 for (size_t y = 0, l = 0; y < y_dim_len; ++y)
438 for (size_t x = 0; x < x_count_local; ++x)
439 for (size_t n = 0; n < crn_dim_len; ++n, ++l)
440 reorder_idx[x + y * x_count_local + n * x_count_local * y_dim_len] = l;
441
442 // generate global cell ids
443 cell_ids = xmalloc(num_cells * sizeof(*cell_ids));
444 for (size_t y = 0, k = 0; y < y_dim_len; ++y)
445 for (size_t x = 0; x < x_count_local; ++x, ++k)
446 cell_ids[k] = y * x_dim_len + x + x_start_local;
447
448
450 (y_dim_len * x_dim_len * crn_dim_len) <= XT_INT_MAX,
451 "ERROR(yac_read_part_scrip_basic_grid_information): "
452 "number of vertices exceed maximum global id (%zu)",
453 (size_t)XT_INT_MAX);
454
455 // generate global vertex ids
456 vertex_ids = xmalloc(num_cells * crn_dim_len * sizeof(*vertex_ids));
457 for (size_t n = 0, l = 0; n < crn_dim_len; ++n)
458 for (size_t i = 0; i < num_cells; ++i, ++l)
459 vertex_ids[l] = cell_ids[i] * crn_dim_len + n;
460 }
461
462 // remove duplicated vertices
463 int * cell_to_vertex = xmalloc(num_vertices * sizeof(*cell_to_vertex));
465 &x_vertices, &y_vertices, &num_vertices,
466 cell_to_vertex, &vertex_ids, comm);
467
468 // we have to reorder cell_to_vertex
470 reorder_idx, num_cells * max_num_vertices_per_cell, cell_to_vertex);
471 free(reorder_idx);
472
473 // determine number of vertices per cell and compact cell_to_vertex
474 int * num_vertices_per_cell =
475 xmalloc(num_cells * sizeof(*num_vertices_per_cell));
476 size_t total_num_cell_vertices = 0;
477 int * to_vertices = cell_to_vertex;
478 int * from_vertices = cell_to_vertex;
479 for (size_t i = 0; i < num_cells; ++i) {
480 size_t curr_num_vertices = 0;
481 if (cell_mask[i]) {
482 int prev_vertex = from_vertices[max_num_vertices_per_cell-1];
483 for (size_t j = 0; j < max_num_vertices_per_cell; ++j, ++from_vertices) {
484 int curr_vertex = *from_vertices;
485 if (prev_vertex != curr_vertex) {
486 prev_vertex = curr_vertex;
487 if (to_vertices != from_vertices) *to_vertices = curr_vertex;
488 ++curr_num_vertices;
489 ++to_vertices;
490 }
491 }
492 } else {
493 from_vertices += max_num_vertices_per_cell;
494 }
495 num_vertices_per_cell[i] = (int)curr_num_vertices;
496 total_num_cell_vertices += curr_num_vertices;
497 }
498
499 if (total_num_cell_vertices != num_cells * max_num_vertices_per_cell)
500 cell_to_vertex =
501 xrealloc(
502 cell_to_vertex, total_num_cell_vertices * sizeof(*cell_to_vertex));
503
504 // detect duplicated cells
505 size_t * duplicated_cell_idx;
506 yac_int * orig_cell_ids;
507 size_t nbr_duplicated_cells;
508
510 cell_to_vertex, num_vertices_per_cell, cell_mask, cell_ids, num_cells,
511 vertex_ids, comm, &duplicated_cell_idx, &orig_cell_ids,
512 &nbr_duplicated_cells);
513
514 // mask out duplicated cells
515 for (size_t i = 0; i < nbr_duplicated_cells; ++i)
516 cell_mask[duplicated_cell_idx[i]] = 0;
517
518 if ((duplicated_cell_idx_ != NULL) && (orig_cell_ids_ != NULL) &&
519 (nbr_duplicated_cells_ != NULL)) {
520 *duplicated_cell_idx_ = duplicated_cell_idx;
521 *orig_cell_ids_ = orig_cell_ids;
522 *nbr_duplicated_cells_ = nbr_duplicated_cells;
523 } else {
524 free(duplicated_cell_idx);
525 free(orig_cell_ids);
526 }
527
528 *num_vertices_ = num_vertices;
529 *num_cells_ = num_cells;
530 *num_vertices_per_cell_ = num_vertices_per_cell;
531 *cell_to_vertex_ = cell_to_vertex;
532 *x_vertices_ = x_vertices;
533 *y_vertices_ = y_vertices;
534 *vertex_ids_ = vertex_ids;
535 if ((x_cells_ != NULL) && (y_cells_ != NULL)) {
536 *x_cells_ = x_cells;
537 *y_cells_ = y_cells;
538 }
539 *cell_ids_ = cell_ids;
540}
541
543 const char * filename, const char * grid_name, size_t * num_cells_,
544 int ** cell_mask) {
545
546 size_t grid_name_len = strlen(grid_name) + 1;
547 char msk_var_name[4 + grid_name_len];
548
549 snprintf(msk_var_name, 4 + grid_name_len, "%s.msk", grid_name);
550
551 int ncid;
552 yac_nc_open(filename, NC_NOWRITE, &ncid);
553
554 // get variable id
555 int msk_var_id;
556 yac_nc_inq_varid(ncid, msk_var_name, &msk_var_id);
557
558 // get dimension ids
559 int ndims;
560 int dim_ids[NC_MAX_VAR_DIMS];
562 nc_inq_var(ncid, msk_var_id, NULL, NULL, &ndims, dim_ids, NULL));
564 ndims == 2,
565 "ERROR(yac_read_scrip_mask_information): "
566 "invalid number of dimensions for variable \"%s\" (has to be 2 not %d)",
567 msk_var_name, ndims);
568 int x_dim_id = dim_ids[1];
569 int y_dim_id = dim_ids[0];
570
571 // get dimension length
572 size_t x_dim_len;
573 size_t y_dim_len;
574 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, x_dim_id, &x_dim_len));
575 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, y_dim_id, &y_dim_len));
576
577 size_t num_cells = x_dim_len * y_dim_len;
578
579 // allocate variable
580 *cell_mask = xmalloc(num_cells * sizeof(**cell_mask));
581
582 //read variable
583 YAC_HANDLE_ERROR(nc_get_var_int(ncid, msk_var_id, *cell_mask));
584
585 YAC_HANDLE_ERROR(nc_close(ncid));
586
587 *num_cells_ = num_cells;
588}
589
591 const char * filename, const char * grid_name, size_t * num_cells_,
592 int ** cell_mask, int io_rank_idx, int num_io_ranks) {
593
594 if ((io_rank_idx >= 0) && (io_rank_idx < num_io_ranks)) {
595
596 size_t grid_name_len = strlen(grid_name) + 1;
597 char msk_var_name[4 + grid_name_len];
598
599 snprintf(msk_var_name, 4 + grid_name_len, "%s.msk", grid_name);
600
601 int ncid;
602 yac_nc_open(filename, NC_NOWRITE, &ncid);
603
604 // get variable id
605 int msk_var_id;
606 yac_nc_inq_varid(ncid, msk_var_name, &msk_var_id);
607
608 // get dimension ids
609 int ndims;
610 int dim_ids[NC_MAX_VAR_DIMS];
612 nc_inq_var(ncid, msk_var_id, NULL, NULL, &ndims, dim_ids, NULL));
614 ndims == 2,
615 "ERROR(yac_read_part_scrip_mask_information): "
616 "invalid number of dimensions for variable \"%s\" (has to be 2 not %d)",
617 msk_var_name, ndims);
618 int x_dim_id = dim_ids[1];
619 int y_dim_id = dim_ids[0];
620
621 // get dimension length
622 size_t x_dim_len;
623 size_t y_dim_len;
624 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, x_dim_id, &x_dim_len));
625 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, y_dim_id, &y_dim_len));
626
627 // decompose x dimension among io ranks
628 size_t x_start_local =
629 (size_t)
630 (((unsigned long)x_dim_len * (unsigned long)io_rank_idx) /
631 (unsigned long)num_io_ranks);
632 size_t x_count_local =
633 (size_t)
634 (((unsigned long)x_dim_len * (io_rank_idx+1)) /
635 (unsigned long)num_io_ranks) - x_start_local;
636
637 size_t num_cells = x_count_local * y_dim_len;
638
639 // allocate variable
640 *cell_mask = xmalloc(num_cells * sizeof(**cell_mask));
641
642 //read variable
643 size_t start[2] = {0, x_start_local};
644 size_t count[2] = {y_dim_len, x_count_local};
646 nc_get_vara_int(ncid, msk_var_id, start, count, *cell_mask));
647
648 YAC_HANDLE_ERROR(nc_close(ncid));
649
650 *num_cells_ = num_cells;
651
652 } else {
653 *cell_mask = NULL;
654 *num_cells_ = 0;
655 }
656}
657
659 char const * grid_filename, char const * mask_filename,
660 char const * grid_name, int valid_mask_value,
661 size_t * num_vertices, size_t * num_cells, int ** num_vertices_per_cell,
662 double ** x_vertices, double ** y_vertices,
663 double ** x_cells, double ** y_cells,
664 int ** cell_to_vertex, int ** cell_core_mask, size_t ** duplicated_cell_idx,
665 size_t ** orig_cell_idx, size_t * nbr_duplicated_cells) {
666
667 size_t cell_mask_size;
668 int * cell_mask;
670 mask_filename, grid_name, &cell_mask_size, &cell_mask);
671
672 for (size_t i = 0; i < cell_mask_size; ++i)
673 cell_mask[i] = cell_mask[i] == valid_mask_value;
674
676 grid_filename, grid_name, cell_mask_size, cell_mask,
677 num_vertices, num_cells, num_vertices_per_cell, cell_to_vertex,
678 x_vertices, y_vertices, x_cells, y_cells, duplicated_cell_idx,
679 orig_cell_idx, nbr_duplicated_cells);
680
681 for (size_t i = 0; i < *num_vertices; ++i) {
682 (*x_vertices)[i] *= YAC_RAD;
683 (*y_vertices)[i] *= YAC_RAD;
684 }
685 if ((x_cells != NULL) && (y_cells != NULL)) {
686 for (size_t i = 0; i < *num_cells; ++i) {
687 (*x_cells)[i] *= YAC_RAD;
688 (*y_cells)[i] *= YAC_RAD;
689 }
690 }
691
693 *num_vertices <= INT_MAX,
694 "ERROR(yac_read_scrip_grid_information): "
695 "too man vertices in grid")
697 *num_cells <= INT_MAX,
698 "ERROR(yac_read_scrip_grid_information): "
699 "too man cells in grid")
700
701 *cell_core_mask = cell_mask;
702}
703
705 char const * grid_filename, char const * mask_filename,
706 MPI_Comm comm, char const * grid_name, int valid_mask_value,
707 size_t * num_vertices, size_t * num_cells, int ** num_vertices_per_cell,
708 double ** x_vertices, double ** y_vertices, yac_int ** vertex_ids,
709 double ** x_cells, double ** y_cells, yac_int ** cell_ids,
710 int ** cell_to_vertex, int ** cell_core_mask, size_t ** duplicated_cell_idx,
711 yac_int ** orig_cell_ids, size_t * nbr_duplicated_cells) {
712
713 // if no communicator was provided
714 if (comm == MPI_COMM_NULL) {
715
716 size_t * orig_cell_idx = NULL;
717
719 grid_filename, mask_filename, grid_name, valid_mask_value,
720 num_vertices, num_cells, num_vertices_per_cell,
721 x_vertices, y_vertices, x_cells, y_cells, cell_to_vertex,
722 cell_core_mask, duplicated_cell_idx,
723 (orig_cell_ids != NULL)?&orig_cell_idx:NULL,
724 nbr_duplicated_cells);
725
726 *vertex_ids = xmalloc(*num_vertices * sizeof(**vertex_ids));
727 for (size_t i = 0; i < *num_vertices; ++i) (*vertex_ids)[i] = i;
728
729 *cell_ids = xmalloc(*num_cells * sizeof(**cell_ids));
730 for (size_t i = 0; i < *num_cells; ++i) (*cell_ids)[i] = i;
731
732 if (orig_cell_ids != NULL) {
733 *orig_cell_ids =
734 xmalloc(*nbr_duplicated_cells * sizeof(**orig_cell_ids));
735 for (size_t i = 0; i < *nbr_duplicated_cells; ++i)
736 (*orig_cell_ids)[i] = (*cell_ids)[orig_cell_idx[i]];
737 free(orig_cell_idx);
738 }
739
740 } else {
741
742 // get io configuration
743 int local_is_io, * io_ranks, num_io_ranks;
744 yac_get_io_ranks(comm, &local_is_io, &io_ranks, &num_io_ranks);
745
746 int comm_rank;
747 yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);
748
749 int io_rank_idx = 0;
750 while ((io_rank_idx < num_io_ranks) &&
751 (comm_rank != io_ranks[io_rank_idx]))
752 ++io_rank_idx;
754 !local_is_io || (io_rank_idx < num_io_ranks),
755 "ERROR(yac_read_part_scrip_grid_information): "
756 "unable to determine io_rank_idx");
757
758 int * cell_mask;
759 size_t cell_mask_size;
761 mask_filename, grid_name, &cell_mask_size, &cell_mask,
762 io_rank_idx, num_io_ranks);
763
764 for (size_t i = 0; i < cell_mask_size; ++i)
765 cell_mask[i] = cell_mask[i] == valid_mask_value;
766
768 grid_filename, grid_name, cell_mask_size, cell_mask,
769 num_vertices, num_cells, num_vertices_per_cell, cell_to_vertex,
770 x_vertices, y_vertices, vertex_ids,
771 x_cells, y_cells, cell_ids, duplicated_cell_idx,
772 orig_cell_ids, nbr_duplicated_cells,
773 io_rank_idx, num_io_ranks, comm);
774
775 for (size_t i = 0; i < *num_vertices; ++i) {
776 (*x_vertices)[i] *= YAC_RAD;
777 (*y_vertices)[i] *= YAC_RAD;
778 }
779 if ((x_cells != NULL) && (y_cells != NULL)) {
780 for (size_t i = 0; i < *num_cells; ++i) {
781 (*x_cells)[i] *= YAC_RAD;
782 (*y_cells)[i] *= YAC_RAD;
783 }
784 }
785
787 *num_vertices <= INT_MAX,
788 "ERROR(yac_read_part_scrip_grid_information): "
789 "too man vertices in grid")
791 *num_cells <= INT_MAX,
792 "ERROR(yac_read_part_scrip_grid_information): "
793 "too man cells in grid")
795 cell_mask_size == *num_cells,
796 "ERROR(yac_read_part_scrip_grid_information): "
797 "mask and grid size do not match "
798 "(mask_file: \"%s\" grid_file: \"%s\" grid_name: \"%s\"",
799 mask_filename, grid_filename, grid_name)
800
801 *cell_core_mask = cell_mask;
802
803 free(io_ranks);
804 }
805}
806
808 char const * grid_filename, char const * mask_filename,
809 MPI_Comm comm, char const * grid_name, int valid_mask_value, int use_ll,
810 double ** x_cells, double ** y_cells, size_t ** duplicated_cell_idx,
811 yac_int ** orig_cell_global_ids, size_t * nbr_duplicated_cells) {
812
813 size_t num_vertices;
814 size_t num_cells;
815
817 double * x_vertices;
818 double * y_vertices;
821 int * cell_to_vertex;
822 int * cell_core_mask;
823
825 grid_filename, mask_filename, comm, grid_name, valid_mask_value,
827 &x_vertices, &y_vertices, &vertex_ids, x_cells, y_cells, &cell_ids,
828 &cell_to_vertex, &cell_core_mask, duplicated_cell_idx,
829 orig_cell_global_ids, nbr_duplicated_cells);
830
831 // check whether any cell has more then four vertices
832 // --> write out a warning and switch to gc edges
833 if (use_ll) {
834 for (size_t i = 0; (i < num_cells) && use_ll; ++i)
835 use_ll = !cell_core_mask[i] || (num_vertices_per_cell[i] <= 4);
836 int comm_rank = 0;
837 if (comm != MPI_COMM_NULL) {
839 MPI_Allreduce(
840 MPI_IN_PLACE, &use_ll, 1, MPI_INT, MPI_MIN, comm), comm);
841 yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);
842 }
843 if (!use_ll && (comm_rank == 0))
844 fprintf(
845 stderr, "WARNING(yac_read_scrip_basic_grid_data_): "
846 "grid \"%s\" from \"%s\": required LL but stored as GC "
847 "(>4 vertices per cell)\n",
848 grid_name, grid_filename);
849 }
850
852 (*generate_basic_grid_data_unstruct_ptr)(
853 size_t, size_t, int *, double *, double *, int *) =
854 (use_ll)?
857
858 struct yac_basic_grid_data grid_data =
859 generate_basic_grid_data_unstruct_ptr(
861 x_vertices, y_vertices, cell_to_vertex);
862
864 free(x_vertices);
865 free(y_vertices);
866 free(cell_to_vertex);
867
868 free(grid_data.core_cell_mask);
869 grid_data.core_cell_mask = cell_core_mask;
870
871 free(grid_data.vertex_ids);
872 grid_data.vertex_ids = vertex_ids;
873 free(grid_data.cell_ids);
874 grid_data.cell_ids = cell_ids;
875
876 return grid_data;
877}
878
880 char const * grid_filename, char const * mask_filename,
881 char const * grid_name, int valid_mask_value, int use_ll_edges) {
882
883 return
885 grid_filename, mask_filename, MPI_COMM_NULL,
886 grid_name, valid_mask_value, use_ll_edges,
887 NULL, NULL, NULL, NULL, NULL);
888}
889
891 char const * grid_filename, char const * mask_filename,
892 MPI_Comm comm, char const * grid_name, int valid_mask_value,
893 char const * name, int use_ll_edges, size_t * cell_coord_idx,
894 size_t ** duplicated_cell_idx, yac_int ** orig_cell_global_ids,
895 size_t * nbr_duplicated_cells) {
896
897 // generate basic grid
898 double * x_cells, * y_cells;
899 struct yac_basic_grid * basic_grid =
901 name,
903 grid_filename, mask_filename, comm, grid_name, valid_mask_value,
904 use_ll_edges, (cell_coord_idx != NULL)?&x_cells:NULL,
905 (cell_coord_idx != NULL)?&y_cells:NULL, duplicated_cell_idx,
906 orig_cell_global_ids, nbr_duplicated_cells));
907
908 if (cell_coord_idx != NULL) {
909
910 // generate 3d cell center coordinates
911 size_t num_cells =
913 yac_coordinate_pointer cell_coords =
914 xmalloc(num_cells * sizeof(*cell_coords));
915 for (size_t i = 0; i < num_cells; ++i)
916 LLtoXYZ(x_cells[i], y_cells[i], cell_coords[i]);
917 free(x_cells);
918 free(y_cells);
919
920 // register cell center coordinates in basic grid
921 *cell_coord_idx =
923 basic_grid, YAC_LOC_CELL, cell_coords);
924 }
925
926 return basic_grid;
927}
928
930 char const * grid_filename, char const * mask_filename,
931 MPI_Fint comm, char const * grid_name, int valid_mask_value,
932 char const * name, int use_ll_edges, size_t * cell_coord_idx,
933 size_t ** duplicated_cell_idx, yac_int ** orig_cell_global_ids,
934 size_t * nbr_duplicated_cells) {
935
936 return
938 grid_filename, mask_filename, MPI_Comm_f2c(comm), grid_name,
939 valid_mask_value, name, use_ll_edges, cell_coord_idx,
940 duplicated_cell_idx, orig_cell_global_ids, nbr_duplicated_cells);
941}
942
944 char const * grid_filename, char const * mask_filename,
945 char const * grid_name, int valid_mask_value, char const * name,
946 int use_ll_edges, size_t * cell_coord_idx,
947 size_t ** duplicated_cell_idx, yac_int ** orig_cell_global_ids,
948 size_t * nbr_duplicated_cells) {
949
950 return
952 grid_filename, mask_filename, MPI_COMM_NULL, grid_name, valid_mask_value,
953 name, use_ll_edges, cell_coord_idx, duplicated_cell_idx,
954 orig_cell_global_ids, nbr_duplicated_cells);
955}
956
957/* ---------------------------------------------------------------- */
958
960 const void * a,const void * b) {
961
962 int ret = ((struct cell_vertices_with_index const *)a)->num_vertices -
963 ((struct cell_vertices_with_index const *)b)->num_vertices;
964 if (ret) return ret;
965 int count = ((struct cell_vertices_with_index const *)a)->num_vertices;
966 for (int i = 0; (i < count) && !ret; ++i)
967 ret = ((struct cell_vertices_with_index const *)a)->cell_to_vertex[i] -
968 ((struct cell_vertices_with_index const *)b)->cell_to_vertex[i];
969 return ret;
970}
971
973 double ** vertex_lon, double ** vertex_lat, yac_int ** vertex_ids,
974 size_t * nbr_vertices, int * old_to_new_id) {
975
976 struct point_with_index * sort_array =
977 xmalloc(*nbr_vertices * sizeof(*sort_array));
978
979 double const scale = (double)(2 << 20);
980 int32_t const periode = (int32_t)(360.0 * scale);
981
982 for (size_t i = 0; i < *nbr_vertices; ++i) {
983
984 int32_t curr_lon = (int32_t)round((*vertex_lon)[i] * scale);
985 int32_t curr_lat = (int32_t)round((*vertex_lat)[i] * scale);
986
987 if ((curr_lat == (int32_t)(90.0 * scale)) ||
988 (curr_lat == (int32_t)(-90.0 * scale))) {
989
990 curr_lon = 0;
991
992 } else {
993 if (curr_lon <= 0)
994 curr_lon = curr_lon - (((curr_lon - periode) / periode) * periode);
995 else if (curr_lon > periode)
996 curr_lon = curr_lon - ((curr_lon / periode) * periode);
997 }
998
999 sort_array[i].lon = curr_lon;
1000 sort_array[i].lat = curr_lat;
1001 sort_array[i].dlon = (*vertex_lon)[i];
1002 sort_array[i].dlat = (*vertex_lat)[i];
1003 sort_array[i].id = (vertex_ids != NULL)?(*vertex_ids)[i]:XT_INT_MAX;
1004
1005 sort_array[i].i = i;
1006 }
1007
1008 yac_mergesort(sort_array, *nbr_vertices, sizeof(*sort_array),
1010
1011 struct point_with_index dummy =
1012 {.lon = INT32_MAX, .lat = INT32_MAX, .i = SIZE_MAX};
1013 struct point_with_index * prev = &dummy, * curr = sort_array;
1014 size_t new_nbr_vertices = 0;
1015 for (size_t i = 0; i < *nbr_vertices; ++i, ++curr) {
1016
1017 if (compare_point_with_index_coord(prev, curr)) {
1018
1019 (*vertex_lon)[new_nbr_vertices] = curr->dlon;
1020 (*vertex_lat)[new_nbr_vertices] = curr->dlat;
1021 if (vertex_ids != NULL)
1022 (*vertex_ids)[new_nbr_vertices] = curr->id;
1023 sort_array[new_nbr_vertices] = *curr;
1024 new_nbr_vertices++;
1025 prev = curr;
1026 }
1027 old_to_new_id[curr->i] = (int)(new_nbr_vertices - 1);
1028 }
1029
1030 (*vertex_lon) = xrealloc(*vertex_lon, new_nbr_vertices * sizeof(**vertex_lon));
1031 (*vertex_lat) = xrealloc(*vertex_lat, new_nbr_vertices * sizeof(**vertex_lat));
1032 if (vertex_ids != NULL)
1033 (*vertex_ids) = xrealloc(*vertex_ids, new_nbr_vertices * sizeof(**vertex_ids));
1034 *nbr_vertices = new_nbr_vertices;
1035
1036 return sort_array;
1037}
1038
1040 double ** vertex_lon, double ** vertex_lat,
1041 size_t * nbr_vertices, int * old_to_new_id) {
1042
1043 free(
1045 vertex_lon, vertex_lat, NULL, nbr_vertices, old_to_new_id));
1046}
1047
1048// djb2 hash function
1049unsigned long djb2_hash(unsigned char * values, size_t count) {
1050
1051 unsigned long hash = 5381;
1052
1053 for (size_t i = 0; i < count; ++i) {
1054 unsigned int value = values[i];
1055 hash = ((hash << 5) + hash) + value; /* hash * 33 + value */
1056 }
1057
1058 return hash;
1059}
1060
1062 void * data, size_t data_size, int comm_size) {
1063
1064 return
1065 djb2_hash((unsigned char*)data, data_size) % (unsigned long)comm_size;
1066}
1067
1068static MPI_Datatype get_point_with_index_mpi_datatype(MPI_Comm comm) {
1069
1070 struct point_with_index dummy;
1071 MPI_Datatype point_with_index_dt;
1072 int array_of_blocklengths[] = {1, 1, 1, 1, 1};
1073 const MPI_Aint array_of_displacements[] =
1074 {(MPI_Aint)(intptr_t)(const void *)&(dummy.lon) -
1075 (MPI_Aint)(intptr_t)(const void *)&dummy,
1076 (MPI_Aint)(intptr_t)(const void *)&(dummy.lat) -
1077 (MPI_Aint)(intptr_t)(const void *)&dummy,
1078 (MPI_Aint)(intptr_t)(const void *)&(dummy.id) -
1079 (MPI_Aint)(intptr_t)(const void *)&dummy,
1080 (MPI_Aint)(intptr_t)(const void *)&(dummy.dlon) -
1081 (MPI_Aint)(intptr_t)(const void *)&dummy,
1082 (MPI_Aint)(intptr_t)(const void *)&(dummy.dlat) -
1083 (MPI_Aint)(intptr_t)(const void *)&dummy};
1084 const MPI_Datatype array_of_types[] =
1085 {MPI_INT32_T, MPI_INT32_T, yac_int_dt, MPI_DOUBLE, MPI_DOUBLE};
1087 MPI_Type_create_struct(5, array_of_blocklengths, array_of_displacements,
1088 array_of_types, &point_with_index_dt), comm);
1089 return yac_create_resized(point_with_index_dt, sizeof(dummy), comm);
1090}
1091
1093 double ** vertex_lon, double ** vertex_lat,
1094 size_t * nbr_vertices_, int * old_to_new_id, yac_int ** vertex_ids,
1095 MPI_Comm comm) {
1096
1097 // remove the locally duplicated vertices
1098 struct point_with_index * sort_array =
1100 vertex_lon, vertex_lat, vertex_ids, nbr_vertices_, old_to_new_id);
1101
1102 sort_array = xrealloc(sort_array, *nbr_vertices_ * sizeof(*sort_array));
1103
1104 size_t nbr_vertices = *nbr_vertices_;
1105
1106 size_t * sendcounts, * recvcounts, * sdispls, * rdispls;
1108 1, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
1109
1110 int comm_rank, comm_size;
1111 yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);
1112 yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
1113
1114 // set up counts and displs for redistribution of points
1115 for (size_t i = 0; i < nbr_vertices; ++i) {
1116 int rank =
1118 &(sort_array[i].lon), 2 * sizeof(sort_array[i].lon), comm_size);
1119 sendcounts[rank]++;
1120 sort_array[i].rank = rank;
1121 sort_array[i].i = i;
1122 }
1123
1124 // sort points by bucket ranks
1125 yac_mergesort(sort_array, nbr_vertices, sizeof(*sort_array),
1127
1129 1, sendcounts, recvcounts, sdispls, rdispls, comm);
1130
1131 size_t recv_count = recvcounts[comm_size - 1] + rdispls[comm_size - 1];
1132 struct point_with_index * dist_sort_array =
1133 xmalloc(recv_count * sizeof(*dist_sort_array));
1134
1135 // redistribute point information
1136 MPI_Datatype point_with_index_dt =
1138 yac_mpi_call(MPI_Type_commit(&point_with_index_dt), comm);
1140 sort_array, sendcounts, sdispls+1,
1141 dist_sort_array, recvcounts, rdispls,
1142 sizeof(*sort_array), point_with_index_dt, comm);
1143
1144 for (size_t i = 0; i < recv_count; ++i) dist_sort_array[i].i = i;
1145
1146 // sort received point data based on coordinates
1147 yac_mergesort(dist_sort_array, recv_count, sizeof(*dist_sort_array),
1149
1150 struct point_with_index dummy =
1151 {.lon = INT32_MAX, .lat = INT32_MAX};
1152 struct point_with_index * prev = &dummy, * curr = dist_sort_array;
1153 for (size_t i = 0; i < recv_count; ++i, ++curr) {
1154 if (compare_point_with_index_coord(prev, curr)) {
1155
1156 prev = curr;
1157
1158 } else {
1159
1160 curr->id = prev->id;
1161 curr->dlon = prev->dlon;
1162 curr->dlat = prev->dlat;
1163 }
1164 }
1165
1166 // bring array into original order
1167 yac_mergesort(dist_sort_array, recv_count, sizeof(*dist_sort_array),
1169
1170 // redistribute point information
1172 dist_sort_array, recvcounts, rdispls,
1173 sort_array, sendcounts, sdispls+1,
1174 sizeof(*sort_array), point_with_index_dt, comm);
1175 free(dist_sort_array);
1176 yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
1177 yac_mpi_call(MPI_Type_free(&point_with_index_dt), comm);
1178
1179 // bring vertex indices into original order
1180 for (size_t i = 0; i < nbr_vertices; ++i) {
1181
1182 (*vertex_lon)[sort_array[i].i] = sort_array[i].dlon;
1183 (*vertex_lat)[sort_array[i].i] = sort_array[i].dlat;
1184 (*vertex_ids)[sort_array[i].i] = sort_array[i].id;
1185 }
1186
1187 free(sort_array);
1188}
1189
1191 int * cell_to_vertex, int * num_vertices_per_cell, int * cell_mask,
1192 size_t nbr_cells, size_t ** duplicated_cell_idx,
1193 size_t ** orig_cell_idx, size_t * nbr_duplicated_cells) {
1194
1195 // count number of unmasked cells
1196 size_t nbr_unmasked_cells = 0;
1197 for (size_t i = 0; i < nbr_cells; ++i)
1198 if (cell_mask[i]) ++nbr_unmasked_cells;
1199
1200 struct cell_vertices_with_index * sort_array =
1201 xmalloc(nbr_unmasked_cells * sizeof(*sort_array));
1202
1203 // get all unmasked cells
1204 for (size_t i = 0, offset = 0, j = 0; i < nbr_cells; ++i) {
1205
1206 if (cell_mask[i]) {
1207
1208 sort_array[j].num_vertices = num_vertices_per_cell[i];
1209 sort_array[j].cell_to_vertex = cell_to_vertex + offset;
1210 sort_array[j].i = i;
1211 ++j;
1212 }
1213
1214 offset += (size_t)(num_vertices_per_cell[i]);
1215 }
1216
1217 // sort cells
1218 yac_mergesort(sort_array, nbr_unmasked_cells, sizeof(*sort_array),
1220
1221 // count number of duplicated cells
1222 *nbr_duplicated_cells = 0;
1223 for (size_t i = 1; i < nbr_unmasked_cells; ++i)
1224 if (!compare_cell_vertices_with_index(sort_array + (i-1), sort_array + i))
1225 ++*nbr_duplicated_cells;
1226
1227 // get duplicated cells
1228 *duplicated_cell_idx =
1229 xmalloc(*nbr_duplicated_cells * sizeof(**duplicated_cell_idx));
1230 *orig_cell_idx =
1231 xmalloc(*nbr_duplicated_cells * sizeof(**orig_cell_idx));
1232 struct cell_vertices_with_index *prev = sort_array, *curr = sort_array + 1;
1233 for (size_t i = 1, j = 0; i < nbr_unmasked_cells; ++i, ++curr) {
1234 if (compare_cell_vertices_with_index(prev, curr)) {
1235 prev = curr;
1236 } else {
1237 (*duplicated_cell_idx)[j] = (size_t)(curr->i);
1238 (*orig_cell_idx)[j] = (size_t)(prev->i);
1239 ++j;
1240 }
1241 }
1242
1243 free(sort_array);
1244}
1245
1247 MPI_Comm comm, size_t count) {
1248
1249 MPI_Datatype cell_to_vertex_ids_dt;
1251 MPI_Type_contiguous(
1252 (int)count, yac_int_dt, &cell_to_vertex_ids_dt), comm);
1253 return cell_to_vertex_ids_dt;
1254}
1255
1257 int * cell_to_vertex, int * num_vertices_per_cell, int * cell_mask,
1258 yac_int * cell_ids, size_t num_cells, yac_int * vertex_ids,
1259 MPI_Comm comm, size_t ** duplicated_cell_idx_,
1260 yac_int ** orig_cell_ids_, size_t * num_duplicated_cells_) {
1261
1262 int comm_size;
1263 yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
1264
1265 // determine the local maximum number of vertices per cell and count
1266 // number of unmasked cells
1267 size_t max_num_vertices_per_cell = 0;
1268 size_t local_num_cells = 0;
1269 for (size_t i = 0; i < num_cells; ++i) {
1270 if (cell_mask[i]) {
1271 if ((size_t)(num_vertices_per_cell[i]) > max_num_vertices_per_cell)
1272 max_num_vertices_per_cell = (size_t)(num_vertices_per_cell[i]);
1273 local_num_cells++;
1274 }
1275 }
1276
1277 // determine the global maximum number of vertices per cell
1279 MPI_Allreduce(
1280 MPI_IN_PLACE, &max_num_vertices_per_cell, 1, YAC_MPI_SIZE_T, MPI_MAX,
1281 comm), comm);
1282
1283 // generate local cell to vertex ids array
1284 yac_int * local_cell_ids =
1285 xmalloc(local_num_cells * (1 + max_num_vertices_per_cell) *
1286 sizeof(*local_cell_ids));
1287 yac_int * local_cell_to_vertex_ids = local_cell_ids + local_num_cells;
1288
1289 // fill local cell to vertex ids array
1290 yac_int * curr_local_cell_ids = local_cell_ids;
1291 yac_int * curr_local_cell_to_vertex_ids = local_cell_to_vertex_ids;
1292 int * curr_cell_to_vertex = cell_to_vertex;
1293 for (size_t i = 0; i < num_cells; ++i) {
1294
1295 size_t curr_num_vertices_per_cell = (size_t)(num_vertices_per_cell[i]);
1296 if (cell_mask[i]) {
1297
1298 *curr_local_cell_ids = cell_ids[i];
1299
1300 // get global vertex ids for current cell
1301 size_t j = 0;
1302 for (; j < curr_num_vertices_per_cell; ++j)
1303 curr_local_cell_to_vertex_ids[j] =
1304 vertex_ids[curr_cell_to_vertex[j]];
1305
1306 // sort global vertex ids
1307 qsort(curr_local_cell_to_vertex_ids, curr_num_vertices_per_cell,
1308 sizeof(*curr_local_cell_to_vertex_ids), compare_yac_int);
1309
1310 // fill remaining entries
1311 for (; j < max_num_vertices_per_cell; ++j)
1312 curr_local_cell_to_vertex_ids[j] = XT_INT_MAX;
1313
1314 curr_local_cell_ids++;
1315 curr_local_cell_to_vertex_ids += max_num_vertices_per_cell;
1316 }
1317 curr_cell_to_vertex += curr_num_vertices_per_cell;
1318 }
1319
1320 size_t * sendcounts, * recvcounts, * sdispls, * rdispls;
1322 1, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
1323
1324 // determine distributed owners for all cells
1325 int * dist_cell_ranks = xmalloc(local_num_cells * sizeof(*dist_cell_ranks));
1326 for (size_t i = 0; i < local_num_cells; ++i) {
1327 int rank =
1329 local_cell_to_vertex_ids + i * max_num_vertices_per_cell,
1330 max_num_vertices_per_cell * sizeof(*local_cell_to_vertex_ids),
1331 comm_size);
1332 sendcounts[rank]++;
1333 dist_cell_ranks[i] = rank;
1334 }
1335
1337 1, sendcounts, recvcounts, sdispls, rdispls, comm);
1338
1339 size_t dist_num_cells = recvcounts[comm_size-1] + rdispls[comm_size-1];
1340 yac_int * yac_int_buffer =
1341 xmalloc((local_num_cells + dist_num_cells) *
1342 (1 + max_num_vertices_per_cell) * sizeof(*yac_int_buffer));
1343 yac_int * dist_cell_ids = yac_int_buffer;
1344 yac_int * dist_cell_vertex_ids = yac_int_buffer + dist_num_cells;
1345 yac_int * send_local_cell_ids =
1346 yac_int_buffer + dist_num_cells * (1 + max_num_vertices_per_cell);
1347 yac_int * send_local_cell_to_vertex_ids =
1348 yac_int_buffer + local_num_cells +
1349 dist_num_cells * (1 + max_num_vertices_per_cell);
1350
1351 // pack send data
1352 for (size_t i = 0; i < local_num_cells; ++i) {
1353 size_t pos = sdispls[dist_cell_ranks[i] + 1]++;
1354 send_local_cell_ids[pos] = local_cell_ids[i];
1355 memcpy(
1356 send_local_cell_to_vertex_ids + pos * max_num_vertices_per_cell,
1357 local_cell_to_vertex_ids + i * max_num_vertices_per_cell,
1358 max_num_vertices_per_cell * sizeof(*send_local_cell_to_vertex_ids));
1359 }
1360 local_cell_ids =
1361 xrealloc(local_cell_ids, local_num_cells * sizeof(*local_cell_ids));
1362
1363 // redistribute cell ids and cell vertex ids
1364 MPI_Datatype cell_to_vertex_ids_dt =
1365 get_cell_to_vertex_ids_mpi_datatype(comm, max_num_vertices_per_cell);
1366 yac_mpi_call(MPI_Type_commit(&cell_to_vertex_ids_dt), comm);
1367 yac_alltoallv_yac_int_p2p(
1368 send_local_cell_ids, sendcounts, sdispls,
1369 dist_cell_ids, recvcounts, rdispls, comm);
1371 send_local_cell_to_vertex_ids, sendcounts, sdispls,
1372 dist_cell_vertex_ids, recvcounts, rdispls,
1373 max_num_vertices_per_cell * sizeof(*dist_cell_vertex_ids),
1374 cell_to_vertex_ids_dt, comm);
1375 yac_mpi_call(MPI_Type_free(&cell_to_vertex_ids_dt), comm);
1376
1377 // free some memory
1378 yac_int_buffer =
1379 xrealloc(
1380 yac_int_buffer,
1381 dist_num_cells * (1 + max_num_vertices_per_cell) *
1382 sizeof(*yac_int_buffer));
1383 dist_cell_ids = yac_int_buffer;
1384 dist_cell_vertex_ids = yac_int_buffer + dist_num_cells;
1385
1386 // generate sortable data structure
1387 struct cell_vertex_ids_with_index * sort_array =
1388 xmalloc(dist_num_cells * sizeof(*sort_array));
1389 for (size_t i = 0; i < dist_num_cells; ++i) {
1390 sort_array[i].cell_id = dist_cell_ids[i];
1391 sort_array[i].vertex_ids =
1392 dist_cell_vertex_ids + i * max_num_vertices_per_cell;
1393 sort_array[i].num_vertices = max_num_vertices_per_cell;
1394 sort_array[i].i = i;
1395 }
1396
1397 // sort cells by vertex ids
1399 sort_array, dist_num_cells, sizeof(*sort_array),
1401
1402 // determine duplicated cells
1403 struct cell_vertex_ids_with_index * prev = sort_array;
1404 struct cell_vertex_ids_with_index * curr = sort_array + 1;
1405 for (size_t i = 1; i < dist_num_cells; ++i, ++curr) {
1406
1408 prev = curr;
1409 } else {
1410 curr->cell_id = prev->cell_id;
1411 }
1412 }
1413
1414 // pack cell ids
1415 yac_int_buffer =
1416 xrealloc(
1417 yac_int_buffer,
1418 (local_num_cells + dist_num_cells) * sizeof(*yac_int_buffer));
1419 yac_int * recv_local_cell_ids = yac_int_buffer;
1420 dist_cell_ids = yac_int_buffer + local_num_cells;
1421 for (size_t i = 0; i < dist_num_cells; ++i)
1422 dist_cell_ids[sort_array[i].i] = sort_array[i].cell_id;
1423 free(sort_array);
1424
1425 // redistribute cell ids
1426 yac_alltoallv_yac_int_p2p(
1427 dist_cell_ids, recvcounts, rdispls,
1428 recv_local_cell_ids, sendcounts, sdispls, comm);
1429 yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
1430
1431 // bring received local cell ids into original order and
1432 // count number of duplicated cells
1433 size_t num_duplicated_cells = 0;
1434 for (size_t i = 0; i < local_num_cells; ++i) {
1435
1436 size_t pos = sdispls[dist_cell_ranks[i]]++;
1437
1438 if (local_cell_ids[i] != recv_local_cell_ids[pos]) {
1439 num_duplicated_cells++;
1440 local_cell_ids[i] = recv_local_cell_ids[pos];
1441 } else {
1442 local_cell_ids[i] = XT_INT_MAX;
1443 }
1444 }
1445 free(yac_int_buffer);
1446 free(dist_cell_ranks);
1447
1448 size_t * duplicated_cell_idx =
1449 xmalloc(num_duplicated_cells * sizeof(*duplicated_cell_idx));
1450 yac_int * orig_cell_ids =
1451 xmalloc(num_duplicated_cells * sizeof(*orig_cell_ids));
1452
1453 // get duplicated cells
1454 for (size_t i = 0, j = 0, k = 0; i < num_cells; ++i) {
1455
1456 if (cell_mask[i]) {
1457
1458 if (local_cell_ids[j] != XT_INT_MAX) {
1459 duplicated_cell_idx[k] = i;
1460 orig_cell_ids[k] = local_cell_ids[j];
1461 ++k;
1462 }
1463 ++j;
1464 }
1465 }
1466 free(local_cell_ids);
1467
1468 *duplicated_cell_idx_ = duplicated_cell_idx;
1469 *orig_cell_ids_ = orig_cell_ids;
1470 *num_duplicated_cells_ = num_duplicated_cells;
1471}
1472
1473#else
1474
1476 char const * grid_filename, char const * mask_filename,
1477 char const * grid_name, int valid_mask_value, int use_ll_edges) {
1478
1479 UNUSED(grid_filename);
1480 UNUSED(mask_filename);
1481 UNUSED(grid_name);
1482 UNUSED(valid_mask_value);
1483 UNUSED(use_ll_edges);
1484 die(
1485 "ERROR(yac_read_scrip_basic_grid_data): "
1486 "YAC is built without the NetCDF support");
1487
1488 return
1490 (size_t[]){0,0}, (int[]){0,0}, NULL, NULL);
1491}
1492
1494 char const * grid_filename, char const * mask_filename,
1495 char const * grid_name, int valid_mask_value, char const * name,
1496 int use_ll_edges, size_t * cell_coord_idx,
1497 size_t ** duplicated_cell_idx, yac_int ** orig_cell_global_ids,
1498 size_t * nbr_duplicated_cells) {
1499
1500 UNUSED(grid_filename);
1501 UNUSED(mask_filename);
1502 UNUSED(grid_name);
1503 UNUSED(valid_mask_value);
1504 UNUSED(name);
1505 UNUSED(use_ll_edges);
1506 UNUSED(cell_coord_idx);
1507 UNUSED(duplicated_cell_idx);
1508 UNUSED(orig_cell_global_ids);
1509 UNUSED(nbr_duplicated_cells);
1510 die(
1511 "ERROR(yac_read_scrip_basic_grid): "
1512 "YAC is built without the NetCDF support");
1513
1514 return NULL;
1515}
1516
1518 char const * grid_filename, char const * mask_filename,
1519 MPI_Comm comm, char const * grid_name, int valid_mask_value,
1520 char const * name, int use_ll_edges, size_t * cell_coord_idx,
1521 size_t ** duplicated_cell_idx, yac_int ** orig_cell_global_ids,
1522 size_t * nbr_duplicated_cells) {
1523
1524 UNUSED(grid_filename);
1525 UNUSED(mask_filename);
1526 UNUSED(comm);
1527 UNUSED(grid_name);
1528 UNUSED(valid_mask_value);
1529 UNUSED(name);
1530 UNUSED(use_ll_edges);
1531 UNUSED(cell_coord_idx);
1532 UNUSED(duplicated_cell_idx);
1533 UNUSED(orig_cell_global_ids);
1534 UNUSED(nbr_duplicated_cells);
1535 die(
1536 "ERROR(yac_read_scrip_basic_grid_parallel): "
1537 "YAC is built without the NetCDF support");
1538
1539 return NULL;
1540}
1541
1543 char const * grid_filename, char const * mask_filename,
1544 char const * grid_name, int valid_mask_value,
1545 size_t * num_vertices, size_t * num_cells, int ** num_vertices_per_cell,
1546 double ** x_vertices, double ** y_vertices,
1547 double ** x_cells, double ** y_cells,
1548 int ** cell_to_vertex, int ** cell_core_mask, size_t ** duplicated_cell_idx,
1549 size_t ** orig_cell_idx, size_t * nbr_duplicated_cells) {
1550
1551 UNUSED(grid_filename);
1552 UNUSED(mask_filename);
1553 UNUSED(grid_name);
1554 UNUSED(valid_mask_value);
1555 UNUSED(num_vertices);
1556 UNUSED(num_cells);
1557 UNUSED(num_vertices_per_cell);
1558 UNUSED(x_vertices);
1559 UNUSED(y_vertices);
1560 UNUSED(x_cells);
1561 UNUSED(y_cells);
1562 UNUSED(cell_to_vertex);
1563 UNUSED(cell_core_mask);
1564 UNUSED(duplicated_cell_idx);
1565 UNUSED(orig_cell_idx);
1566 UNUSED(nbr_duplicated_cells);
1567 die(
1568 "ERROR(yac_read_scrip_grid_information): "
1569 "YAC is built without the NetCDF support");
1570}
1571
1572#endif // YAC_NETCDF_ENABLED
struct yac_basic_grid * yac_basic_grid_new(char const *name, struct yac_basic_grid_data grid_data)
Definition basic_grid.c:45
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:203
size_t yac_basic_grid_get_data_size(struct yac_basic_grid *grid, enum yac_location location)
Definition basic_grid.c:142
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_ll(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell, double *x_vertices, double *y_vertices, int *cell_to_vertex)
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:71
#define YAC_RAD
Definition geometry.h:30
static void LLtoXYZ(double lon, double lat, double p_out[])
Definition geometry.h:287
struct @8::@9 value
void yac_get_io_ranks(MPI_Comm comm, int *local_is_io_, int **io_ranks_, int *num_io_ranks_)
Definition io_utils.c:303
void yac_nc_inq_varid(int ncid, char const *name, int *varidp)
Definition io_utils.c:405
void yac_nc_open(const char *path, int omode, int *ncidp)
Definition io_utils.c:344
#define YAC_HANDLE_ERROR(exp)
Definition io_utils.h:30
@ YAC_LOC_CELL
Definition location.h:14
void yac_mergesort(void *base, size_t num, size_t size, int(*compar)(const void *, const void *))
Definition mergesort.c:134
#define xrealloc(ptr, size)
Definition ppm_xfuncs.h:67
#define xmalloc(size)
Definition ppm_xfuncs.h:66
struct yac_basic_grid * yac_read_scrip_basic_grid_parallel(char const *grid_filename, char const *mask_filename, MPI_Comm comm, char const *grid_name, int valid_mask_value, char const *name, int use_ll_edges, size_t *cell_coord_idx, size_t **duplicated_cell_idx, yac_int **orig_cell_global_ids, size_t *nbr_duplicated_cells)
static void yac_read_part_scrip_grid_information(char const *grid_filename, char const *mask_filename, MPI_Comm comm, char const *grid_name, int valid_mask_value, size_t *num_vertices, size_t *num_cells, int **num_vertices_per_cell, double **x_vertices, double **y_vertices, yac_int **vertex_ids, double **x_cells, double **y_cells, yac_int **cell_ids, int **cell_to_vertex, int **cell_core_mask, size_t **duplicated_cell_idx, yac_int **orig_cell_ids, size_t *nbr_duplicated_cells)
unsigned long djb2_hash(unsigned char *values, size_t count)
static void yac_read_part_scrip_basic_grid_information(const char *filename, const char *grid_name, size_t cell_mask_size, int *cell_mask, size_t *num_vertices_, size_t *num_cells_, int **num_vertices_per_cell_, int **cell_to_vertex_, double **x_vertices_, double **y_vertices_, yac_int **vertex_ids_, double **x_cells_, double **y_cells_, yac_int **cell_ids_, size_t **duplicated_cell_idx_, yac_int **orig_cell_ids_, size_t *nbr_duplicated_cells_, int io_rank_idx, int num_io_ranks, MPI_Comm comm)
struct yac_basic_grid_data yac_read_scrip_basic_grid_data(char const *grid_filename, char const *mask_filename, char const *grid_name, int valid_mask_value, int use_ll_edges)
static int compare_point_with_index_coord_id(const void *a_, const void *b_)
static void yac_read_part_scrip_mask_information(const char *filename, const char *grid_name, size_t *num_cells_, int **cell_mask, int io_rank_idx, int num_io_ranks)
void yac_read_scrip_grid_information(char const *grid_filename, char const *mask_filename, char const *grid_name, int valid_mask_value, size_t *num_vertices, size_t *num_cells, int **num_vertices_per_cell, double **x_vertices, double **y_vertices, double **x_cells, double **y_cells, int **cell_to_vertex, int **cell_core_mask, size_t **duplicated_cell_idx, size_t **orig_cell_idx, size_t *nbr_duplicated_cells)
static void yac_read_scrip_basic_grid_information(const char *filename, const char *grid_name, size_t cell_mask_size, int *cell_mask, size_t *num_vertices_, size_t *num_cells_, int **num_vertices_per_cell_, int **cell_to_vertex_, double **x_vertices, double **y_vertices, double **x_cells, double **y_cells, size_t **duplicated_cell_idx_, size_t **orig_cell_idx_, size_t *nbr_duplicated_cells_)
static void yac_read_scrip_mask_information(const char *filename, const char *grid_name, size_t *num_cells_, int **cell_mask)
static int compare_cell_vertices_with_index(const void *a, const void *b)
static void detect_duplicated_cells_parallel(int *cell_to_vertex, int *num_vertices_per_cell, int *cell_mask, yac_int *cell_ids, size_t num_cells, yac_int *vertex_ids, MPI_Comm comm, size_t **duplicated_cell_idx, yac_int **orig_cell_ids, size_t *num_duplicated_cells_)
static MPI_Datatype get_point_with_index_mpi_datatype(MPI_Comm comm)
static void remove_duplicated_vertices_parallel(double **vertex_lon, double **vertex_lat, size_t *nbr_vertices, int *old_to_new_id, yac_int **vertex_ids, MPI_Comm comm)
static struct point_with_index * remove_duplicated_vertices_(double **vertex_lon, double **vertex_lat, yac_int **vertex_ids, size_t *nbr_vertices, int *old_to_new_id)
static int compare_cell_vertex_ids_with_index_ids_id(const void *a_, const void *b_)
static int compute_bucket(void *data, size_t data_size, int comm_size)
struct yac_basic_grid * yac_read_scrip_basic_grid_parallel_f2c(char const *grid_filename, char const *mask_filename, MPI_Fint comm, char const *grid_name, int valid_mask_value, char const *name, int use_ll_edges, size_t *cell_coord_idx, size_t **duplicated_cell_idx, yac_int **orig_cell_global_ids, size_t *nbr_duplicated_cells)
static struct yac_basic_grid_data yac_read_scrip_basic_grid_data_(char const *grid_filename, char const *mask_filename, MPI_Comm comm, char const *grid_name, int valid_mask_value, int use_ll, double **x_cells, double **y_cells, size_t **duplicated_cell_idx, yac_int **orig_cell_global_ids, size_t *nbr_duplicated_cells)
static int compare_yac_int(const void *a, const void *b)
static int compare_point_with_index_rank(const void *a_, const void *b_)
static void detect_duplicated_cells(int *cell_to_vertex, int *num_vertices_per_cell, int *cell_mask, size_t nbr_cells, size_t **duplicated_cell_idx, size_t **orig_cell_idx, size_t *nbr_duplicated_cells)
static MPI_Datatype get_cell_to_vertex_ids_mpi_datatype(MPI_Comm comm, size_t count)
static int compare_cell_vertex_ids_with_index_ids(const void *a_, const void *b_)
static int compare_point_with_index_i(const void *a_, const void *b_)
static int compare_int(const void *a, const void *b)
static void remove_duplicated_vertices(double **vertex_lon, double **vertex_lat, size_t *nbr_vertices, int *old_to_new_id)
struct yac_basic_grid * yac_read_scrip_basic_grid(char const *grid_filename, char const *mask_filename, char const *grid_name, int valid_mask_value, char const *name, int use_ll_edges, size_t *cell_coord_idx, size_t **duplicated_cell_idx, yac_int **orig_cell_global_ids, size_t *nbr_duplicated_cells)
static int compare_point_with_index_coord(const void *a_, const void *b_)
void yac_quicksort_index_size_t_int(size_t *a, size_t n, int *idx)
#define YAC_ASSERT_F(exp, format,...)
Definition yac_assert.h:19
#define die(msg)
Definition yac_assert.h:12
#define YAC_ASSERT(exp, msg)
Definition yac_assert.h:16
void yac_generate_alltoallv_args(int count, size_t const *sendcounts, size_t *recvcounts, size_t *sdispls, size_t *rdispls, MPI_Comm comm)
Definition yac_mpi.c:571
void yac_free_comm_buffers(size_t *sendcounts, size_t *recvcounts, size_t *sdispls, size_t *rdispls)
Definition yac_mpi.c:627
void yac_get_comm_buffers(int count, size_t **sendcounts, size_t **recvcounts, size_t **sdispls, size_t **rdispls, MPI_Comm comm)
Definition yac_mpi.c:596
MPI_Datatype yac_create_resized(MPI_Datatype dt, size_t new_size, MPI_Comm comm)
Definition yac_mpi.c:550
void yac_alltoallv_p2p(void const *send_buffer, size_t const *sendcounts, size_t const *sdispls, void *recv_buffer, size_t const *recvcounts, size_t const *rdispls, size_t dt_size, MPI_Datatype dt, MPI_Comm comm)
Definition yac_mpi.c:131
#define yac_mpi_call(call, comm)
#define YAC_MPI_SIZE_T
Xt_int yac_int
Definition yac_types.h:15
#define yac_int_dt
Definition yac_types.h:16
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19