YAC: src/core/basic_grid.c Source File

17struct yac_basic_grid {

18 char * name;

19 int is_empty;

20 struct yac_field_data_set * field_data_set;

21 struct yac_basic_grid_data data;

22};

24static struct yac_basic_grid_data yac_basic_grid_data_empty = {

25 .vertex_coordinates = NULL,

26 .cell_ids = NULL,

27 .vertex_ids = NULL,

28 .edge_ids = NULL,

29 .num_cells = 0,

30 .num_vertices = 0,

31 .num_edges = 0,

32 .core_cell_mask = NULL,

33 .core_vertex_mask = NULL,

34 .core_edge_mask = NULL,

35 .num_vertices_per_cell = NULL,

36 .num_cells_per_vertex = NULL,

37 .cell_to_vertex = NULL,

38 .cell_to_vertex_offsets = NULL,

39 .cell_to_edge = NULL,

40 .cell_to_edge_offsets = NULL,

41 .vertex_to_cell = NULL,

42 .vertex_to_cell_offsets = NULL,

43 .edge_to_vertex = NULL,

44 .edge_type = NULL,

45 .num_total_cells = 0,

46 .num_total_vertices = 0,

47 .num_total_edges = 0

48};

50struct yac_basic_grid * yac_basic_grid_new(

51 char const * name, struct yac_basic_grid_data grid_data) {

52

53 struct yac_basic_grid * grid = xmalloc(1 * sizeof(*grid));

54

55 grid->name = strdup(name);

56 grid->is_empty = 0;

57 grid->field_data_set = yac_field_data_set_empty_new();

58 grid->data = grid_data;

59

60 return grid;

61}

63struct yac_basic_grid * yac_basic_grid_empty_new(char const * name) {

64 struct yac_basic_grid * grid =

65 yac_basic_grid_new(name, yac_basic_grid_data_empty);

66 grid->is_empty = 1;

67 return grid;

68}

70void yac_basic_grid_delete(struct yac_basic_grid * grid) {

71

72 YAC_ASSERT(

73 grid, "ERROR(yac_basic_grid_delete): "

74 "NULL is not a valid value for argument grid")

75 free(grid->name);

76 yac_field_data_set_delete(grid->field_data_set);

77 yac_basic_grid_data_free(grid->data);

78 free(grid);

79}

81yac_const_coordinate_pointer yac_basic_grid_get_field_coordinates(

82 struct yac_basic_grid * grid, struct yac_interp_field field) {

83

84 YAC_ASSERT(

85 grid, "ERROR(yac_basic_grid_get_field_coordinates): "

86 "NULL is not a valid value for argument grid")

87

88 yac_const_coordinate_pointer coords =

89 (field.coordinates_idx != SIZE_MAX)?

90 yac_field_data_get_coordinates_data(

91 yac_basic_grid_get_field_data(grid, field.location),

92 field.coordinates_idx):NULL;

93

94 // if no field coordinates are defined, but the location is at the corners of

95 // of the grid cells, return coordinates of them

96 return

97 ((coords != NULL) || (field.location != YAC_LOC_CORNER))?

98 coords:((yac_const_coordinate_pointer)(grid->data.vertex_coordinates));

99}

101int const * yac_basic_grid_get_core_mask(

102 struct yac_basic_grid * grid, enum yac_location location) {

103

104 YAC_ASSERT(

105 (location == YAC_LOC_CELL) ||

106 (location == YAC_LOC_CORNER) ||

107 (location == YAC_LOC_EDGE),

108 "ERROR(yac_basic_grid_get_core_mask): invalid location")

109

110 switch (location) {

111 default:

112 case(YAC_LOC_CELL): return grid->data.core_cell_mask;

113 case(YAC_LOC_CORNER): return grid->data.core_vertex_mask;

114 case(YAC_LOC_EDGE): return grid->data.core_edge_mask;

115 };

116}

118int const * yac_basic_grid_get_field_mask(

119 struct yac_basic_grid * grid, struct yac_interp_field field) {

120

121 if (field.masks_idx == SIZE_MAX) return NULL;

122

123 return

124 yac_field_data_get_mask_data(

125 yac_basic_grid_get_field_data(grid, field.location), field.masks_idx);

126}

128char const * yac_basic_grid_get_name(struct yac_basic_grid * grid) {

129

130 YAC_ASSERT(

131 grid, "ERROR(yac_basic_grid_get_name): "

132 "NULL is not a valid value for argument grid")

133

134 return grid->name;

135}

137struct yac_basic_grid_data * yac_basic_grid_get_data(

138 struct yac_basic_grid * grid) {

139

140 YAC_ASSERT(

141 grid, "ERROR(yac_basic_grid_get_data): "

142 "NULL is not a valid value for argument grid")

143

144 return &(grid->data);

145}

147size_t yac_basic_grid_get_data_size(

148 struct yac_basic_grid * grid, enum yac_location location) {

149

150 YAC_ASSERT(

151 grid, "ERROR(yac_basic_grid_get_data_size): "

152 "NULL is not a valid value for argument grid")

153 YAC_ASSERT(

154 (location == YAC_LOC_CELL) ||

155 (location == YAC_LOC_CORNER) ||

156 (location == YAC_LOC_EDGE),

157 "ERROR(yac_basic_grid_get_data_size): invalid location")

158

159 switch (location) {

160 default:

161 case (YAC_LOC_CELL):

162 return grid->data.num_cells;

163 case (YAC_LOC_CORNER):

164 return grid->data.num_vertices;

165 case (YAC_LOC_EDGE):

166 return grid->data.num_edges;

167 };

168}

170size_t yac_basic_grid_get_data_size_f2c(

171 struct yac_basic_grid * grid, int location) {

172

173 return

174 yac_basic_grid_get_data_size(grid, yac_get_location(location));

175}

177size_t yac_basic_grid_get_named_mask_idx(

178 struct yac_basic_grid * grid, enum yac_location location,

179 char const * mask_name) {

180

181 if (mask_name == NULL) return SIZE_MAX;

182

183 struct yac_field_data * data =

184 yac_basic_grid_get_field_data(grid, location);

185

186 if (data == NULL) return SIZE_MAX;

187

188 size_t mask_idx = SIZE_MAX;

189 size_t masks_count = yac_field_data_get_masks_count(data);

190

191 for (size_t i = 0; (i < masks_count) && (mask_idx == SIZE_MAX); ++i) {

192 char const * curr_mask_name =

193 yac_field_data_get_mask_name(data, i);

194 if ((curr_mask_name != NULL) &&

195 (!strcmp(mask_name, curr_mask_name)))

196 mask_idx = i;

197 }

198

199 YAC_ASSERT_F(

200 mask_idx != SIZE_MAX,

201 "ERROR(yac_basic_grid_get_named_mask_idx): grid \"%s\" does not contain "

202 "%s-mask with the name \"%s\"",

203 grid->name, yac_loc2str(location), mask_name)

204

205 return mask_idx;

206}

208size_t yac_basic_grid_add_coordinates_nocpy(

209 struct yac_basic_grid * grid,

210 enum yac_location location, yac_coordinate_pointer coordinates) {

211

212 YAC_ASSERT(

213 grid, "ERROR(yac_basic_grid_add_coordinates_nocpy): "

214 "NULL is not a valid value for argument grid")

215 YAC_ASSERT_F(

216 !grid->is_empty, "ERROR(yac_basic_grid_add_coordinates_nocpy): "

217 "grid \"%s\" is an empty grid", grid->name)

218

219 return

220 yac_field_data_set_add_coordinates_nocpy(

221 grid->field_data_set, location, coordinates);

222}

224size_t yac_basic_grid_add_coordinates_nocpy_f2c(

225 struct yac_basic_grid * grid, int location, double * coordinates) {

226

227 return

228 yac_basic_grid_add_coordinates_nocpy(

229 grid, yac_get_location(location), (yac_coordinate_pointer)coordinates);

230}

232size_t yac_basic_grid_add_coordinates(

233 struct yac_basic_grid * grid, enum yac_location location,

234 yac_coordinate_pointer coordinates, size_t count) {

235

236 YAC_ASSERT(

237 grid, "ERROR(yac_basic_grid_add_coordinates): "

238 "NULL is not a valid value for argument grid")

239 YAC_ASSERT_F(

240 !grid->is_empty, "ERROR(yac_basic_grid_add_coordinates): "

241 "grid \"%s\" is an empty grid", grid->name)

242

243 return

244 yac_field_data_set_add_coordinates(

245 grid->field_data_set, location, coordinates, count);

246}

248size_t yac_basic_grid_add_coordinates_f2c(

249 struct yac_basic_grid * grid, int location,

250 double * coordinates, size_t count) {

251

252 return

253 yac_basic_grid_add_coordinates(

254 grid, yac_get_location(location),

255 (yac_coordinate_pointer)coordinates, count);

256}

258size_t yac_basic_grid_add_mask_nocpy(

259 struct yac_basic_grid * grid,

260 enum yac_location location, int const * mask,

261 char const * mask_name) {

262

263 YAC_ASSERT(

264 grid, "ERROR(yac_basic_grid_add_mask_nocpy): "

265 "NULL is not a valid value for argument grid")

266 YAC_ASSERT_F(

267 !grid->is_empty, "ERROR(yac_basic_grid_add_mask_nocpy): "

268 "grid \"%s\" is an empty grid", grid->name)

269

270 return

271 yac_field_data_set_add_mask_nocpy(

272 grid->field_data_set, location, mask, mask_name);

273}

275size_t yac_basic_grid_add_mask_nocpy_f2c(

276 struct yac_basic_grid * grid, int location,

277 int const * mask, char const * mask_name) {

278

279 return

280 yac_basic_grid_add_mask_nocpy(

281 grid, yac_get_location(location), mask, mask_name);

282}

284size_t yac_basic_grid_add_mask(

285 struct yac_basic_grid * grid, enum yac_location location,

286 int const * mask, size_t count, char const * mask_name) {

287

288 YAC_ASSERT(

289 grid, "ERROR(yac_basic_grid_add_mask): "

290 "NULL is not a valid value for argument grid")

291 YAC_ASSERT_F(

292 !grid->is_empty, "ERROR(yac_basic_grid_add_mask): "

293 "grid \"%s\" is an empty grid", grid->name)

294

295 return

296 yac_field_data_set_add_mask(

297 grid->field_data_set, location, mask, count, mask_name);

298}

300size_t yac_basic_grid_add_mask_f2c(

301 struct yac_basic_grid * grid, int location,

302 int const * mask, size_t count, char const * mask_name) {

303

304 return

305 yac_basic_grid_add_mask(

306 grid, yac_get_location(location), mask, count, mask_name);

307}

309struct yac_field_data * yac_basic_grid_get_field_data(

310 struct yac_basic_grid * grid, enum yac_location location) {

311

312 YAC_ASSERT(

313 grid, "ERROR(yac_basic_grid_get_field_data): "

314 "NULL is not a valid value for argument grid")

315

316 if (grid->is_empty)

317 return NULL;

318 else

319 return

320 yac_field_data_set_get_field_data(

321 grid->field_data_set, location);

322}

324struct yac_basic_grid * yac_basic_grid_reg_2d_new(

325 char const * name, size_t nbr_vertices[2], int cyclic[2],

326 double *lon_vertices, double *lat_vertices) {

327

328 return

329 yac_basic_grid_new(

330 name,

331 yac_generate_basic_grid_data_reg_2d(

332 nbr_vertices, cyclic, lon_vertices, lat_vertices));

333}

335struct yac_basic_grid * yac_basic_grid_reg_2d_deg_new(

336 char const * name, size_t nbr_vertices[2], int cyclic[2],

337 double *lon_vertices, double *lat_vertices) {

338

339 return

340 yac_basic_grid_new(

341 name,

342 yac_generate_basic_grid_data_reg_2d_deg(

343 nbr_vertices, cyclic, lon_vertices, lat_vertices));

344}

346struct yac_basic_grid * yac_basic_grid_curve_2d_new(

347 char const * name, size_t nbr_vertices[2], int cyclic[2],

348 double *lon_vertices, double *lat_vertices) {

349

350 return

351 yac_basic_grid_new(

352 name,

353 yac_generate_basic_grid_data_curve_2d(

354 nbr_vertices, cyclic, lon_vertices, lat_vertices));

355}

357struct yac_basic_grid * yac_basic_grid_curve_2d_deg_new(

358 char const * name, size_t nbr_vertices[2], int cyclic[2],

359 double *lon_vertices, double *lat_vertices) {

360

361 return

362 yac_basic_grid_new(

363 name,

364 yac_generate_basic_grid_data_curve_2d_deg(

365 nbr_vertices, cyclic, lon_vertices, lat_vertices));

366}

368struct yac_basic_grid * yac_basic_grid_unstruct_new(

369 char const * name, size_t nbr_vertices, size_t nbr_cells,

370 int *num_vertices_per_cell, double *x_vertices, double *y_vertices,

371 int *cell_to_vertex) {

372

373 return

374 yac_basic_grid_new(

375 name,

376 yac_generate_basic_grid_data_unstruct(

377 nbr_vertices, nbr_cells, num_vertices_per_cell,

378 x_vertices, y_vertices, cell_to_vertex));

379}

381struct yac_basic_grid * yac_basic_grid_unstruct_deg_new(

382 char const * name, size_t nbr_vertices, size_t nbr_cells,

383 int *num_vertices_per_cell, double *x_vertices, double *y_vertices,

384 int *cell_to_vertex) {

385

386 return

387 yac_basic_grid_new(

388 name,

389 yac_generate_basic_grid_data_unstruct_deg(

390 nbr_vertices, nbr_cells, num_vertices_per_cell,

391 x_vertices, y_vertices, cell_to_vertex));

392}

394struct yac_basic_grid * yac_basic_grid_unstruct_ll_new(

395 char const * name, size_t nbr_vertices, size_t nbr_cells,

396 int *num_vertices_per_cell, double *x_vertices, double *y_vertices,

397 int *cell_to_vertex) {

398

399 return

400 yac_basic_grid_new(

401 name,

402 yac_generate_basic_grid_data_unstruct_ll(

403 nbr_vertices, nbr_cells, num_vertices_per_cell,

404 x_vertices, y_vertices, cell_to_vertex));

405}

407struct yac_basic_grid * yac_basic_grid_unstruct_ll_deg_new(

408 char const * name, size_t nbr_vertices, size_t nbr_cells,

409 int *num_vertices_per_cell, double *x_vertices, double *y_vertices,

410 int *cell_to_vertex) {

411

412 return

413 yac_basic_grid_new(

414 name,

415 yac_generate_basic_grid_data_unstruct_ll_deg(

416 nbr_vertices, nbr_cells, num_vertices_per_cell,

417 x_vertices, y_vertices, cell_to_vertex));

418}

420struct yac_basic_grid * yac_basic_grid_unstruct_edge_new(

421 char const * name, size_t nbr_vertices, size_t nbr_cells, size_t nbr_edges,

422 int *num_edges_per_cell, double *x_vertices, double *y_vertices,

423 int *cell_to_edge, int *edge_to_vertex) {

424

425 return

426 yac_basic_grid_new(

427 name,

428 yac_generate_basic_grid_data_unstruct_edge(

429 nbr_vertices, nbr_cells, nbr_edges, num_edges_per_cell,

430 x_vertices, y_vertices, cell_to_edge, edge_to_vertex));

431}

433struct yac_basic_grid * yac_basic_grid_unstruct_edge_deg_new(

434 char const * name, size_t nbr_vertices, size_t nbr_cells, size_t nbr_edges,

435 int *num_edges_per_cell, double *x_vertices, double *y_vertices,

436 int *cell_to_edge, int *edge_to_vertex) {

437

438 return

439 yac_basic_grid_new(

440 name,

441 yac_generate_basic_grid_data_unstruct_edge_deg(

442 nbr_vertices, nbr_cells, nbr_edges, num_edges_per_cell,

443 x_vertices, y_vertices, cell_to_edge, edge_to_vertex));

444}

446struct yac_basic_grid * yac_basic_grid_unstruct_edge_ll_new(

447 char const * name, size_t nbr_vertices, size_t nbr_cells, size_t nbr_edges,

448 int *num_edges_per_cell, double *x_vertices, double *y_vertices,

449 int *cell_to_edge, int *edge_to_vertex) {

450

451 return

452 yac_basic_grid_new(

453 name,

454 yac_generate_basic_grid_data_unstruct_edge_ll(

455 nbr_vertices, nbr_cells, nbr_edges, num_edges_per_cell,

456 x_vertices, y_vertices, cell_to_edge, edge_to_vertex));

457}

459struct yac_basic_grid * yac_basic_grid_unstruct_edge_ll_deg_new(

460 char const * name, size_t nbr_vertices, size_t nbr_cells, size_t nbr_edges,

461 int *num_edges_per_cell, double *x_vertices, double *y_vertices,

462 int *cell_to_edge, int *edge_to_vertex) {

463

464 return

465 yac_basic_grid_new(

466 name,

467 yac_generate_basic_grid_data_unstruct_edge_ll_deg(

468 nbr_vertices, nbr_cells, nbr_edges, num_edges_per_cell,

469 x_vertices, y_vertices, cell_to_edge, edge_to_vertex));

470}

472struct yac_basic_grid * yac_basic_grid_cloud_new(

473 char const * name, size_t nbr_points, double *x_points, double *y_points) {

474

475 return

476 yac_basic_grid_new(

477 name,

478 yac_generate_basic_grid_data_cloud(nbr_points, x_points, y_points));

479}

481struct yac_basic_grid * yac_basic_grid_cloud_deg_new(

482 char const * name, size_t nbr_points, double *x_points, double *y_points) {

483

484 return

485 yac_basic_grid_new(

486 name,

487 yac_generate_basic_grid_data_cloud_deg(nbr_points, x_points, y_points));

488}

490struct yac_basic_grid * yac_basic_grid_reg_2d_rot_new(

491 char const * name, size_t nbr_vertices[2], int cyclic[2],

492 double *lon_vertices, double *lat_vertices,

493 double north_pole_lon, double north_pole_lat) {

494

495 return

496 yac_basic_grid_new(

497 name,

498 yac_generate_basic_grid_data_reg_2d_rot(

499 nbr_vertices, cyclic, lon_vertices, lat_vertices,

500 north_pole_lon, north_pole_lat));

501}

503struct yac_basic_grid * yac_basic_grid_reg_2d_rot_deg_new(

504 char const * name, size_t nbr_vertices[2], int cyclic[2],

505 double *lon_vertices, double *lat_vertices,

506 double north_pole_lon, double north_pole_lat) {

507

508 return

509 yac_basic_grid_new(

510 name,

511 yac_generate_basic_grid_data_reg_2d_rot_deg(

512 nbr_vertices, cyclic, lon_vertices, lat_vertices,

513 north_pole_lon, north_pole_lat));

514}

518static int def_dim(

519 char const * filename, int ncid, char const * name, size_t dim_len,

520 int file_is_new) {

521

522 int create_dim = file_is_new;

523 int dim_id;

524

525 // if the netcdf file already exists

526 if (!file_is_new) {

527

528 // check whether the dimension already exists in the file

529 int status = nc_inq_dimid(ncid, name, &dim_id);

530

531 // if the dimension already exists, check for a consistent dimension size

532 if (!((create_dim = (status == NC_EBADDIM)))) {

533

534 YAC_HANDLE_ERROR(status);

535

536 size_t temp_dim_len;

537 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dim_id, &temp_dim_len));

538

539 YAC_ASSERT_F(

540 dim_len == temp_dim_len,

541 "ERROR(def_dim): file \"%s\" already contains dimension \"%s\", "

542 "but it has a different size %zu != %zu",

543 filename, name, dim_len, temp_dim_len);

544 }

545 }

546

547 if (create_dim) YAC_HANDLE_ERROR(nc_def_dim(ncid, name, dim_len, &dim_id));

548

549 return dim_id;

550}

552static void def_var(

553 char const * filename, int ncid, char const * name, nc_type type,

554 int ndims, const int * dimids, char const * att_name, char const * att_value,

555 int file_is_new) {

556

557 int create_var = file_is_new;

558 int var_id;

559

560 // if the netcdf file already exists

561 if (!file_is_new) {

562

563 // check whether the variable already exists in the file

564 int status = nc_inq_varid(ncid, name, &var_id);

565

566 // if the variable already exists

567 if (!((create_var = (status == NC_ENOTVAR)))) {

568

569 YAC_HANDLE_ERROR(status);

570

571 // check the consistency of the variable definition

572 nc_type temp_type;

573 int temp_ndims;

574 int temp_dimids[NC_MAX_VAR_DIMS];

575 YAC_HANDLE_ERROR(

576 nc_inq_var(

577 ncid, var_id, NULL, &temp_type, &temp_ndims, temp_dimids, NULL));

578 YAC_ASSERT_F(

579 type == temp_type,

580 "ERROR(def_var): file \"%s\" already contains variable \"%s\", "

581 "but it has a different data type %d != %d",

582 filename, name, type, temp_type);

583 YAC_ASSERT_F(

584 ndims == temp_ndims,

585 "ERROR(def_var): file \"%s\" already contains variable \"%s\", "

586 "but it has a different number of dimensions %d != %d",

587 filename, name, ndims, temp_ndims);

588 for (int i = 0; i < ndims; ++i) {

589 YAC_ASSERT_F(

590 dimids[i] == temp_dimids[i],

591 "ERROR(def_var): file \"%s\" already contains variable \"%s\", "

592 "but it has a different dimensions dimid[%d] %d != %d",

593 filename, name, i, dimids[i], temp_dimids[i])

594 }

595

596 if (att_name != NULL) {

597

598 size_t att_len;

599 YAC_HANDLE_ERROR(nc_inq_attlen(ncid, var_id, att_name, &att_len));

600

601 char * temp_att_value = xcalloc((att_len + 1), sizeof(temp_att_value));

602 YAC_HANDLE_ERROR(nc_get_att_text(ncid, var_id, att_name, temp_att_value));

603

604 YAC_ASSERT_F(

605 !strcmp(att_value, temp_att_value),

606 "ERROR(def_var): file \"%s\" already contains variable \"%s\", "

607 "but it has a different attribute value (name: \"%s\") "

608 "\"%s\" != \"%s\"",

609 filename, name, att_name, att_value, temp_att_value);

610 free(temp_att_value);

611 }

612 }

613 }

614

615 if (create_var) {

616 YAC_HANDLE_ERROR(nc_def_var(ncid, name, type, ndims, dimids, &var_id));

617 if (att_name != NULL)

618 YAC_HANDLE_ERROR(

619 nc_put_att_text(ncid, var_id, att_name, strlen(att_value), att_value));

620 }

621}

623static void create_grid_file(

624 char const * filename, char const * grid_name, size_t num_cells,

625 int num_vertices_per_cell, int cell_center_coords_avaiable,

626 int cell_global_ids_available, int core_cell_mask_available,

627 int vertex_global_ids_available, int core_vertex_mask_available,

628 int edge_global_ids_available, int core_edge_mask_available) {

629

630 size_t grid_name_len = strlen(grid_name) + 1;

631

632 char nv_dim_name[3 + grid_name_len];

633 char nc_dim_name[3 + grid_name_len];

634

635 char cla_var_name[4 + grid_name_len];

636 char clo_var_name[4 + grid_name_len];

637 char lat_var_name[4 + grid_name_len];

638 char lon_var_name[4 + grid_name_len];

639 char rnk_var_name[4 + grid_name_len];

640 char cmk_var_name[4 + grid_name_len];

641 char gid_var_name[4 + grid_name_len];

642 char vcmk_var_name[5 + grid_name_len];

643 char vgid_var_name[5 + grid_name_len];

644 char ecmk_var_name[5 + grid_name_len];

645 char egid_var_name[5 + grid_name_len];

646 char const * unit_att = "units";

647 char const * coord_unit = "degree";

648

649 snprintf(nv_dim_name, 3 + grid_name_len, "nv_%s", grid_name);

650 snprintf(nc_dim_name, 3 + grid_name_len, "nc_%s", grid_name);

651

652 snprintf(cla_var_name, 4 + grid_name_len, "%s.cla", grid_name);

653 snprintf(clo_var_name, 4 + grid_name_len, "%s.clo", grid_name);

654 snprintf(lat_var_name, 4 + grid_name_len, "%s.lat", grid_name);

655 snprintf(lon_var_name, 4 + grid_name_len, "%s.lon", grid_name);

656 snprintf(rnk_var_name, 4 + grid_name_len, "%s.rnk", grid_name);

657 snprintf(cmk_var_name, 4 + grid_name_len, "%s.cmk", grid_name);

658 snprintf(gid_var_name, 4 + grid_name_len, "%s.gid", grid_name);

659 snprintf(vcmk_var_name, 5 + grid_name_len, "%s.vcmk", grid_name);

660 snprintf(vgid_var_name, 5 + grid_name_len, "%s.vgid", grid_name);

661 snprintf(ecmk_var_name, 5 + grid_name_len, "%s.ecmk", grid_name);

662 snprintf(egid_var_name, 5 + grid_name_len, "%s.egid", grid_name);

663

664 int ncid;

665

666 // create/open file

667 int file_is_new = !yac_file_exists(filename);

668 if (file_is_new) {

669 yac_nc_create(filename, NC_CLOBBER | NC_64BIT_OFFSET, &ncid);

670 } else {

671 yac_nc_open(filename, NC_WRITE | NC_SHARE, &ncid);

672 nc_redef(ncid);

673 }

674

675 // define dimensions

676 int nv_dim_id =

677 def_dim(filename, ncid, nv_dim_name, (size_t)num_vertices_per_cell, file_is_new);

678 int nc_dim_id =

679 def_dim(filename, ncid, nc_dim_name, (size_t)num_cells, file_is_new);

680

681 // define variables

682 int corner_dims[2] = {nc_dim_id, nv_dim_id};

683 int cell_dims[1] = {nc_dim_id};

684 def_var(

685 filename, ncid, cla_var_name, NC_DOUBLE, 2, corner_dims,

686 unit_att, coord_unit, file_is_new);

687 def_var(

688 filename, ncid, clo_var_name, NC_DOUBLE, 2, corner_dims,

689 unit_att, coord_unit, file_is_new);

690 if (cell_center_coords_avaiable) {

691 def_var(

692 filename, ncid, lat_var_name, NC_DOUBLE, 1, cell_dims,

693 unit_att, coord_unit, file_is_new);

694 def_var(

695 filename, ncid, lon_var_name, NC_DOUBLE, 1, cell_dims,

696 unit_att, coord_unit, file_is_new);

697 }

698 if (cell_global_ids_available)

699 def_var(

700 filename, ncid, gid_var_name, NC_INT, 1, cell_dims, NULL, NULL, file_is_new);

701 if (core_cell_mask_available)

702 def_var(

703 filename, ncid, cmk_var_name, NC_INT, 1, cell_dims, NULL, NULL, file_is_new);

704 if (vertex_global_ids_available)

705 def_var(

706 filename, ncid, vgid_var_name, NC_INT, 2, corner_dims, NULL, NULL, file_is_new);

707 if (core_vertex_mask_available)

708 def_var(

709 filename, ncid, vcmk_var_name, NC_INT, 2, corner_dims, NULL, NULL, file_is_new);

710 if (edge_global_ids_available)

711 def_var(

712 filename, ncid, egid_var_name, NC_INT, 2, corner_dims, NULL, NULL, file_is_new);

713 if (core_edge_mask_available)

714 def_var(

715 filename, ncid, ecmk_var_name, NC_INT, 2, corner_dims, NULL, NULL, file_is_new);

716 def_var(

717 filename, ncid, rnk_var_name, NC_INT, 1, cell_dims, NULL, NULL, file_is_new);

718

719 time_t now = time(NULL);

720 char str_now_UTC[32];

721 strftime(

722 str_now_UTC, sizeof(str_now_UTC), "%Y-%m-%dT%H:%M:%SZ",

723 gmtime(&now));

724 char const * yac_version = YAC_VERSION;

725 char const * created_by = "Created by YAC";

726 char const * grid_type = "curvilinear";

727 YAC_HANDLE_ERROR(

728 nc_put_att_text(ncid, NC_GLOBAL, "YAC", strlen(yac_version), yac_version));

729 YAC_HANDLE_ERROR(

730 nc_put_att_text(ncid, NC_GLOBAL, "title", strlen(created_by), created_by));

731 YAC_HANDLE_ERROR(

732 nc_put_att_text(ncid, NC_GLOBAL, "description", strlen(created_by), created_by));

733 YAC_HANDLE_ERROR(

734 nc_put_att_text(ncid, NC_GLOBAL, "grid", strlen(grid_type), grid_type));

735 YAC_HANDLE_ERROR(

736 nc_put_att_text(ncid, NC_GLOBAL, "timeStamp", strlen(str_now_UTC), str_now_UTC));

737

738 // end definition

739 YAC_HANDLE_ERROR(nc_enddef(ncid));

740

741 // close file

742 YAC_HANDLE_ERROR(nc_close(ncid));

743}

745static void put_vara(

746 int ncid, char const * grid_name, char const * var_ext,

747 size_t * start, size_t * count, void const * buffer) {

748

749 if (count[0] == 0) return;

750

751 char * var_name =

752 xmalloc((strlen(grid_name) + strlen(var_ext) + 2) * sizeof(*var_name));

753

754 sprintf(var_name, "%s.%s", grid_name, var_ext);

755

756 int var_id;

757

758 yac_nc_inq_varid(ncid, var_name, &var_id);

759

760 free(var_name);

761

762 YAC_HANDLE_ERROR(

763 nc_put_vara(ncid, var_id, start, count, buffer));

764}

768void yac_basic_grid_to_file_parallel(

769 struct yac_basic_grid * grid, char const * filename, MPI_Comm comm) {

770

771#ifndef YAC_NETCDF_ENABLED

772

773 UNUSED(grid);

774 UNUSED(filename);

775 UNUSED(comm);

776

777 die(

778 "ERROR(yac_basic_grid_to_file_parallel): "

779 "YAC is built without the NetCDF support");

780#else

781

782 int comm_rank, comm_size;

783 yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);

784 yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);

785

786 { // consistency check

787 int filename_len = (int)strlen(filename) + 1;

788 yac_mpi_call(MPI_Bcast(&filename_len, 1, MPI_INT, 0, comm), comm);

789 char * filename_buffer =

790 xmalloc((size_t)filename_len * sizeof(*filename_buffer));

791 if (comm_rank == 0) strcpy(filename_buffer, filename);

792 yac_mpi_call(

793 MPI_Bcast(filename_buffer, filename_len, MPI_CHAR, 0, comm), comm);

794

795 YAC_ASSERT_F(

796 !strcmp(filename, filename_buffer),

797 "ERROR(yac_basic_grid_to_file_parallel): "

798 "inconsistent filename (\"%s\" on rank 0 != \"%s\" on rank %d)",

799 filename_buffer, filename, comm_rank);

800 free(filename_buffer);

801 }

802

803 struct yac_field_data * cell_field =

804 yac_basic_grid_get_field_data(grid, YAC_LOC_CELL);

805 yac_const_coordinate_pointer cell_center_coords =

806 ((cell_field != NULL) &&

807 (yac_field_data_get_coordinates_count(cell_field) > 0))?

808 yac_field_data_get_coordinates_data(cell_field, 0):NULL;

809 struct yac_basic_grid_data * grid_data = yac_basic_grid_get_data(grid);

810

811 uint64_t local_cell_count = (uint64_t)grid->data.num_cells;

812 uint64_t * global_cell_counts =

813 xmalloc((size_t)comm_size * sizeof(*global_cell_counts));

814 int max_num_vertices_per_cell = 0;

815 int cell_center_coords_avaiable =

816 (grid_data->num_cells == 0) || (cell_center_coords != NULL);

817 int cell_global_ids_available =

818 (grid_data->num_cells == 0) || (grid_data->cell_ids != NULL);

819 int core_cell_mask_available =

820 (grid_data->num_cells == 0) || (grid_data->core_cell_mask != NULL);

821 int vertex_global_ids_available =

822 (grid_data->num_vertices == 0) || (grid_data->vertex_ids != NULL);

823 int core_vertex_mask_available =

824 (grid_data->num_vertices == 0) || (grid_data->core_vertex_mask != NULL);

825 int edge_global_ids_available =

826 (grid_data->num_edges == 0) || (grid_data->edge_ids != NULL);

827 int core_edge_mask_available =

828 (grid_data->num_edges == 0) || (grid_data->core_edge_mask != NULL);

829

830 for (size_t i = 0; i < grid->data.num_cells; ++i)

831 if (grid->data.num_vertices_per_cell[i] > max_num_vertices_per_cell)

832 max_num_vertices_per_cell = grid->data.num_vertices_per_cell[i];

833

834 yac_mpi_call(

835 MPI_Allgather(

836 &local_cell_count, 1, MPI_UINT64_T,

837 global_cell_counts, 1, MPI_UINT64_T, comm), comm);

838 yac_mpi_call(

839 MPI_Allreduce(

840 MPI_IN_PLACE, &max_num_vertices_per_cell, 1, MPI_INT, MPI_MAX, comm),

841 comm);

842 int ints[7] = {cell_center_coords_avaiable,

843 cell_global_ids_available,

844 core_cell_mask_available,

845 vertex_global_ids_available,

846 core_vertex_mask_available,

847 edge_global_ids_available,

848 core_edge_mask_available};

849 yac_mpi_call(

850 MPI_Allreduce(

851 MPI_IN_PLACE, ints, sizeof(ints)/sizeof(ints[0]),

852 MPI_INT, MPI_MIN, comm), comm);

853 cell_center_coords_avaiable = ints[0];

854 cell_global_ids_available = ints[1];

855 core_cell_mask_available = ints[2];

856 vertex_global_ids_available = ints[3];

857 core_vertex_mask_available = ints[4];

858 edge_global_ids_available = ints[5];

859 core_edge_mask_available = ints[6];

860

861 size_t global_cell_count = 0;

862 for (int i = 0; i < comm_size; ++i)

863 global_cell_count += (size_t)(global_cell_counts[i]);

864

865 YAC_ASSERT_F(

866 global_cell_count > 0,

867 "ERROR(yac_basic_grid_to_file_parallel): grid \"%s\" has no cells",

868 grid->name);

869

870 // create the grid file

871 if (comm_rank == 0)

872 create_grid_file(

873 filename, grid->name, global_cell_count, max_num_vertices_per_cell,

874 cell_center_coords_avaiable,

875 cell_global_ids_available, core_cell_mask_available,

876 vertex_global_ids_available, core_vertex_mask_available,

877 edge_global_ids_available, core_edge_mask_available);

878 yac_mpi_call(MPI_Barrier(comm), comm);

879

880 // determine processes that will do output

881 int io_flag;

882 int * io_ranks;

883 int num_io_ranks;

884 yac_get_io_ranks(comm, &io_flag, &io_ranks, &num_io_ranks);

885

886 size_t recv_count = 0;

887 size_t io_start_idx = SIZE_MAX;

888

889 if (io_flag) {

890

891 int io_rank_idx;

892 for (io_rank_idx = 0; io_rank_idx < num_io_ranks; ++io_rank_idx)

893 if (io_ranks[io_rank_idx] == comm_rank) break;

894

895 io_start_idx =

896 (size_t)(

897 ((long long)(io_rank_idx) * (long long)global_cell_count)/

898 (long long)num_io_ranks);

899 recv_count =

900 (size_t)(

901 ((long long)(io_rank_idx + 1) * (long long)global_cell_count)/

902 (long long)num_io_ranks - (long long)io_start_idx);

903 }

904

905 // generate transfer positions

906 struct Xt_com_pos * com_pos_buffer =

907 xmalloc(5 * (size_t)comm_size * sizeof(*com_pos_buffer));

908 struct Xt_com_pos * cell_src_com = com_pos_buffer;

909 struct Xt_com_pos * cell_dst_com = com_pos_buffer + (size_t)comm_size;

910 struct Xt_com_pos * vertex_src_com = com_pos_buffer + 2 * (size_t)comm_size;

911 struct Xt_com_pos * vertex_dst_com = com_pos_buffer + 3 * (size_t)comm_size;

912 struct Xt_com_pos * edge_src_com = com_pos_buffer + 4 * (size_t)comm_size;

913 struct Xt_com_pos * edge_dst_com = vertex_dst_com;

914 int num_src_msg = 0, num_dst_msg = 0;

915 int * transfer_pos_buffer =

916 xmalloc(

917 (size_t)(2 * max_num_vertices_per_cell + 1) *

918 (grid->data.num_cells + recv_count) * sizeof(*transfer_pos_buffer));

919 int * cell_send_pos = transfer_pos_buffer;

920 int * cell_recv_pos = transfer_pos_buffer + grid->data.num_cells;

921 int * vertex_send_pos =

922 transfer_pos_buffer + grid->data.num_cells + recv_count;

923 int * vertex_recv_pos =

924 transfer_pos_buffer + grid->data.num_cells + recv_count +

925 (size_t)max_num_vertices_per_cell * grid->data.num_cells;

926 int * edge_send_pos =

927 transfer_pos_buffer + grid->data.num_cells + recv_count +

928 (size_t)max_num_vertices_per_cell * grid->data.num_cells +

929 (size_t)max_num_vertices_per_cell * recv_count;

930

931 for (size_t i = 0; i < grid->data.num_cells; ++i) {

932

933 size_t * cell_to_vertex =

934 grid->data.cell_to_vertex + grid->data.cell_to_vertex_offsets[i];

935 size_t * cell_to_edge =

936 grid->data.cell_to_edge + grid->data.cell_to_edge_offsets[i];

937 int num_vertices = grid->data.num_vertices_per_cell[i];

938

939 cell_send_pos[i] = (int)i;

940 for (int j = 0; j < num_vertices; ++j)

941 vertex_send_pos[i * (size_t)max_num_vertices_per_cell + (size_t)j] =

942 cell_to_vertex[j];

943 for (int j = num_vertices; j < max_num_vertices_per_cell; ++j)

944 vertex_send_pos[i * (size_t)max_num_vertices_per_cell + (size_t)j] = 0;

945 for (int j = 0; j < num_vertices; ++j)

946 edge_send_pos[i * (size_t)max_num_vertices_per_cell + (size_t)j] =

947 cell_to_edge[j];

948 for (int j = num_vertices; j < max_num_vertices_per_cell; ++j)

949 edge_send_pos[i * (size_t)max_num_vertices_per_cell + (size_t)j] = 0;

950 }

951

952 for (size_t i = 0, k = 0; i < recv_count; ++i) {

953 cell_recv_pos[i] = (int)i;

954 for (int j = 0; j < max_num_vertices_per_cell; ++j, ++k)

955 vertex_recv_pos[k] = k;

956 }

957

958 size_t io_end_idx = 0;

959 size_t global_count_accu = 0;

960 for (int io_rank_idx = 0, curr_rank = 0; io_rank_idx < num_io_ranks;

961 ++io_rank_idx) {

962

963 int curr_io_rank = io_ranks[io_rank_idx];

964 size_t io_start_idx = io_end_idx;

965 io_end_idx =

966 (size_t)(

967 ((long long)(io_rank_idx + 1) * (long long)global_cell_count)/

968 (long long)num_io_ranks);

969 size_t io_size = io_end_idx - io_start_idx;

970 size_t curr_recv_size = 0;

971

972 while(global_count_accu < io_end_idx) {

973

974 size_t count =

975 MIN(

976 (global_count_accu + global_cell_counts[curr_rank]) -

977 (io_start_idx + curr_recv_size), io_size - curr_recv_size);

978

979 if (curr_rank == comm_rank) {

980 cell_src_com[num_src_msg].transfer_pos = cell_send_pos;

981 cell_src_com[num_src_msg].num_transfer_pos = (int)count;

982 cell_src_com[num_src_msg].rank = curr_io_rank;

983 cell_send_pos += count;

984

985 vertex_src_com[num_src_msg].transfer_pos = vertex_send_pos;

986 vertex_src_com[num_src_msg].num_transfer_pos =

987 (int)count * max_num_vertices_per_cell;

988 vertex_src_com[num_src_msg].rank = curr_io_rank;

989 vertex_send_pos +=

990 count * (size_t)max_num_vertices_per_cell;

991

992 edge_src_com[num_src_msg].transfer_pos = edge_send_pos;

993 edge_src_com[num_src_msg].num_transfer_pos =

994 (int)count * max_num_vertices_per_cell;

995 edge_src_com[num_src_msg].rank = curr_io_rank;

996 edge_send_pos +=

997 count * (size_t)max_num_vertices_per_cell;

998

999 num_src_msg++;

1000 }

1001

1002 if (curr_io_rank == comm_rank) {

1003 cell_dst_com[num_dst_msg].transfer_pos = cell_recv_pos;

1004 cell_dst_com[num_dst_msg].num_transfer_pos = (int)count;

1005 cell_dst_com[num_dst_msg].rank = curr_rank;

1006 cell_recv_pos += count;

1007

1008 vertex_dst_com[num_dst_msg].transfer_pos = vertex_recv_pos;

1009 vertex_dst_com[num_dst_msg].num_transfer_pos =

1010 (int)count * max_num_vertices_per_cell;

1011 vertex_dst_com[num_dst_msg].rank = curr_rank;

1012 vertex_recv_pos += count * (size_t)max_num_vertices_per_cell;

1013

1014 num_dst_msg++;

1015 }

1016

1017 if ((global_count_accu + global_cell_counts[curr_rank]) <=

1018 io_end_idx) {

1019 global_count_accu += global_cell_counts[curr_rank];

1020 curr_rank++;

1021 }

1022

1023 curr_recv_size += count;

1024

1025 if (curr_recv_size >= io_size) break;

1026 }

1027 }

1028

1029 free(global_cell_counts);

1030

1031 // open grid file

1032 int ncid;

1033 size_t start[2], count[2];

1034 if (io_flag) {

1035

1036 int io_rank_idx;

1037 for (io_rank_idx = 0; io_rank_idx < num_io_ranks; ++io_rank_idx)

1038 if (io_ranks[io_rank_idx] == comm_rank) break;

1039

1040 start[0] = io_start_idx;

1041 start[1] = 0;

1042 count[0] = recv_count;

1043 count[1] = max_num_vertices_per_cell;

1044

1045 yac_nc_open(filename, NC_WRITE | NC_SHARE, &ncid);

1046 } else {

1047 count[0] = 0;

1048 count[1] = 0;

1049 }

1050

1051 void * recv_buffer =

1052 xmalloc(

1053 (size_t)max_num_vertices_per_cell * recv_count *

1054 MAX(MAX(sizeof(double), sizeof(int)), sizeof(yac_int)));

1055

1056 //----------------------------------------------------------------------------

1057 // redistribute cell based data

1058 //----------------------------------------------------------------------------

1059

1060 Xt_xmap cell_xmap =

1061 xt_xmap_intersection_pos_new(

1062 num_src_msg, cell_src_com, num_dst_msg, cell_dst_com, comm);

1063

1064 Xt_redist cell_redist_int = xt_redist_p2p_new(cell_xmap, MPI_INT);

1065

1066 // number of vertices per cell

1067 int * num_vertices_per_cell =

1068 xmalloc(recv_count * sizeof(*num_vertices_per_cell));

1069 xt_redist_s_exchange1(

1070 cell_redist_int, grid->data.num_vertices_per_cell, num_vertices_per_cell);

1071

1072 // cell center coordinates (if available)

1073 if (cell_center_coords_avaiable) {

1074

1075 // generate cell center lon/lat data

1076 double * send_buffer_dble =

1077 xmalloc(2 * grid->data.num_cells * sizeof(*send_buffer_dble));;

1078 double * send_buffer_lon = send_buffer_dble;

1079 double * send_buffer_lat = send_buffer_dble + grid->data.num_cells;

1080 for (size_t i = 0; i < grid->data.num_cells; ++i) {

1081 XYZtoLL(

1082 cell_center_coords[i], send_buffer_lon + i, send_buffer_lat + i);

1083 send_buffer_lon[i] /= YAC_RAD;

1084 send_buffer_lat[i] /= YAC_RAD;

1085 }

1086

1087 // exchange cell center lon/lat data and write it to file

1088 Xt_redist cell_redist_dble = xt_redist_p2p_new(cell_xmap, MPI_DOUBLE);

1089 xt_redist_s_exchange1(cell_redist_dble, send_buffer_lon, recv_buffer);

1090 put_vara(ncid, grid->name, "lon", start, count, recv_buffer);

1091 xt_redist_s_exchange1(cell_redist_dble, send_buffer_lat, recv_buffer);

1092 put_vara(ncid, grid->name, "lat", start, count, recv_buffer);

1093 xt_redist_delete(cell_redist_dble);

1094 free(send_buffer_dble);

1095 }

1096

1097 int * cell_send_buffer_int =

1098 xmalloc(grid->data.num_cells * sizeof(*cell_send_buffer_int));

1099

1100 // cell ranks

1101 {

1102 // generate cell owner rank data

1103 for (size_t i = 0; i < grid->data.num_cells; ++i)

1104 cell_send_buffer_int[i] = comm_rank;

1105

1106 // exchange cell owner ranks and write it to file

1107 xt_redist_s_exchange1(cell_redist_int, cell_send_buffer_int, recv_buffer);

1108 put_vara(ncid, grid->name, "rnk", start, count, recv_buffer);

1109 }

1110

1111 // cell core mask (if available)

1112 if (core_cell_mask_available) {

1113

1114 // exchange core cell mask and write to file

1115 xt_redist_s_exchange1(

1116 cell_redist_int, grid->data.core_cell_mask, recv_buffer);

1117 put_vara(ncid, grid->name, "cmk", start, count, recv_buffer);

1118 }

1119

1120 // cell global ids (if available)

1121 if (cell_global_ids_available) {

1122

1123 // convert global ids to int

1124 for (size_t i = 0; i < grid->data.num_cells; ++i)

1125 cell_send_buffer_int[i] = (int)(grid->data.cell_ids[i]);

1126

1127 // exchange cell global ids and write to file

1128 xt_redist_s_exchange1(

1129 cell_redist_int, cell_send_buffer_int, recv_buffer);

1130 put_vara(ncid, grid->name, "gid", start, count, recv_buffer);

1131 }

1132

1133 xt_redist_delete(cell_redist_int);

1134 free(cell_send_buffer_int);

1135

1136 xt_xmap_delete(cell_xmap);

1137

1138 //----------------------------------------------------------------------------

1139 // redistribute vertex based data

1140 //----------------------------------------------------------------------------

1141

1142 Xt_xmap vertex_xmap =

1143 xt_xmap_intersection_pos_new(

1144 num_src_msg, vertex_src_com, num_dst_msg, vertex_dst_com, comm);

1145

1146 // vertex coordinates

1147 {

1148 // generate vertex lon/lat data

1149 double * send_buffer_dble =

1150 xmalloc(2 * grid->data.num_vertices * sizeof(*send_buffer_dble));

1151 double * send_buffer_lon = send_buffer_dble;

1152 double * send_buffer_lat = send_buffer_dble + grid->data.num_vertices;

1153 for (size_t i = 0; i < grid->data.num_vertices; ++i) {

1154 XYZtoLL(

1155 grid->data.vertex_coordinates[i],

1156 send_buffer_lon + i, send_buffer_lat + i);

1157 send_buffer_lon[i] /= YAC_RAD;

1158 send_buffer_lat[i] /= YAC_RAD;

1159 }

1160

1161 // exchange vertex lon/lat data and write it to file

1162 Xt_redist vertex_redist_dble = xt_redist_p2p_new(vertex_xmap, MPI_DOUBLE);

1163 xt_redist_s_exchange1(vertex_redist_dble, send_buffer_lon, recv_buffer);

1164 for (size_t i = 0, k = 0; i < recv_count;

1165 ++i, k += (size_t)max_num_vertices_per_cell)

1166 for (size_t j = (size_t)num_vertices_per_cell[i];

1167 j < (size_t)max_num_vertices_per_cell; ++j)

1168 ((double*)recv_buffer)[k+j] = DBL_MAX;

1169 put_vara(ncid, grid->name, "clo", start, count, recv_buffer);

1170 xt_redist_s_exchange1(vertex_redist_dble, send_buffer_lat, recv_buffer);

1171 free(send_buffer_dble);

1172 for (size_t i = 0, k = 0; i < recv_count;

1173 ++i, k += (size_t)max_num_vertices_per_cell)

1174 for (size_t j = (size_t)num_vertices_per_cell[i];

1175 j < (size_t)max_num_vertices_per_cell; ++j)

1176 ((double*)recv_buffer)[k+j] = DBL_MAX;

1177 put_vara(ncid, grid->name, "cla", start, count, recv_buffer);

1178 xt_redist_delete(vertex_redist_dble);

1179 }

1180

1181 Xt_redist vertex_redist_int = xt_redist_p2p_new(vertex_xmap, MPI_INT);

1182 int * vertex_send_buffer_int =

1183 xmalloc(grid->data.num_vertices * sizeof(*vertex_send_buffer_int));

1184

1185 // core vertex mask (if available)

1186 if (core_vertex_mask_available) {

1187

1188 // exchange core vertex mask and write to file

1189 xt_redist_s_exchange1(

1190 vertex_redist_int, grid->data.core_vertex_mask, recv_buffer);

1191 for (size_t i = 0, k = 0; i < recv_count;

1192 ++i, k += (size_t)max_num_vertices_per_cell)

1193 for (size_t j = (size_t)num_vertices_per_cell[i];

1194 j < (size_t)max_num_vertices_per_cell; ++j)

1195 ((int*)recv_buffer)[k+j] = INT_MAX;

1196 put_vara(ncid, grid->name, "vcmk", start, count, recv_buffer);

1197 }

1198

1199 // vertex global ids (if available)

1200 if (vertex_global_ids_available) {

1201

1202 // convert global ids to int

1203 for (size_t i = 0; i < grid->data.num_vertices; ++i)

1204 vertex_send_buffer_int[i] = (int)(grid->data.vertex_ids[i]);

1205

1206 // exchange vertex global ids and write to file

1207 xt_redist_s_exchange1(

1208 vertex_redist_int, vertex_send_buffer_int, recv_buffer);

1209 for (size_t i = 0, k = 0; i < recv_count;

1210 ++i, k += (size_t)max_num_vertices_per_cell)

1211 for (size_t j = (size_t)num_vertices_per_cell[i];

1212 j < (size_t)max_num_vertices_per_cell; ++j)

1213 ((int*)recv_buffer)[k+j] = INT_MAX;

1214 put_vara(ncid, grid->name, "vgid", start, count, recv_buffer);

1215 }

1216

1217 free(vertex_send_buffer_int);

1218 xt_redist_delete(vertex_redist_int);

1219

1220 xt_xmap_delete(vertex_xmap);

1221

1222 //----------------------------------------------------------------------------

1223 // redistribute edge based data

1224 //----------------------------------------------------------------------------

1225

1226 if (core_edge_mask_available || edge_global_ids_available) {

1227

1228 Xt_xmap edge_xmap =

1229 xt_xmap_intersection_pos_new(

1230 num_src_msg, edge_src_com, num_dst_msg, edge_dst_com, comm);

1231

1232 Xt_redist edge_redist_int = xt_redist_p2p_new(edge_xmap, MPI_INT);

1233 int * edge_send_buffer_int =

1234 xmalloc(grid->data.num_edges * sizeof(*edge_send_buffer_int));

1235

1236 // core edge mask (if available)

1237 if (core_edge_mask_available) {

1238

1239 // exchange core edge mask and write to file

1240 xt_redist_s_exchange1(

1241 edge_redist_int, grid->data.core_edge_mask, recv_buffer);

1242 for (size_t i = 0, k = 0; i < recv_count;

1243 ++i, k += (size_t)max_num_vertices_per_cell)

1244 for (size_t j = (size_t)num_vertices_per_cell[i];

1245 j < (size_t)max_num_vertices_per_cell; ++j)

1246 ((int*)recv_buffer)[k+j] = INT_MAX;

1247 put_vara(ncid, grid->name, "ecmk", start, count, recv_buffer);

1248 }

1249

1250 // edge global ids (if available)

1251 if (edge_global_ids_available) {

1252

1253 // convert global ids to int

1254 for (size_t i = 0; i < grid->data.num_edges; ++i)

1255 edge_send_buffer_int[i] = (int)(grid->data.edge_ids[i]);

1256

1257 // exchange edge global ids and write to file

1258 xt_redist_s_exchange1(

1259 edge_redist_int, edge_send_buffer_int, recv_buffer);

1260 for (size_t i = 0, k = 0; i < recv_count;

1261 ++i, k += (size_t)max_num_vertices_per_cell)

1262 for (size_t j = (size_t)num_vertices_per_cell[i];

1263 j < (size_t)max_num_vertices_per_cell; ++j)

1264 ((int*)recv_buffer)[k+j] = INT_MAX;

1265 put_vara(ncid, grid->name, "egid", start, count, recv_buffer);

1266 }

1267

1268 free(edge_send_buffer_int);

1269 xt_redist_delete(edge_redist_int);

1270

1271 xt_xmap_delete(edge_xmap);

1272 }

1273

1274 // close file

1275 if (io_flag) YAC_HANDLE_ERROR(nc_close(ncid));

1276

1277 free(num_vertices_per_cell);

1278 free(recv_buffer);

1279 free(com_pos_buffer);

1280 free(transfer_pos_buffer);

1281 free(io_ranks);

1282

1283 // ensure that the writing of the weight file is complete

1284 yac_mpi_call(MPI_Barrier(comm), comm);

1285

1286#endif // YAC_NETCDF_ENABLED

1287}

1289void yac_basic_grid_compute_cell_areas(

1290 struct yac_basic_grid * grid, double * cell_areas) {

1291

1292 yac_basic_grid_data_compute_cell_areas(

1293 grid->data, cell_areas);

1294}