1 // CGNS Assumptions: 2 // * All boundary conditions are listed as Family nodes at the "top" level. 3 // * Single Base. 4 // * ZoneGridConnectivity is 1to1 with point lists for unstructured 5 6 // Copyright(C) 1999-2020 National Technology & Engineering Solutions 7 // of Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with 8 // NTESS, the U.S. Government retains certain rights in this software. 9 // 10 // See packages/seacas/LICENSE for details 11 12 #include <cgns/Iocgns_Defines.h> 13 14 #include <Ioss_CodeTypes.h> 15 #include <Ioss_Sort.h> 16 #include <Ioss_Utils.h> 17 #include <cassert> 18 #include <cgns/Iocgns_ParallelDatabaseIO.h> 19 #include <cgns/Iocgns_Utils.h> 20 #include <cstddef> 21 #include <ctime> 22 #include <fmt/ostream.h> 23 #include <fstream> 24 #include <iostream> 25 #include <mpi.h> 26 #include <numeric> 27 #include <string> 28 #include <sys/select.h> 29 #include <vector> 30 31 #include <vtk_cgns.h> // xxx(kitware) 32 #ifdef SEACAS_HAVE_MPI 33 #include VTK_CGNS(pcgnslib.h) 34 #else 35 #include VTK_CGNS(cgnslib.h) 36 #endif 37 38 #if !defined(CGNSLIB_H) 39 #error "Could not include cgnslib.h" 40 #endif 41 42 #include "Ioss_CommSet.h" 43 #include "Ioss_DBUsage.h" 44 #include "Ioss_DatabaseIO.h" 45 #include "Ioss_ElementBlock.h" 46 #include "Ioss_ElementTopology.h" 47 #include "Ioss_EntityType.h" 48 #include "Ioss_Field.h" 49 #include "Ioss_FileInfo.h" 50 #include "Ioss_IOFactory.h" 51 #include "Ioss_NodeBlock.h" 52 #include "Ioss_ParallelUtils.h" 53 #include "Ioss_Property.h" 54 #include "Ioss_Region.h" 55 #include "Ioss_SideBlock.h" 56 #include "Ioss_SideSet.h" 57 #include "Ioss_State.h" 58 #include "Ioss_StructuredBlock.h" 59 #include "Ioss_Utils.h" 60 #include "Ioss_VariableType.h" 61 62 using GL_IdVector = std::vector<std::pair<int, int>>; 63 64 namespace { cgns_mpi_type()65 MPI_Datatype cgns_mpi_type() 66 { 67 #if CG_SIZEOF_SIZE == 8 68 return MPI_LONG_LONG_INT; 69 #else 70 return MPI_INT; 71 #endif 72 } 73 map_local_to_global_implicit(CGNSIntVector & data,ssize_t count,const CGNSIntVector & global_implicit_map)74 void map_local_to_global_implicit(CGNSIntVector &data, ssize_t count, 75 const CGNSIntVector &global_implicit_map) 76 { 77 for (ssize_t i = 0; i < count; i++) { 78 data[i] = global_implicit_map[data[i] - 1]; 79 } 80 } 81 82 GL_IdVector gather_nodes_to_proc0(Ioss::Map &global_id_map, int processor, int64_t offset, 83 const Ioss::ParallelUtils &util, size_t min_id, 84 size_t max_id = std::numeric_limits<size_t>::max()) 85 { 86 GL_IdVector I_nodes; 87 GL_IdVector I_nodes_recv; 88 for (size_t i = 0; i < global_id_map.size(); i++) { 89 auto global_id = global_id_map.map()[i + 1]; 90 if ((size_t)global_id >= min_id && (size_t)global_id <= max_id) { 91 I_nodes.emplace_back((int)global_id, (int)i + 1 + offset); 92 } 93 } 94 95 std::vector<int> recv_count; 96 util.gather(2 * (int)I_nodes.size(), recv_count); 97 std::vector<int> recv_off(recv_count); 98 99 auto count = std::accumulate(recv_count.begin(), recv_count.end(), 0); 100 if (processor == 0) { 101 I_nodes_recv.resize(count / 2); 102 Ioss::Utils::generate_index(recv_off); 103 } 104 105 MPI_Gatherv(I_nodes.data(), (int)I_nodes.size() * 2, MPI_INT, I_nodes_recv.data(), 106 recv_count.data(), recv_off.data(), MPI_INT, 0, util.communicator()); 107 108 if (processor == 0) { 109 Ioss::sort(I_nodes_recv.begin(), I_nodes_recv.end()); 110 } 111 return I_nodes_recv; 112 } 113 intersect(const GL_IdVector & I,const GL_IdVector & J)114 GL_IdVector intersect(const GL_IdVector &I, const GL_IdVector &J) 115 { 116 // Find all common global_ids (entry.first) between 'I' and 'J'. 117 // When found, store the proc-zone-local position (entry.second) 118 // in 'common' 119 // PRECONDITION: 'I' and 'J' are sorted. 120 121 auto min_size = std::min(I.size(), J.size()); 122 GL_IdVector common; 123 common.reserve(min_size); 124 125 auto II = I.begin(); 126 auto JJ = J.begin(); 127 while (II != I.end() && JJ != J.end()) { 128 if ((*II).first < (*JJ).first) { 129 ++II; 130 } 131 else { 132 if (!((*JJ).first < (*II).first)) { 133 common.emplace_back((*II).second, (*JJ).second); 134 ++II; 135 } 136 ++JJ; 137 } 138 } 139 common.shrink_to_fit(); 140 return common; 141 } 142 } // namespace 143 144 namespace Iocgns { 145 ParallelDatabaseIO(Ioss::Region * region,const std::string & filename,Ioss::DatabaseUsage db_usage,MPI_Comm communicator,const Ioss::PropertyManager & props)146 ParallelDatabaseIO::ParallelDatabaseIO(Ioss::Region *region, const std::string &filename, 147 Ioss::DatabaseUsage db_usage, MPI_Comm communicator, 148 const Ioss::PropertyManager &props) 149 : Ioss::DatabaseIO(region, filename, db_usage, communicator, props) 150 { 151 usingParallelIO = true; 152 dbState = Ioss::STATE_UNKNOWN; 153 154 #if IOSS_DEBUG_OUTPUT 155 if (myProcessor == 0) { 156 fmt::print("CGNS ParallelDatabaseIO using {}-bit integers.\n" 157 " using the parallel CGNS library and API.\n", 158 CG_SIZEOF_SIZE); 159 } 160 #endif 161 if (!is_input()) { 162 if (properties.exists("FLUSH_INTERVAL")) { 163 m_flushInterval = properties.get("FLUSH_INTERVAL").get_int(); 164 } 165 166 { 167 bool file_per_state = false; 168 Ioss::Utils::check_set_bool_property(properties, "FILE_PER_STATE", file_per_state); 169 if (file_per_state) { 170 set_file_per_state(true); 171 } 172 } 173 } 174 175 Ioss::DatabaseIO::openDatabase__(); 176 } 177 ~ParallelDatabaseIO()178 ParallelDatabaseIO::~ParallelDatabaseIO() 179 { 180 for (auto >b : m_globalToBlockLocalNodeMap) { 181 delete gtb.second; 182 } 183 try { 184 closeBaseDatabase__(); 185 closeDatabase__(); 186 } 187 catch (...) { 188 } 189 } 190 get_file_pointer()191 int ParallelDatabaseIO::get_file_pointer() const 192 { 193 if (m_cgnsFilePtr < 0) { 194 openDatabase__(); 195 } 196 return m_cgnsFilePtr; 197 } 198 openDatabase__()199 void ParallelDatabaseIO::openDatabase__() const 200 { 201 if (m_cgnsFilePtr < 0) { 202 int mode = is_input() ? CG_MODE_READ : CG_MODE_WRITE; 203 if (!is_input()) { 204 if (m_cgnsFilePtr == -2) { 205 // Writing multiple steps with a "flush" (cg_close() / cg_open()) 206 mode = CG_MODE_MODIFY; 207 } 208 else { 209 // Check whether appending to existing file... 210 if (open_create_behavior() == Ioss::DB_APPEND || 211 open_create_behavior() == Ioss::DB_MODIFY) { 212 // Append to file if it already exists -- See if the file exists. 213 Ioss::FileInfo file = Ioss::FileInfo(decoded_filename()); 214 if (file.exists()) { 215 mode = CG_MODE_MODIFY; 216 } 217 } 218 } 219 } 220 221 bool do_timer = false; 222 Ioss::Utils::check_set_bool_property(properties, "IOSS_TIME_FILE_OPEN_CLOSE", do_timer); 223 double t_begin = (do_timer ? Ioss::Utils::timer() : 0); 224 225 CGCHECKM(cg_set_file_type(CG_FILE_HDF5)); 226 227 #if CGNS_VERSION > 3320 228 CGCHECKM(cgp_mpi_comm(util().communicator())); 229 #else 230 // Older versions of cgp_mpi_comm returned an internal NO_ERROR 231 // value which is equal to -1. 232 cgp_mpi_comm(util().communicator()); 233 #endif 234 CGCHECKM(cgp_pio_mode(CGP_COLLECTIVE)); 235 Ioss::DatabaseIO::openDatabase__(); 236 int ierr = cgp_open(get_dwname().c_str(), mode, &m_cgnsFilePtr); 237 238 if (do_timer) { 239 double t_end = Ioss::Utils::timer(); 240 double duration = util().global_minmax(t_end - t_begin, Ioss::ParallelUtils::DO_MAX); 241 if (myProcessor == 0) { 242 fmt::print(Ioss::DEBUG(), "{} File Open Time = {}\n", is_input() ? "Input" : "Output", 243 duration); 244 } 245 } 246 247 if (ierr != CG_OK) { 248 // NOTE: Code will not continue past this call... 249 std::ostringstream errmsg; 250 fmt::print(errmsg, "ERROR: Problem opening file '{}' for {} access. CGNS Error: '{}'", 251 get_filename(), (is_input() ? "read" : "write"), cg_get_error()); 252 IOSS_ERROR(errmsg); 253 } 254 255 if (properties.exists("INTEGER_SIZE_API")) { 256 int isize = properties.get("INTEGER_SIZE_API").get_int(); 257 if (isize == 8) { 258 set_int_byte_size_api(Ioss::USE_INT64_API); 259 } 260 if (isize == 4) { 261 set_int_byte_size_api(Ioss::USE_INT32_API); 262 } 263 } 264 else if (CG_SIZEOF_SIZE == 64) { 265 set_int_byte_size_api(Ioss::USE_INT64_API); 266 } 267 268 if (mode == CG_MODE_MODIFY && get_region() != nullptr) { 269 Utils::update_db_zone_property(m_cgnsFilePtr, get_region(), myProcessor, true, true); 270 } 271 #if 0 272 // This isn't currently working since CGNS currently has chunking 273 // disabled for HDF5 files and compression requires chunking. 274 if (!is_input()) { 275 if (properties.exists("COMPRESSION_LEVEL")) { 276 int comp = properties.get("COMPRESSION_LEVEL").get_int(); 277 cg_configure(CG_CONFIG_HDF5_COMPRESS, (void*)comp); 278 } 279 } 280 #endif 281 } 282 assert(m_cgnsFilePtr >= 0); 283 } 284 closeBaseDatabase__()285 void ParallelDatabaseIO::closeBaseDatabase__() const 286 { 287 if (m_cgnsBasePtr > 0) { 288 bool do_timer = false; 289 Ioss::Utils::check_set_bool_property(properties, "IOSS_TIME_FILE_OPEN_CLOSE", do_timer); 290 double t_begin = (do_timer ? Ioss::Utils::timer() : 0); 291 292 CGCHECKM(cg_close(m_cgnsBasePtr)); 293 m_cgnsBasePtr = -1; 294 295 if (do_timer) { 296 double t_end = Ioss::Utils::timer(); 297 double duration = util().global_minmax(t_end - t_begin, Ioss::ParallelUtils::DO_MAX); 298 if (myProcessor == 0) { 299 fmt::print(Ioss::DEBUG(), "{} Base File Close Time = {}\n", 300 is_input() ? "Input" : "Output", duration); 301 } 302 } 303 } 304 } 305 closeDatabase__()306 void ParallelDatabaseIO::closeDatabase__() const 307 { 308 if (m_cgnsFilePtr > 0) { 309 bool do_timer = false; 310 Ioss::Utils::check_set_bool_property(properties, "IOSS_TIME_FILE_OPEN_CLOSE", do_timer); 311 double t_begin = (do_timer ? Ioss::Utils::timer() : 0); 312 313 CGCHECKM(cgp_close(m_cgnsFilePtr)); 314 315 if (do_timer) { 316 double t_end = Ioss::Utils::timer(); 317 double duration = util().global_minmax(t_end - t_begin, Ioss::ParallelUtils::DO_MAX); 318 if (myProcessor == 0) { 319 fmt::print(Ioss::DEBUG(), "{} File Close Time = {}\n", is_input() ? "Input" : "Output", 320 duration); 321 } 322 } 323 closeDW(); 324 m_cgnsFilePtr = -1; 325 } 326 } 327 finalize_database()328 void ParallelDatabaseIO::finalize_database() const 329 { 330 if (is_input()) { 331 return; 332 } 333 334 if (m_timesteps.empty()) { 335 return; 336 } 337 338 if (!m_dbFinalized) { 339 int file_ptr; 340 if (get_file_per_state()) { 341 file_ptr = m_cgnsBasePtr; 342 } 343 else { 344 file_ptr = get_file_pointer(); 345 } 346 Utils::finalize_database(file_ptr, m_timesteps, get_region(), myProcessor, true); 347 m_dbFinalized = true; 348 } 349 } 350 release_memory__()351 void ParallelDatabaseIO::release_memory__() 352 { 353 nodeMap.release_memory(); 354 elemMap.release_memory(); 355 try { 356 decomp.reset(); 357 } 358 catch (...) { 359 } 360 } 361 node_global_to_local__(int64_t global,bool)362 int64_t ParallelDatabaseIO::node_global_to_local__(int64_t global, bool /* must_exist */) const 363 { 364 // TODO: Fix 365 return global; 366 } 367 element_global_to_local__(int64_t global)368 int64_t ParallelDatabaseIO::element_global_to_local__(int64_t global) const 369 { 370 // TODO: Fix 371 return global; 372 } 373 read_meta_data__()374 void ParallelDatabaseIO::read_meta_data__() 375 { 376 openDatabase__(); 377 378 // Determine the number of bases in the grid. 379 // Currently only handle 1. 380 int n_bases = 0; 381 CGCHECKM(cg_nbases(get_file_pointer(), &n_bases)); 382 if (n_bases != 1) { 383 std::ostringstream errmsg; 384 fmt::print(errmsg, 385 "CGNS: Too many bases ({}); only support files with a single bases at this time", 386 n_bases); 387 IOSS_ERROR(errmsg); 388 } 389 390 get_step_times__(); 391 392 if (open_create_behavior() == Ioss::DB_APPEND) { 393 return; 394 } 395 396 m_meshType = Utils::check_mesh_type(get_file_pointer()); 397 398 // In CGNS, there are duplicated nodes at block boundaries. 399 // We typically only want to retain one copy of these and ignore the other. 400 properties.add(Ioss::Property("RETAIN_FREE_NODES", "NO")); 401 402 if (int_byte_size_api() == 8) { 403 decomp = std::unique_ptr<DecompositionDataBase>( 404 new DecompositionData<int64_t>(properties, util().communicator())); 405 } 406 else { 407 decomp = std::unique_ptr<DecompositionDataBase>( 408 new DecompositionData<int>(properties, util().communicator())); 409 } 410 assert(decomp != nullptr); 411 decomp->decompose_model(get_file_pointer(), m_meshType); 412 413 // ======================================================================== 414 // Get the number of assemblies in the mesh... 415 // Will be the 'families' that contain nodes of 'FamVC_*' 416 Utils::add_assemblies(get_file_pointer(), this); 417 418 if (m_meshType == Ioss::MeshType::STRUCTURED) { 419 handle_structured_blocks(); 420 } 421 else if (m_meshType == Ioss::MeshType::UNSTRUCTURED) { 422 handle_unstructured_blocks(); 423 } 424 #if IOSS_ENABLE_HYBRID 425 else if (mesh_type == Ioss::MeshType::HYBRID) { 426 } 427 #endif 428 else { 429 std::ostringstream errmsg; 430 fmt::print(errmsg, "ERROR: CGNS: Mesh is not Unstructured or Structured which are the only " 431 "types currently supported"); 432 IOSS_ERROR(errmsg); 433 } 434 435 Utils::add_transient_variables(get_file_pointer(), m_timesteps, get_region(), 436 get_field_recognition(), get_field_separator(), myProcessor, 437 true); 438 } 439 handle_unstructured_blocks()440 void ParallelDatabaseIO::handle_unstructured_blocks() 441 { 442 get_region()->property_add( 443 Ioss::Property("global_node_count", (int64_t)decomp->global_node_count())); 444 get_region()->property_add( 445 Ioss::Property("global_element_count", (int64_t)decomp->global_elem_count())); 446 447 nodeCount = decomp->ioss_node_count(); 448 elementCount = decomp->ioss_elem_count(); 449 450 // ======================================================================== 451 // Get the number of families in the mesh... 452 // Will treat these as sidesets if they are of the type "FamilyBC_t" 453 Utils::add_sidesets(get_file_pointer(), this); 454 455 // ======================================================================== 456 // Get the number of zones (element blocks) in the mesh... 457 int base = 1; 458 int i = 0; 459 for (auto &block : decomp->m_elementBlocks) { 460 std::string element_topo = block.topologyType; 461 auto *eblock = new Ioss::ElementBlock(this, block.name(), element_topo, block.ioss_count()); 462 eblock->property_add(Ioss::Property("base", base)); 463 eblock->property_add(Ioss::Property("zone", block.zone())); 464 eblock->property_add(Ioss::Property("id", block.zone())); 465 eblock->property_add(Ioss::Property("guid", util().generate_guid(block.zone()))); 466 eblock->property_add(Ioss::Property("section", block.section())); 467 eblock->property_add(Ioss::Property("original_block_order", i++)); 468 get_region()->add(eblock); 469 #if IOSS_DEBUG_OUTPUT 470 fmt::print(Ioss::DEBUG(), "Added block {}, IOSS topology = '{}' with {} element.\n", 471 block.name(), element_topo, block.ioss_count()); 472 #endif 473 // See if this zone/block is a member of any assemblies... 474 Utils::add_to_assembly(get_file_pointer(), get_region(), eblock, base, block.zone()); 475 } 476 477 // ======================================================================== 478 // Have sidesets, now create sideblocks for each sideset... 479 int id = 0; 480 for (auto &sset : decomp->m_sideSets) { 481 // See if there is an Ioss::SideSet with a matching name... 482 Ioss::SideSet *ioss_sset = get_region()->get_sideset(sset.ss_name()); 483 if (ioss_sset != nullptr) { 484 auto zone = decomp->m_zones[sset.zone()]; 485 std::string block_name = fmt::format("{}/{}", zone.m_name, sset.name()); 486 std::string face_topo = sset.topologyType; 487 #if IOSS_DEBUG_OUTPUT 488 fmt::print(Ioss::DEBUG(), 489 "Processor {}: Added sideblock '{}' of topo {} with {} faces to sset '{}'\n", 490 myProcessor, block_name, face_topo, sset.ioss_count(), ioss_sset->name()); 491 #endif 492 const auto &block = decomp->m_elementBlocks[sset.parentBlockIndex]; 493 494 std::string parent_topo = block.topologyType; 495 Ioss::SideBlock *sblk = 496 new Ioss::SideBlock(this, block_name, face_topo, parent_topo, sset.ioss_count()); 497 sblk->property_add(Ioss::Property("id", id)); 498 sblk->property_add(Ioss::Property("guid", util().generate_guid(id + 1))); 499 sblk->property_add(Ioss::Property("base", 1)); 500 sblk->property_add(Ioss::Property("zone", sset.zone())); 501 sblk->property_add(Ioss::Property("section", sset.section())); 502 Ioss::ElementBlock *eblock = get_region()->get_element_block(block.name()); 503 if (eblock != nullptr) { 504 sblk->set_parent_element_block(eblock); 505 } 506 ioss_sset->add(sblk); 507 } 508 id++; // Really just index into m_sideSets list. 509 } 510 511 auto *nblock = new Ioss::NodeBlock(this, "nodeblock_1", nodeCount, 3); 512 513 nblock->property_add(Ioss::Property("base", base)); 514 get_region()->add(nblock); 515 516 // Create a single node commset 517 auto *commset = 518 new Ioss::CommSet(this, "commset_node", "node", decomp->get_commset_node_size()); 519 commset->property_add(Ioss::Property("id", 1)); 520 commset->property_add(Ioss::Property("guid", util().generate_guid(1))); 521 522 get_region()->add(commset); 523 } 524 finalize_structured_blocks()525 size_t ParallelDatabaseIO::finalize_structured_blocks() 526 { 527 // If there are any Structured blocks, need to iterate them and their 1-to-1 connections 528 // and update the donor_zone id for zones that had not yet been processed at the time of 529 // definition... 530 const auto &blocks = get_region()->get_structured_blocks(); 531 for (auto &block : blocks) { 532 int64_t guid = block->get_property("guid").get_int(); 533 for (auto &conn : block->m_zoneConnectivity) { 534 if (conn.m_donorZone < 0) { 535 auto donor_iter = m_zoneNameMap.find(conn.m_donorName); 536 assert(donor_iter != m_zoneNameMap.end()); 537 conn.m_donorZone = (*donor_iter).second; 538 } 539 conn.m_donorGUID = util().generate_guid(conn.m_donorZone, conn.m_donorProcessor); 540 conn.m_ownerGUID = guid; 541 } 542 } 543 544 size_t num_nodes = Utils::resolve_nodes(*get_region(), myProcessor, true); 545 return num_nodes; 546 } 547 handle_structured_blocks()548 void ParallelDatabaseIO::handle_structured_blocks() 549 { 550 int base = 1; 551 552 Utils::add_sidesets(get_file_pointer(), this); 553 554 char basename[CGNS_MAX_NAME_LENGTH + 1]; 555 int cell_dimension = 0; 556 int phys_dimension = 0; 557 CGCHECKM(cg_base_read(get_file_pointer(), base, basename, &cell_dimension, &phys_dimension)); 558 559 // Iterate all structured blocks and set the intervals to zero 560 // if the m_proc field does not match current processor... 561 const auto &zones = decomp->m_structuredZones; 562 563 for (auto &zone : zones) { 564 if (zone->m_adam == zone) { 565 // This is a "root" zone from the undecomposed mesh... 566 // Now see if there are any non-empty blocks with 567 // this m_adam on this processor. If exists, then create 568 // a StructuredBlock; otherwise, create an empty block. 569 auto block_name = zone->m_name; 570 571 Ioss::StructuredBlock *block = nullptr; 572 Ioss::IJK_t zeros{{0, 0, 0}}; 573 for (auto &pzone : zones) { 574 if (pzone->m_proc == myProcessor && pzone->m_adam == zone) { 575 // Create a non-empty structured block on this processor... 576 block = new Ioss::StructuredBlock(this, block_name, phys_dimension, pzone->m_ordinal, 577 pzone->m_offset, pzone->m_adam->m_ordinal); 578 579 for (auto &zgc : pzone->m_zoneConnectivity) { 580 // Update donor_zone to point to adam zone instead of child. 581 auto dz = zones[zgc.m_donorZone - 1]; 582 assert(dz->m_zone == zgc.m_donorZone); 583 auto oz = zones[zgc.m_ownerZone - 1]; 584 assert(oz->m_zone == zgc.m_ownerZone); 585 zgc.m_donorZone = dz->m_adam->m_zone; 586 zgc.m_ownerZone = oz->m_adam->m_zone; 587 block->m_zoneConnectivity.push_back(zgc); 588 } 589 break; 590 } 591 } 592 if (block == nullptr) { 593 // There is no block on this processor corresponding to the m_adam 594 // block. Create an empty block... 595 block = new Ioss::StructuredBlock(this, block_name, phys_dimension, zeros, zeros, 596 zone->m_adam->m_ordinal); 597 for (auto zgc : zone->m_zoneConnectivity) { 598 zgc.m_isActive = false; 599 // Update donor_zone to point to adam zone instead of child. 600 auto dz = zones[zgc.m_donorZone - 1]; 601 assert(dz->m_zone == zgc.m_donorZone); 602 auto oz = zones[zgc.m_ownerZone - 1]; 603 assert(oz->m_zone == zgc.m_ownerZone); 604 zgc.m_donorZone = dz->m_adam->m_zone; 605 zgc.m_ownerZone = oz->m_adam->m_zone; 606 block->m_zoneConnectivity.push_back(zgc); 607 } 608 } 609 assert(block != nullptr); 610 get_region()->add(block); 611 612 block->property_add(Ioss::Property("base", base)); 613 block->property_add(Ioss::Property("zone", zone->m_adam->m_zone)); 614 block->property_add(Ioss::Property("db_zone", zone->m_adam->m_zone)); 615 block->property_add(Ioss::Property("id", zone->m_adam->m_zone)); 616 int64_t guid = util().generate_guid(zone->m_adam->m_zone); 617 block->property_add(Ioss::Property("guid", guid)); 618 619 // See if this zone/block is a member of any assemblies... 620 Utils::add_to_assembly(get_file_pointer(), get_region(), block, base, zone->m_adam->m_zone); 621 622 #if IOSS_DEBUG_OUTPUT 623 fmt::print(Ioss::DEBUG(), "Added block {}, Structured with ID = {}, GUID = {}\n", 624 block_name, zone->m_adam->m_zone, guid); 625 #endif 626 } 627 } 628 629 // ======================================================================== 630 // Iterate each StructuredBlock, get its zone. For that zone, get the number of 631 // boundary conditions and then iterate those and create sideblocks in the 632 // corresponding sideset. 633 const auto &sbs = get_region()->get_structured_blocks(); 634 for (const auto &block : sbs) { 635 // Handle boundary conditions... 636 Utils::add_structured_boundary_conditions(get_file_pointer(), block, true); 637 } 638 639 size_t node_count = finalize_structured_blocks(); 640 auto * nblock = new Ioss::NodeBlock(this, "nodeblock_1", node_count, phys_dimension); 641 nblock->property_add(Ioss::Property("base", base)); 642 get_region()->add(nblock); 643 } 644 resolve_zone_shared_nodes(const CGNSIntVector & nodes,CGNSIntVector & connectivity_map,size_t & owned_node_count,size_t & owned_node_offset)645 void ParallelDatabaseIO::resolve_zone_shared_nodes(const CGNSIntVector &nodes, 646 CGNSIntVector & connectivity_map, 647 size_t & owned_node_count, 648 size_t & owned_node_offset) const 649 { 650 // Determine number of processors that have nodes for this zone. 651 // Avoid mpi_comm_split call if possible. 652 int have_nodes = nodes.empty() ? 0 : 1; 653 int shared_zone_proc_count = 0; 654 MPI_Allreduce(&have_nodes, &shared_zone_proc_count, 1, Ioss::mpi_type(int(0)), MPI_SUM, 655 util().communicator()); 656 657 if (shared_zone_proc_count <= 1) { 658 // There are no shared nodes in this zone. 659 owned_node_count = nodes.size(); 660 std::iota(connectivity_map.begin(), connectivity_map.end(), 1); 661 return; 662 } 663 664 have_nodes = have_nodes == 0 ? MPI_UNDEFINED : 1; 665 MPI_Comm pcomm; 666 MPI_Comm_split(util().communicator(), have_nodes, myProcessor, &pcomm); 667 668 if (have_nodes == 1) { 669 // From here on down, only processors that have nodes are involved... 670 // This zone has nodes/cells on two or more processors. Need to determine 671 // which nodes are shared. 672 673 Ioss::ParallelUtils pm(pcomm); 674 size_t proc_count = pm.parallel_size(); 675 assert((int)proc_count == shared_zone_proc_count); 676 677 // Distribute each node to an "owning" processor based on its id 678 // and assuming a linear distribution (e.g., if on 3 processors, each 679 // proc will "own" 1/3 of the id range. 680 // nodes is sorted. 681 int64_t min = nodes[0]; 682 int64_t max = nodes.back(); 683 min = pm.global_minmax(min, Ioss::ParallelUtils::DO_MIN); 684 max = pm.global_minmax(max, Ioss::ParallelUtils::DO_MAX); 685 686 std::vector<int> p_count(proc_count); 687 size_t per_proc = (max - min + proc_count) / proc_count; 688 size_t proc = 0; 689 int64_t top = min + per_proc; 690 691 // NOTE: nodes is sorted... 692 for (auto node : nodes) { 693 while (node >= top) { 694 top += per_proc; 695 proc++; 696 } 697 assert(proc < proc_count); 698 p_count[proc]++; 699 } 700 701 // Tell each processor how many nodes it will be getting... 702 // Each processor will be responsible for a subsetted range of the overall range. 703 // This processor, should range from min + my_proc*per_proc to min + (my_proc+1)*per_proc. 704 std::vector<int> r_count(proc_count); 705 MPI_Alltoall(p_count.data(), 1, Ioss::mpi_type(int(0)), r_count.data(), 1, 706 Ioss::mpi_type(int(0)), pcomm); 707 708 std::vector<int> recv_disp(proc_count); 709 std::vector<int> send_disp(proc_count); 710 size_t sums = 0; 711 size_t sumr = 0; 712 for (size_t p = 0; p < proc_count; p++) { 713 recv_disp[p] = sumr; 714 sumr += r_count[p]; 715 716 send_disp[p] = sums; 717 sums += p_count[p]; 718 } 719 CGNSIntVector r_nodes(sumr); 720 Ioss::MY_Alltoallv(nodes, p_count, send_disp, r_nodes, r_count, recv_disp, pcomm); 721 722 // Iterate r_nodes list to find duplicate nodes... 723 auto delta = min + pm.parallel_rank() * per_proc; 724 std::vector<int> dup_nodes(per_proc); 725 for (auto &r_node : r_nodes) { 726 auto n = r_node - delta; 727 assert(n < per_proc); 728 dup_nodes[n]++; 729 if (dup_nodes[n] > 1) { 730 r_node = 0; 731 } 732 } 733 734 // Send filtered list back to original processors -- store in 'u_nodes' 735 // This is set of unique block nodes owned by this processor. 736 // If an entry in 'u_nodes' is 0, then that is a non-owned shared node. 737 CGNSIntVector u_nodes(nodes.size()); 738 Ioss::MY_Alltoallv(r_nodes, r_count, recv_disp, u_nodes, p_count, send_disp, pcomm); 739 740 // Count non-zero entries in u_nodes... 741 int64_t local_node_count = 742 std::count_if(u_nodes.cbegin(), u_nodes.cend(), [](int64_t i) { return i > 0; }); 743 owned_node_count = local_node_count; // Calling code wants to know this 744 745 // Determine offset into the zone node block for each processors "chunk" 746 int64_t local_node_offset = 0; 747 MPI_Exscan(&local_node_count, &local_node_offset, 1, Ioss::mpi_type(local_node_count), 748 MPI_SUM, pcomm); 749 owned_node_offset = local_node_offset; // Calling code wants to know this 750 751 // This generates the position of each owned node in this zone consistent 752 // over all processors that this zone is active on. 753 for (auto &u_node : u_nodes) { 754 if (u_node > 0) { 755 u_node = ++local_node_offset; // 1-based local node id for all owned nodes. 756 } 757 } 758 759 // u_nodes now contains the global -> block-local node map for all owned nodes 760 // on the processor. 761 // The zeroes in u_nodes are shared nodes on the processor boundary. 762 // Resend nodes and u_nodes so can resolve the ids of the shared nodes. 763 CGNSIntVector g_to_zone_local(sumr); 764 Ioss::MY_Alltoallv(nodes, p_count, send_disp, r_nodes, r_count, recv_disp, pcomm); 765 Ioss::MY_Alltoallv(u_nodes, p_count, send_disp, g_to_zone_local, r_count, recv_disp, pcomm); 766 767 // Iterate g_to_zone_local to find a zero entry. 768 for (size_t i = 0; i < g_to_zone_local.size(); i++) { 769 if (g_to_zone_local[i] == 0) { 770 // The global id is r_nodes[i] which must also appear earlier in the list... 771 for (size_t j = 0; j < i; j++) { 772 if (r_nodes[j] == r_nodes[i]) { 773 g_to_zone_local[i] = g_to_zone_local[j]; 774 break; 775 } 776 } 777 } 778 } 779 780 // Now, send updated g_to_zone_local back to original processors... 781 Ioss::MY_Alltoallv(g_to_zone_local, r_count, recv_disp, u_nodes, p_count, send_disp, pcomm); 782 783 // At this point: 784 // 'nodes' contains the global node ids that are referenced in this zones connectivity. 785 // 'u_nodes' contains the zone-local 1-based position of that node in this zones node list. 786 // 787 // 788 #ifndef NDEBUG 789 for (auto u_node : u_nodes) { 790 assert(u_node > 0); 791 } 792 #endif 793 std::swap(connectivity_map, u_nodes); 794 MPI_Comm_free(&pcomm); 795 } 796 } 797 write_meta_data()798 void ParallelDatabaseIO::write_meta_data() 799 { 800 int num_zones = get_region()->get_property("element_block_count").get_int() + 801 get_region()->get_property("structured_block_count").get_int(); 802 m_bcOffset.resize(num_zones + 1); // use 1-based zones... 803 m_zoneOffset.resize(num_zones + 1); // use 1-based zones... 804 805 elementCount = 806 Utils::common_write_meta_data(get_file_pointer(), *get_region(), m_zoneOffset, true); 807 } 808 get_step_times__()809 void ParallelDatabaseIO::get_step_times__() 810 { 811 Utils::get_step_times(get_file_pointer(), m_timesteps, get_region(), timeScaleFactor, 812 myProcessor); 813 } 814 write_adjacency_data()815 void ParallelDatabaseIO::write_adjacency_data() 816 { 817 // Determine adjacency information between unstructured blocks. 818 // If block I is adjacent to block J, then they will share at 819 // least 1 "side" (face 3D or edge 2D). 820 // Currently, assuming they are adjacent if they share at least one node... 821 822 // TODO: All calculations are done on processor 0 instead of being distributed. 823 // this will not scale well... 824 825 const auto &blocks = get_region()->get_element_blocks(); 826 if (blocks.size() <= 1) { 827 return; // No adjacent blocks if only one block... 828 } 829 830 // ================= 831 // Determine the minimum and maximum global node id for each zone. 832 // This will be used When determining whether 2 zones are 833 // connected by checking whether the global id node ranges overlap 834 std::vector<int64_t> zone_min_id(blocks.size() + 1, std::numeric_limits<int64_t>::max()); 835 std::vector<int64_t> zone_max_id(blocks.size() + 1, std::numeric_limits<int64_t>::min()); 836 837 for (const auto &block : blocks) { 838 int zone = block->get_property("zone").get_int(); 839 assert((size_t)zone < blocks.size() + 1); 840 841 const auto &I_map = m_globalToBlockLocalNodeMap[zone]; 842 843 // Get min and max global node id for each zone... 844 if (I_map->size() > 0) { 845 auto min_max = std::minmax_element(std::next(I_map->map().begin()), I_map->map().end()); 846 zone_min_id[zone] = *(min_max.first); 847 zone_max_id[zone] = *(min_max.second); 848 } 849 } 850 851 // Get min/max over all processors for each zone... 852 util().global_array_minmax(zone_min_id, Ioss::ParallelUtils::DO_MIN); 853 util().global_array_minmax(zone_max_id, Ioss::ParallelUtils::DO_MAX); 854 // ================= 855 856 auto node_offset = get_processor_zone_node_offset(); 857 858 // Now iterate the blocks again. If the node ranges overlap, then 859 // there is a possibility that there are contiguous nodes; if the 860 // ranges don't overlap, then no possibility... 861 for (auto I = blocks.cbegin(); I != std::prev(blocks.cend()); I++) { 862 int base = (*I)->get_property("base").get_int(); 863 int zone = (*I)->get_property("zone").get_int(); 864 865 // See how many zone I nodes Proc x has that can potentially 866 // overlap with the zones I+1 to end. This will be all nodes 867 // with id > min(zone_min_id[I+1..]) 868 auto min_I = std::min_element(&zone_min_id[zone + 1], &zone_min_id[blocks.size()]); 869 auto I_nodes = gather_nodes_to_proc0(*m_globalToBlockLocalNodeMap[zone], myProcessor, 870 node_offset[zone - 1], util(), *min_I); 871 872 for (auto J = I + 1; J != blocks.end(); J++) { 873 cgsize_t dzone = (*J)->get_property("zone").get_int(); 874 875 if (zone_min_id[dzone] <= zone_max_id[zone] && zone_max_id[dzone] >= zone_min_id[zone]) { 876 // Now gather all zone J nodes that can potentially overlap 877 // with zone I to processor 0... 878 auto J_nodes = gather_nodes_to_proc0(*m_globalToBlockLocalNodeMap[dzone], myProcessor, 879 node_offset[dzone - 1], util(), zone_min_id[zone], 880 zone_max_id[zone]); 881 GL_IdVector common; 882 if (myProcessor == 0) { 883 common = intersect(I_nodes, J_nodes); 884 885 #if IOSS_DEBUG_OUTPUT 886 fmt::print(Ioss::DEBUG(), "Zone {}: {}, Donor Zone {}: {} Common: {}\n\t", zone, 887 I_nodes.size(), dzone, J_nodes.size(), common.size()); 888 889 for (const auto &p : common) { 890 fmt::print(Ioss::DEBUG(), "{}, ", p.first); 891 } 892 fmt::print(Ioss::DEBUG(), "\n\t"); 893 for (const auto &p : common) { 894 fmt::print(Ioss::DEBUG(), "{}, ", p.second); 895 } 896 fmt::print(Ioss::DEBUG(), "\n"); 897 #endif 898 } 899 900 int size = (int)common.size(); 901 MPI_Bcast(&size, 1, MPI_INT, 0, util().communicator()); 902 903 if (size > 0) { 904 // This 'cg_conn_write' should probably be a parallel 905 // function. Since one does not exist, we output the same 906 // data on all processors. Seems to work, but is klugy. 907 908 common.resize(size); 909 MPI_Bcast(common.data(), 2 * size, MPI_INT, 0, util().communicator()); 910 911 CGNSIntVector point_list; 912 CGNSIntVector point_list_donor; 913 point_list.reserve(common.size()); 914 point_list_donor.reserve(common.size()); 915 916 for (auto pnt : common) { 917 point_list.push_back(pnt.first); 918 point_list_donor.push_back(pnt.second); 919 } 920 921 int gc_idx = 0; 922 std::string name = fmt::format("{}_to_{}", (*I)->name(), (*J)->name()); 923 const auto &d1_name = (*J)->name(); 924 925 CGCHECKM(cg_conn_write(get_file_pointer(), base, zone, name.c_str(), CGNS_ENUMV(Vertex), 926 CGNS_ENUMV(Abutting1to1), CGNS_ENUMV(PointList), point_list.size(), 927 point_list.data(), d1_name.c_str(), CGNS_ENUMV(Unstructured), 928 CGNS_ENUMV(PointListDonor), CGNS_ENUMV(DataTypeNull), point_list_donor.size(), 929 point_list_donor.data(), &gc_idx)); 930 931 name = fmt::format("{}_to_{}", (*J)->name(), (*I)->name()); 932 const auto &d2_name = (*I)->name(); 933 934 CGCHECKM(cg_conn_write(get_file_pointer(), base, dzone, name.c_str(), CGNS_ENUMV(Vertex), 935 CGNS_ENUMV(Abutting1to1), CGNS_ENUMV(PointList), point_list_donor.size(), 936 point_list_donor.data(), d2_name.c_str(), CGNS_ENUMV(Unstructured), 937 CGNS_ENUMV(PointListDonor), CGNS_ENUMV(DataTypeNull), point_list.size(), 938 point_list.data(), &gc_idx)); 939 } 940 } 941 } 942 } 943 } 944 begin__(Ioss::State state)945 bool ParallelDatabaseIO::begin__(Ioss::State state) 946 { 947 dbState = state; 948 return true; 949 } 950 free_state_pointer()951 void ParallelDatabaseIO::free_state_pointer() 952 { 953 // If this is the first state file created, then we need to save a reference 954 // to the base CGNS file so we can update the metadata and create links to 955 // the state files (if we are using the file-per-state option) 956 if (m_cgnsBasePtr < 0) { 957 m_cgnsBasePtr = m_cgnsFilePtr; 958 m_cgnsFilePtr = -1; 959 } 960 closeDatabase__(); 961 } 962 open_state_file(int state)963 void ParallelDatabaseIO::open_state_file(int state) 964 { 965 // Close current state file (if any)... 966 free_state_pointer(); 967 968 // Update filename to append state count... 969 decodedFilename.clear(); 970 971 Ioss::FileInfo db(originalDBFilename); 972 std::string new_filename; 973 if (!db.pathname().empty()) { 974 new_filename += db.pathname() + "/"; 975 } 976 977 new_filename += fmt::format("{}-SolutionAtStep{:05}.{}", db.basename(), state, db.extension()); 978 979 DBFilename = new_filename; 980 981 Iocgns::Utils::write_state_meta_data(get_file_pointer(), *get_region(), true); 982 } 983 end__(Ioss::State state)984 bool ParallelDatabaseIO::end__(Ioss::State state) 985 { 986 // Transitioning out of state 'state' 987 switch (state) { 988 case Ioss::STATE_DEFINE_MODEL: 989 if (!is_input() && open_create_behavior() != Ioss::DB_APPEND && 990 open_create_behavior() != Ioss::DB_MODIFY) { 991 write_meta_data(); 992 } 993 if (!is_input() && open_create_behavior() == Ioss::DB_APPEND) { 994 Utils::update_db_zone_property(m_cgnsFilePtr, get_region(), myProcessor, isParallel, true); 995 } 996 break; 997 case Ioss::STATE_MODEL: 998 if (!is_input() && open_create_behavior() != Ioss::DB_APPEND && 999 open_create_behavior() != Ioss::DB_MODIFY) { 1000 write_adjacency_data(); 1001 } 1002 break; 1003 case Ioss::STATE_DEFINE_TRANSIENT: 1004 if (!is_input() && open_create_behavior() != Ioss::DB_APPEND && 1005 open_create_behavior() != Ioss::DB_MODIFY) { 1006 write_results_meta_data(); 1007 } 1008 break; 1009 default: // ignore everything else... 1010 break; 1011 } 1012 return true; 1013 } 1014 begin_state__(int state,double)1015 bool ParallelDatabaseIO::begin_state__(int state, double /* time */) 1016 { 1017 if (is_input()) { 1018 return true; 1019 } 1020 if (get_file_per_state()) { 1021 // Close current state file (if any); create new state file and output metadata... 1022 open_state_file(state); 1023 write_results_meta_data(); 1024 } 1025 Utils::write_flow_solution_metadata(get_file_pointer(), m_cgnsBasePtr, get_region(), state, 1026 &m_currentVertexSolutionIndex, 1027 &m_currentCellCenterSolutionIndex, true); 1028 m_dbFinalized = false; 1029 return true; 1030 } 1031 end_state__(int state,double time)1032 bool ParallelDatabaseIO::end_state__(int state, double time) 1033 { 1034 if (!is_input()) { 1035 m_timesteps.push_back(time); 1036 assert(m_timesteps.size() == (size_t)state); 1037 } 1038 1039 if (!is_input()) { 1040 bool do_flush = true; 1041 if (m_flushInterval != 1) { 1042 if (m_flushInterval == 0 || state % m_flushInterval != 0) { 1043 do_flush = false; 1044 } 1045 } 1046 1047 if (do_flush) { 1048 flush_database__(); 1049 } 1050 } 1051 1052 return true; 1053 } 1054 flush_database__()1055 void ParallelDatabaseIO::flush_database__() const 1056 { 1057 // For HDF5 files, it looks like we need to close the database between 1058 // writes if we want to have a valid database for external access or 1059 // to protect against a crash corrupting the file. 1060 finalize_database(); 1061 closeDatabase__(); 1062 m_cgnsFilePtr = -2; // Tell openDatabase__ that we want to append 1063 } 1064 get_map(entity_type type)1065 const Ioss::Map &ParallelDatabaseIO::get_map(entity_type type) const 1066 { 1067 if (m_meshType == Ioss::MeshType::UNSTRUCTURED) { 1068 switch (type) { 1069 case entity_type::NODE: { 1070 size_t offset = decomp->decomp_node_offset(); 1071 size_t count = decomp->decomp_node_count(); 1072 return get_map(nodeMap, nodeCount, offset, count, entity_type::NODE); 1073 } 1074 case entity_type::ELEM: { 1075 size_t offset = decomp->decomp_elem_offset(); 1076 size_t count = decomp->decomp_elem_count(); 1077 return get_map(elemMap, elementCount, offset, count, entity_type::ELEM); 1078 } 1079 } 1080 } 1081 else { 1082 assert(1 == 0); 1083 } 1084 std::ostringstream errmsg; 1085 fmt::print(errmsg, "INTERNAL ERROR: Invalid map type. " 1086 "Something is wrong in the Iocgns::ParallelDatabaseIO::get_map() function. " 1087 "Please report.\n"); 1088 IOSS_ERROR(errmsg); 1089 } 1090 get_map(Ioss::Map & entity_map,int64_t entity_count,int64_t file_offset,int64_t file_count,entity_type type)1091 const Ioss::Map &ParallelDatabaseIO::get_map(Ioss::Map &entity_map, int64_t entity_count, 1092 int64_t file_offset, int64_t file_count, 1093 entity_type type) const 1094 1095 { 1096 // Allocate space for node number map and read it in... 1097 // Can be called multiple times, allocate 1 time only 1098 if (entity_map.map().empty()) { 1099 entity_map.set_size(entity_count); 1100 1101 if (is_input()) { 1102 Ioss::MapContainer file_data(file_count); 1103 1104 // For cgns, my file_data is just nodes from file_offset to file_offset+file_count 1105 std::iota(file_data.begin(), file_data.end(), file_offset + 1); 1106 1107 if (type == entity_type::NODE) 1108 decomp->communicate_node_data(file_data.data(), &entity_map.map()[1], 1); 1109 else if (type == entity_type::ELEM) 1110 decomp->communicate_element_data(file_data.data(), &entity_map.map()[1], 1); 1111 1112 // Check for sequential node map. 1113 // If not, build the reverse G2L node map... 1114 entity_map.is_sequential(true); 1115 entity_map.build_reverse_map(); 1116 } 1117 else { 1118 // Output database; entity_map.map not set yet... Build a default map. 1119 entity_map.set_default(entity_count); 1120 } 1121 } 1122 return entity_map; 1123 } 1124 get_field_internal(const Ioss::Region *,const Ioss::Field &,void *,size_t)1125 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::Region * /* reg */, 1126 const Ioss::Field & /* field */, void * /* data */, 1127 size_t /* data_size */) const 1128 { 1129 return -1; 1130 } 1131 get_field_internal(const Ioss::NodeBlock * nb,const Ioss::Field & field,void * data,size_t data_size)1132 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::NodeBlock *nb, 1133 const Ioss::Field &field, void *data, 1134 size_t data_size) const 1135 { 1136 // A CGNS DatabaseIO object can have two "types" of NodeBlocks: 1137 // * The normal "all nodes in the model" NodeBlock as used by Exodus 1138 // * A "nodes in a zone" NodeBlock which contains the subset of nodes 1139 // "owned" by a specific StructuredBlock or ElementBlock zone. 1140 // 1141 // Question: How to determine if the NodeBlock is the "global" Nodeblock 1142 // or a "sub" NodeBlock: Use the "is_nonglobal_nodeblock()" function. 1143 if (nb->is_nonglobal_nodeblock()) { 1144 return get_field_internal_sub_nb(nb, field, data, data_size); 1145 } 1146 1147 size_t num_to_get = field.verify(data_size); 1148 1149 Ioss::Field::RoleType role = field.get_role(); 1150 if (role == Ioss::Field::MESH) { 1151 if (field.get_name() == "mesh_model_coordinates_x" || 1152 field.get_name() == "mesh_model_coordinates_y" || 1153 field.get_name() == "mesh_model_coordinates_z" || 1154 field.get_name() == "mesh_model_coordinates") { 1155 decomp->get_node_coordinates(get_file_pointer(), (double *)data, field); 1156 } 1157 1158 else if (field.get_name() == "ids") { 1159 // Map the local ids in this node block 1160 // (1...node_count) to global node ids. 1161 get_map(entity_type::NODE).map_implicit_data(data, field, num_to_get, 0); 1162 } 1163 // The 1..global_node_count id. In a parallel-decomposed run, 1164 // it maps the node back to its implicit position in the serial 1165 // undecomposed mesh file. This is ONLY provided for backward- 1166 // compatibility and should not be used unless absolutely required. 1167 else if (field.get_name() == "implicit_ids") { 1168 size_t offset = decomp->decomp_node_offset(); 1169 size_t count = decomp->decomp_node_count(); 1170 if (int_byte_size_api() == 4) { 1171 std::vector<int> file_ids(count); 1172 std::iota(file_ids.begin(), file_ids.end(), offset + 1); 1173 decomp->communicate_node_data(file_ids.data(), (int *)data, 1); 1174 } 1175 else { 1176 std::vector<int64_t> file_ids(count); 1177 std::iota(file_ids.begin(), file_ids.end(), offset + 1); 1178 decomp->communicate_node_data(file_ids.data(), (int64_t *)data, 1); 1179 } 1180 } 1181 1182 else if (field.get_name() == "connectivity") { 1183 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 1184 } 1185 else if (field.get_name() == "connectivity_raw") { 1186 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 1187 } 1188 else if (field.get_name() == "owning_processor") { 1189 // If parallel, then set the "locally_owned" property on the nodeblocks. 1190 Ioss::CommSet *css = get_region()->get_commset("commset_node"); 1191 if (int_byte_size_api() == 8) { 1192 auto *idata = static_cast<int64_t *>(data); 1193 std::fill(idata, idata + nodeCount, myProcessor); 1194 1195 // Cannot call: 1196 // `css->get_field_data("entity_processor_raw", ent_proc);` 1197 // directly since it will cause a deadlock (in threaded code), 1198 // expand out into corresponding `get_field_internal` call. 1199 Ioss::Field ep_field = css->get_field("entity_processor_raw"); 1200 std::vector<int64_t> ent_proc(ep_field.raw_count() * 1201 ep_field.raw_storage()->component_count()); 1202 size_t ep_data_size = ent_proc.size() * sizeof(int64_t); 1203 get_field_internal(css, ep_field, ent_proc.data(), ep_data_size); 1204 for (size_t i = 0; i < ent_proc.size(); i += 2) { 1205 int64_t node = ent_proc[i + 0]; 1206 int64_t proc = ent_proc[i + 1]; 1207 if (proc < idata[node - 1]) { 1208 idata[node - 1] = proc; 1209 } 1210 } 1211 } 1212 else { 1213 int *idata = static_cast<int *>(data); 1214 std::fill(idata, idata + nodeCount, myProcessor); 1215 1216 Ioss::Field ep_field = css->get_field("entity_processor_raw"); 1217 std::vector<int> ent_proc(ep_field.raw_count() * 1218 ep_field.raw_storage()->component_count()); 1219 size_t ep_data_size = ent_proc.size() * sizeof(int); 1220 get_field_internal(css, ep_field, ent_proc.data(), ep_data_size); 1221 for (size_t i = 0; i < ent_proc.size(); i += 2) { 1222 int node = ent_proc[i + 0]; 1223 int proc = ent_proc[i + 1]; 1224 if (proc < idata[node - 1]) { 1225 idata[node - 1] = proc; 1226 } 1227 } 1228 } 1229 } 1230 else { 1231 num_to_get = Ioss::Utils::field_warning(nb, field, "input"); 1232 } 1233 } 1234 else if (role == Ioss::Field::TRANSIENT) { 1235 // Locate the FlowSolution node corresponding to the correct state/step/time 1236 // TODO: do this at read_meta_data() and store... 1237 int step = get_region()->get_current_state(); 1238 auto var_type = field.transformed_storage(); 1239 int comp_count = var_type->component_count(); 1240 1241 if (comp_count == 1) { 1242 decomp->get_node_field(get_file_pointer(), step, Utils::index(field), (double *)data); 1243 } 1244 else { 1245 std::vector<double> ioss_tmp(num_to_get); 1246 for (int i = 0; i < comp_count; i++) { 1247 decomp->get_node_field(get_file_pointer(), step, Utils::index(field) + i, 1248 ioss_tmp.data()); 1249 1250 size_t index = i; 1251 auto * rdata = static_cast<double *>(data); 1252 for (size_t j = 0; j < num_to_get; j++) { 1253 rdata[index] = ioss_tmp[j]; 1254 index += comp_count; 1255 } 1256 } 1257 } 1258 } 1259 else { 1260 num_to_get = Ioss::Utils::field_warning(nb, field, "input"); 1261 } 1262 return num_to_get; 1263 } 1264 get_field_internal_sub_nb(const Ioss::NodeBlock * nb,const Ioss::Field & field,void * data,size_t data_size)1265 int64_t ParallelDatabaseIO::get_field_internal_sub_nb(const Ioss::NodeBlock *nb, 1266 const Ioss::Field &field, void *data, 1267 size_t data_size) const 1268 { 1269 // Reads field data on a NodeBlock which is a "sub" NodeBlock -- contains the nodes for a 1270 // StructuredBlock instead of for the entire model. 1271 // Currently only TRANSIENT fields are input this way. No valid reason, but that is the current 1272 // use case. 1273 1274 // Get the StructuredBlock that this NodeBlock is contained in: 1275 const Ioss::GroupingEntity *sb = nb->contained_in(); 1276 int base = 1; 1277 int zone = Iocgns::Utils::get_db_zone(sb); 1278 cgsize_t num_to_get = field.verify(data_size); 1279 1280 Ioss::Field::RoleType role = field.get_role(); 1281 if (role == Ioss::Field::TRANSIENT) { 1282 // Locate the FlowSolution node corresponding to the correct state/step/time 1283 // TODO: do this at read_meta_data() and store... 1284 int step = get_region()->get_current_state(); 1285 1286 int solution_index = 1287 Utils::find_solution_index(get_file_pointer(), base, zone, step, CGNS_ENUMV(Vertex)); 1288 1289 auto *rdata = static_cast<double *>(data); 1290 assert(num_to_get == sb->get_property("node_count").get_int()); 1291 cgsize_t rmin[3] = {0, 0, 0}; 1292 cgsize_t rmax[3] = {0, 0, 0}; 1293 if (num_to_get > 0) { 1294 rmin[0] = sb->get_property("offset_i").get_int() + 1; 1295 rmin[1] = sb->get_property("offset_j").get_int() + 1; 1296 rmin[2] = sb->get_property("offset_k").get_int() + 1; 1297 1298 rmax[0] = rmin[0] + sb->get_property("ni").get_int(); 1299 rmax[1] = rmin[1] + sb->get_property("nj").get_int(); 1300 rmax[2] = rmin[2] + sb->get_property("nk").get_int(); 1301 1302 assert(num_to_get == 1303 (rmax[0] - rmin[0] + 1) * (rmax[1] - rmin[1] + 1) * (rmax[2] - rmin[2] + 1)); 1304 } 1305 1306 auto var_type = field.transformed_storage(); 1307 int comp_count = var_type->component_count(); 1308 char field_suffix_separator = get_field_separator(); 1309 1310 if (comp_count == 1) { 1311 CGCHECKM(cg_field_read(get_file_pointer(), base, zone, solution_index, 1312 field.get_name().c_str(), CGNS_ENUMV(RealDouble), rmin, rmax, rdata)); 1313 } 1314 else { 1315 std::vector<double> cgns_data(num_to_get); 1316 for (int i = 0; i < comp_count; i++) { 1317 std::string var_name = 1318 var_type->label_name(field.get_name(), i + 1, field_suffix_separator); 1319 1320 CGCHECKM(cg_field_read(get_file_pointer(), base, zone, solution_index, var_name.c_str(), 1321 CGNS_ENUMV(RealDouble), rmin, rmax, cgns_data.data())); 1322 for (cgsize_t j = 0; j < num_to_get; j++) { 1323 rdata[comp_count * j + i] = cgns_data[j]; 1324 } 1325 } 1326 } 1327 } 1328 // Ignoring all other field role types... 1329 return num_to_get; 1330 } 1331 get_field_internal(const Ioss::EdgeBlock *,const Ioss::Field &,void *,size_t)1332 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::EdgeBlock * /* nb */, 1333 const Ioss::Field & /* field */, void * /* data */, 1334 size_t /* data_size */) const 1335 { 1336 return -1; 1337 } get_field_internal(const Ioss::FaceBlock *,const Ioss::Field &,void *,size_t)1338 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::FaceBlock * /* nb */, 1339 const Ioss::Field & /* field */, void * /* data */, 1340 size_t /* data_size */) const 1341 { 1342 return -1; 1343 } 1344 get_field_internal(const Ioss::StructuredBlock * sb,const Ioss::Field & field,void * data,size_t data_size)1345 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::StructuredBlock *sb, 1346 const Ioss::Field &field, void *data, 1347 size_t data_size) const 1348 { 1349 Ioss::Field::RoleType role = field.get_role(); 1350 cgsize_t base = sb->get_property("base").get_int(); 1351 cgsize_t zone = sb->get_property("zone").get_int(); 1352 1353 cgsize_t num_to_get = field.verify(data_size); 1354 1355 cgsize_t rmin[3] = {0, 0, 0}; 1356 cgsize_t rmax[3] = {0, 0, 0}; 1357 1358 bool cell_field = Utils::is_cell_field(field); 1359 if (cell_field) { 1360 assert(num_to_get == sb->get_property("cell_count").get_int()); 1361 if (num_to_get > 0) { 1362 rmin[0] = sb->get_property("offset_i").get_int() + 1; 1363 rmin[1] = sb->get_property("offset_j").get_int() + 1; 1364 rmin[2] = sb->get_property("offset_k").get_int() + 1; 1365 1366 rmax[0] = rmin[0] + sb->get_property("ni").get_int() - 1; 1367 rmax[1] = rmin[1] + sb->get_property("nj").get_int() - 1; 1368 rmax[2] = rmin[2] + sb->get_property("nk").get_int() - 1; 1369 } 1370 } 1371 else { 1372 // cell nodal field. 1373 assert(num_to_get == sb->get_property("node_count").get_int()); 1374 if (num_to_get > 0) { 1375 rmin[0] = sb->get_property("offset_i").get_int() + 1; 1376 rmin[1] = sb->get_property("offset_j").get_int() + 1; 1377 rmin[2] = sb->get_property("offset_k").get_int() + 1; 1378 1379 rmax[0] = rmin[0] + sb->get_property("ni").get_int(); 1380 rmax[1] = rmin[1] + sb->get_property("nj").get_int(); 1381 rmax[2] = rmin[2] + sb->get_property("nk").get_int(); 1382 } 1383 } 1384 1385 assert(num_to_get == 0 || num_to_get == (rmax[0] - rmin[0] + 1) * (rmax[1] - rmin[1] + 1) * 1386 (rmax[2] - rmin[2] + 1)); 1387 double *rdata = num_to_get > 0 ? static_cast<double *>(data) : nullptr; 1388 1389 if (role == Ioss::Field::MESH) { 1390 1391 if (field.get_name() == "mesh_model_coordinates_x") { 1392 CGCHECKM(cgp_coord_read_data(get_file_pointer(), base, zone, 1, rmin, rmax, rdata)); 1393 } 1394 1395 else if (field.get_name() == "mesh_model_coordinates_y") { 1396 CGCHECKM(cgp_coord_read_data(get_file_pointer(), base, zone, 2, rmin, rmax, rdata)); 1397 } 1398 1399 else if (field.get_name() == "mesh_model_coordinates_z") { 1400 CGCHECKM(cgp_coord_read_data(get_file_pointer(), base, zone, 3, rmin, rmax, rdata)); 1401 } 1402 1403 else if (field.get_name() == "mesh_model_coordinates") { 1404 char basename[CGNS_MAX_NAME_LENGTH + 1]; 1405 int cell_dimension = 0; 1406 int phys_dimension = 0; 1407 CGCHECKM( 1408 cg_base_read(get_file_pointer(), base, basename, &cell_dimension, &phys_dimension)); 1409 1410 // Data required by upper classes store x0, y0, z0, ... xn, 1411 // yn, zn. Data stored in cgns file is x0, ..., xn, y0, 1412 // ..., yn, z0, ..., zn so we have to allocate some scratch 1413 // memory to read in the data and then map into supplied 1414 // 'data' 1415 std::vector<double> coord(num_to_get); 1416 CGCHECKM(cgp_coord_read_data(get_file_pointer(), base, zone, 1, rmin, rmax, coord.data())); 1417 1418 // Map to global coordinate position... 1419 for (cgsize_t i = 0; i < num_to_get; i++) { 1420 rdata[phys_dimension * i + 0] = coord[i]; 1421 } 1422 1423 if (phys_dimension >= 2) { 1424 CGCHECKM( 1425 cgp_coord_read_data(get_file_pointer(), base, zone, 2, rmin, rmax, coord.data())); 1426 1427 // Map to global coordinate position... 1428 for (cgsize_t i = 0; i < num_to_get; i++) { 1429 rdata[phys_dimension * i + 1] = coord[i]; 1430 } 1431 } 1432 1433 if (phys_dimension == 3) { 1434 CGCHECKM( 1435 cgp_coord_read_data(get_file_pointer(), base, zone, 3, rmin, rmax, coord.data())); 1436 1437 // Map to global coordinate position... 1438 for (cgsize_t i = 0; i < num_to_get; i++) { 1439 rdata[phys_dimension * i + 2] = coord[i]; 1440 } 1441 } 1442 } 1443 else if (field.get_name() == "cell_node_ids") { 1444 if (field.get_type() == Ioss::Field::INT64) { 1445 auto *idata = static_cast<int64_t *>(data); 1446 sb->get_cell_node_ids(idata, true); 1447 } 1448 else { 1449 assert(field.get_type() == Ioss::Field::INT32); 1450 int *idata = static_cast<int *>(data); 1451 sb->get_cell_node_ids(idata, true); 1452 } 1453 } 1454 else if (field.get_name() == "cell_ids") { 1455 if (field.get_type() == Ioss::Field::INT64) { 1456 auto *idata = static_cast<int64_t *>(data); 1457 sb->get_cell_ids(idata, true); 1458 } 1459 else { 1460 assert(field.get_type() == Ioss::Field::INT32); 1461 int *idata = static_cast<int *>(data); 1462 sb->get_cell_ids(idata, true); 1463 } 1464 } 1465 else { 1466 num_to_get = Ioss::Utils::field_warning(sb, field, "input"); 1467 } 1468 } 1469 else if (role == Ioss::Field::TRANSIENT) { 1470 auto var_type = field.transformed_storage(); 1471 int comp_count = var_type->component_count(); 1472 1473 int sol_index = 0; 1474 int step = get_region()->get_current_state(); 1475 if (cell_field) { 1476 sol_index = Utils::find_solution_index(get_file_pointer(), base, zone, step, CGNS_ENUMV(CellCenter)); 1477 } 1478 else { 1479 sol_index = Utils::find_solution_index(get_file_pointer(), base, zone, step, CGNS_ENUMV(Vertex)); 1480 } 1481 int field_offset = Utils::index(field); 1482 1483 if (comp_count == 1) { 1484 CGCHECKM(cgp_field_read_data(get_file_pointer(), base, zone, sol_index, field_offset, rmin, 1485 rmax, rdata)); 1486 } 1487 else { 1488 std::vector<double> cgns_data(num_to_get); 1489 for (int i = 0; i < comp_count; i++) { 1490 CGCHECKM(cgp_field_read_data(get_file_pointer(), base, zone, sol_index, field_offset + i, 1491 rmin, rmax, cgns_data.data())); 1492 for (cgsize_t j = 0; j < num_to_get; j++) { 1493 rdata[comp_count * j + i] = cgns_data[j]; 1494 } 1495 } 1496 } 1497 } 1498 else { 1499 num_to_get = Ioss::Utils::field_warning(sb, field, "input"); 1500 } 1501 return num_to_get; 1502 } 1503 get_field_internal(const Ioss::ElementBlock * eb,const Ioss::Field & field,void * data,size_t data_size)1504 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::ElementBlock *eb, 1505 const Ioss::Field &field, void *data, 1506 size_t data_size) const 1507 { 1508 int base = eb->get_property("base").get_int(); 1509 int zone = eb->get_property("zone").get_int(); 1510 size_t num_to_get = field.verify(data_size); 1511 auto role = field.get_role(); 1512 1513 if (role == Ioss::Field::MESH) { 1514 // Handle the MESH fields required for a CGNS file model. 1515 // (The 'genesis' portion) 1516 1517 if (field.get_name() == "connectivity_raw" || field.get_name() == "connectivity") { 1518 1519 // The connectivity is stored in a 1D array. 1520 // The element_node index varies fastest 1521 int order = eb->get_property("original_block_order").get_int(); 1522 decomp->get_block_connectivity(get_file_pointer(), data, order); 1523 if (field.get_type() == Ioss::Field::INT32) { 1524 auto *idata = reinterpret_cast<int *>(data); 1525 Utils::map_cgns_connectivity(eb->topology(), num_to_get, idata); 1526 } 1527 else { 1528 auto *idata = reinterpret_cast<int64_t *>(data); 1529 Utils::map_cgns_connectivity(eb->topology(), num_to_get, idata); 1530 } 1531 } 1532 else if (field.get_name() == "ids" || field.get_name() == "implicit_ids") { 1533 // Map the local ids in this element block 1534 // (1..element_count) to global element ids. 1535 get_map(entity_type::ELEM).map_implicit_data(data, field, num_to_get, eb->get_offset()); 1536 } 1537 else { 1538 num_to_get = Ioss::Utils::field_warning(eb, field, "input"); 1539 } 1540 } 1541 else if (role == Ioss::Field::TRANSIENT) { 1542 // Locate the FlowSolution node corresponding to the correct state/step/time 1543 // TODO: do this at read_meta_data() and store... 1544 int step = get_region()->get_current_state(); 1545 int solution_index = 1546 Utils::find_solution_index(get_file_pointer(), base, zone, step, CGNS_ENUMV(CellCenter)); 1547 1548 int order = eb->get_property("original_block_order").get_int(); 1549 1550 const Ioss::VariableType *var_type = field.raw_storage(); 1551 1552 // Read into a double variable 1553 // TODO: Support other field types... 1554 size_t num_entity = eb->entity_count(); 1555 std::vector<double> temp(num_entity); 1556 1557 // get number of components, cycle through each component 1558 size_t comp_count = var_type->component_count(); 1559 for (size_t i = 0; i < comp_count; i++) { 1560 int field_offset = Utils::index(field) + i; 1561 decomp->get_element_field(get_file_pointer(), solution_index, order, field_offset, 1562 temp.data()); 1563 1564 // Transfer to 'data' array. 1565 size_t k = 0; 1566 assert(field.get_type() == Ioss::Field::REAL); 1567 auto *rvar = static_cast<double *>(data); 1568 for (size_t j = i; j < num_entity * comp_count; j += comp_count) { 1569 rvar[j] = temp[k++]; 1570 } 1571 assert(k == num_entity); 1572 } 1573 } 1574 else { 1575 num_to_get = Ioss::Utils::field_warning(eb, field, "unknown"); 1576 } 1577 return num_to_get; 1578 } 1579 get_field_internal(const Ioss::NodeSet *,const Ioss::Field &,void *,size_t)1580 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::NodeSet * /* ns */, 1581 const Ioss::Field & /* field */, void * /* data */, 1582 size_t /* data_size */) const 1583 { 1584 return -1; 1585 } get_field_internal(const Ioss::EdgeSet *,const Ioss::Field &,void *,size_t)1586 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::EdgeSet * /* ns */, 1587 const Ioss::Field & /* field */, void * /* data */, 1588 size_t /* data_size */) const 1589 { 1590 return -1; 1591 } get_field_internal(const Ioss::FaceSet *,const Ioss::Field &,void *,size_t)1592 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::FaceSet * /* ns */, 1593 const Ioss::Field & /* field */, void * /* data */, 1594 size_t /* data_size */) const 1595 { 1596 return -1; 1597 } get_field_internal(const Ioss::ElementSet *,const Ioss::Field &,void *,size_t)1598 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::ElementSet * /* ns */, 1599 const Ioss::Field & /* field */, void * /* data */, 1600 size_t /* data_size */) const 1601 { 1602 return -1; 1603 } get_field_internal(const Ioss::SideBlock * sb,const Ioss::Field & field,void * data,size_t data_size)1604 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::SideBlock *sb, 1605 const Ioss::Field &field, void *data, 1606 size_t data_size) const 1607 { 1608 int id = sb->get_property("id").get_int(); 1609 const auto &sset = decomp->m_sideSets[id]; 1610 1611 ssize_t num_to_get = field.verify(data_size); 1612 if (num_to_get > 0) { 1613 int64_t entity_count = sb->entity_count(); 1614 if (num_to_get != entity_count) { 1615 std::ostringstream errmsg; 1616 fmt::print(errmsg, "ERROR: Partial field input not yet implemented for side blocks"); 1617 IOSS_ERROR(errmsg); 1618 } 1619 } 1620 1621 Ioss::Field::RoleType role = field.get_role(); 1622 if (role == Ioss::Field::MESH) { 1623 if (field.get_name() == "element_side_raw" || field.get_name() == "element_side") { 1624 1625 decomp->get_sideset_element_side(get_file_pointer(), sset, data); 1626 return num_to_get; 1627 } 1628 else { 1629 num_to_get = Ioss::Utils::field_warning(sb, field, "input"); 1630 } 1631 } 1632 else { 1633 num_to_get = Ioss::Utils::field_warning(sb, field, "input"); 1634 } 1635 return num_to_get; 1636 } 1637 get_field_internal(const Ioss::SideSet *,const Ioss::Field &,void *,size_t)1638 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::SideSet * /* fs */, 1639 const Ioss::Field & /* field */, void * /* data */, 1640 size_t /* data_size */) const 1641 { 1642 return -1; 1643 } get_field_internal(const Ioss::CommSet * cs,const Ioss::Field & field,void * data,size_t data_size)1644 int64_t ParallelDatabaseIO::get_field_internal(const Ioss::CommSet *cs, const Ioss::Field &field, 1645 void *data, size_t data_size) const 1646 { 1647 size_t num_to_get = field.verify(data_size); 1648 1649 // Return the <entity (node or side), processor> pair 1650 if (field.get_name() == "entity_processor" || field.get_name() == "entity_processor_raw") { 1651 1652 // Check type -- node or side 1653 std::string type = cs->get_property("entity_type").get_string(); 1654 1655 if (type == "node") { 1656 1657 bool do_map = field.get_name() == "entity_processor"; 1658 // Convert local node id to global node id and store in 'data' 1659 const Ioss::MapContainer &map = get_map(entity_type::NODE).map(); 1660 if (int_byte_size_api() == 4) { 1661 decomp->get_node_entity_proc_data(static_cast<int *>(data), map, do_map); 1662 } 1663 else { 1664 decomp->get_node_entity_proc_data(static_cast<int64_t *>(data), map, do_map); 1665 } 1666 } 1667 else { 1668 std::ostringstream errmsg; 1669 fmt::print(errmsg, "ERROR: Invalid commset type {}", type); 1670 IOSS_ERROR(errmsg); 1671 } 1672 } 1673 else if (field.get_name() == "ids") { 1674 // Do nothing, just handles an idiosyncrasy of the GroupingEntity 1675 } 1676 else { 1677 num_to_get = Ioss::Utils::field_warning(cs, field, "input"); 1678 } 1679 return num_to_get; 1680 } 1681 handle_element_ids(const Ioss::ElementBlock * eb,void * ids,size_t num_to_get,size_t offset,size_t count)1682 int64_t ParallelDatabaseIO::handle_element_ids(const Ioss::ElementBlock *eb, void *ids, 1683 size_t num_to_get, size_t offset, 1684 size_t count) const 1685 { 1686 bool in_define = (dbState == Ioss::STATE_MODEL) || (dbState == Ioss::STATE_DEFINE_MODEL); 1687 if (in_define) { 1688 if (m_elemGlobalImplicitMap.empty()) { 1689 m_elemGlobalImplicitMap.resize(elementCount); 1690 } 1691 // Build the implicit_global map used to map an elements 1692 // local-implicit position to the global-implicit 1693 // position. Primarily used for sideset elements. 1694 // Elements starting at 'eb_offset' map to the global implicit 1695 // position of 'offset' 1696 int64_t eb_offset = eb->get_offset(); 1697 for (size_t i = 0; i < count; i++) { 1698 m_elemGlobalImplicitMap[eb_offset + i] = offset + i + 1; 1699 } 1700 } 1701 1702 elemMap.set_size(elementCount); 1703 int64_t eb_offset = eb->get_offset(); 1704 if (int_byte_size_api() == 4) { 1705 elemMap.set_map(static_cast<int *>(ids), num_to_get, eb_offset, in_define); 1706 } 1707 else { 1708 elemMap.set_map(static_cast<int64_t *>(ids), num_to_get, eb_offset, in_define); 1709 } 1710 return num_to_get; 1711 } 1712 put_field_internal(const Ioss::Region *,const Ioss::Field &,void *,size_t)1713 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::Region * /* region */, 1714 const Ioss::Field & /* field */, void * /* data */, 1715 size_t /* data_size */) const 1716 { 1717 return -1; 1718 } 1719 put_field_internal(const Ioss::NodeBlock * nb,const Ioss::Field & field,void * data,size_t data_size)1720 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::NodeBlock *nb, 1721 const Ioss::Field &field, void *data, 1722 size_t data_size) const 1723 { 1724 // A CGNS DatabaseIO object can have two "types" of NodeBlocks: 1725 // * The normal "all nodes in the model" NodeBlock as used by Exodus 1726 // * A "nodes in a zone" NodeBlock which contains the subset of nodes 1727 // "owned" by a specific StructuredBlock or ElementBlock zone. 1728 // 1729 // Question: How to determine if the NodeBlock is the "global" Nodeblock 1730 // or a "sub" NodeBlock: Use the "is_nonglobal_nodeblock()" function. 1731 if (nb->is_nonglobal_nodeblock()) { 1732 return put_field_internal_sub_nb(nb, field, data, data_size); 1733 } 1734 1735 Ioss::Field::RoleType role = field.get_role(); 1736 cgsize_t base = 1; 1737 size_t num_to_get = field.verify(data_size); 1738 1739 // Instead of outputting a global nodeblock's worth of data, 1740 // the data is output a "zone" at a time. 1741 // The m_globalToBlockLocalNodeMap[zone] map is used (Ioss::Map pointer) 1742 // This map is built during the output of block connectivity, 1743 // so for cgns unstructured mesh, we need to output ElementBlock connectivity 1744 // prior to outputting nodal coordinates. 1745 for (const auto &z : m_globalToBlockLocalNodeMap) { 1746 if (z.second == nullptr) { 1747 std::ostringstream errmsg; 1748 fmt::print(errmsg, 1749 "ERROR: CGNS: The globalToBlockLocalNodeMap is not defined, so nodal fields " 1750 "cannot be output."); 1751 IOSS_ERROR(errmsg); 1752 } 1753 } 1754 1755 if (role == Ioss::Field::MESH) { 1756 if (field.get_name() == "ids") { 1757 // The ids coming in are the global ids; their position is the 1758 // local id -1 (That is, data[0] contains the global id of local 1759 // node 1) 1760 handle_node_ids(data, num_to_get); 1761 } 1762 else if (field.get_name() == "mesh_model_coordinates" || 1763 field.get_name() == "mesh_model_coordinates_x" || 1764 field.get_name() == "mesh_model_coordinates_y" || 1765 field.get_name() == "mesh_model_coordinates_z") { 1766 auto *rdata = static_cast<double *>(data); 1767 1768 std::vector<int64_t> node_offset = get_processor_zone_node_offset(); 1769 1770 if (field.get_name() == "mesh_model_coordinates") { 1771 int spatial_dim = nb->get_property("component_degree").get_int(); 1772 1773 // Output all coordinates, a zone's worth of data at a time... 1774 1775 for (const auto &block : m_globalToBlockLocalNodeMap) { 1776 auto zone = block.first; 1777 // NOTE: 'block_map' has one more entry than node_count. First entry is for something 1778 // else. But, ->size() returns correct size (ignoring first entry) 1779 // 'block_map' is 1-based. 1780 const auto & block_map = block.second; 1781 std::vector<double> x(block_map->size()); 1782 std::vector<double> y(block_map->size()); 1783 std::vector<double> z(block_map->size()); 1784 1785 for (size_t i = 0; i < block_map->size(); i++) { 1786 auto global = block_map->map()[i + 1]; 1787 auto local = nodeMap.global_to_local(global) - 1; 1788 assert(local >= 0 && local < (int64_t)num_to_get); 1789 1790 x[i] = rdata[local * spatial_dim + 0]; 1791 if (spatial_dim > 1) { 1792 y[i] = rdata[local * spatial_dim + 1]; 1793 } 1794 if (spatial_dim > 2) { 1795 z[i] = rdata[local * spatial_dim + 2]; 1796 } 1797 } 1798 1799 // Create the zone 1800 // Output this zones coordinates... 1801 int crd_idx = 0; 1802 CGCHECKM(cgp_coord_write(get_file_pointer(), base, zone, CGNS_ENUMV(RealDouble), "CoordinateX", 1803 &crd_idx)); 1804 cgsize_t start = node_offset[zone - 1] + 1; 1805 cgsize_t finish = start + block_map->size() - 1; 1806 1807 auto xx = block_map->size() > 0 ? x.data() : nullptr; 1808 CGCHECKM( 1809 cgp_coord_write_data(get_file_pointer(), base, zone, crd_idx, &start, &finish, xx)); 1810 1811 if (spatial_dim > 1) { 1812 CGCHECKM(cgp_coord_write(get_file_pointer(), base, zone, CGNS_ENUMV(RealDouble), "CoordinateY", 1813 &crd_idx)); 1814 auto yy = block_map->size() > 0 ? y.data() : nullptr; 1815 CGCHECKM(cgp_coord_write_data(get_file_pointer(), base, zone, crd_idx, &start, 1816 &finish, yy)); 1817 } 1818 1819 if (spatial_dim > 2) { 1820 CGCHECKM(cgp_coord_write(get_file_pointer(), base, zone, CGNS_ENUMV(RealDouble), "CoordinateZ", 1821 &crd_idx)); 1822 auto zz = block_map->size() > 0 ? z.data() : nullptr; 1823 CGCHECKM(cgp_coord_write_data(get_file_pointer(), base, zone, crd_idx, &start, 1824 &finish, zz)); 1825 } 1826 } 1827 } 1828 else { 1829 // Outputting only a single coordinate value... 1830 for (const auto &block : m_globalToBlockLocalNodeMap) { 1831 auto zone = block.first; 1832 // NOTE: 'block_map' has one more entry than node_count. First entry is for something 1833 // else. 1834 // 'block_map' is 1-based. 1835 const auto & block_map = block.second; 1836 std::vector<double> xyz(block_map->size()); 1837 1838 for (size_t i = 0; i < block_map->size(); i++) { 1839 auto global = block_map->map()[i + 1]; 1840 auto local = nodeMap.global_to_local(global) - 1; 1841 xyz[i] = rdata[local]; 1842 } 1843 1844 std::string cgns_name = "Invalid"; 1845 if (field.get_name() == "mesh_model_coordinates_x") { 1846 cgns_name = "CoordinateX"; 1847 } 1848 else if (field.get_name() == "mesh_model_coordinates_y") { 1849 cgns_name = "CoordinateY"; 1850 } 1851 else if (field.get_name() == "mesh_model_coordinates_z") { 1852 cgns_name = "CoordinateZ"; 1853 } 1854 // Create the zone 1855 // Output this zones coordinates... 1856 int crd_idx = 0; 1857 CGCHECKM(cgp_coord_write(get_file_pointer(), base, zone, CGNS_ENUMV(RealDouble), 1858 cgns_name.c_str(), &crd_idx)); 1859 cgsize_t start = node_offset[zone - 1] + 1; 1860 cgsize_t finish = start + block_map->size() - 1; 1861 auto xx = block_map->size() > 0 ? xyz.data() : nullptr; 1862 CGCHECKM( 1863 cgp_coord_write_data(get_file_pointer(), base, zone, crd_idx, &start, &finish, xx)); 1864 } 1865 } 1866 } 1867 else { 1868 num_to_get = Ioss::Utils::field_warning(nb, field, "output"); 1869 } 1870 } 1871 else if (role == Ioss::Field::TRANSIENT) { 1872 // Instead of outputting a global nodeblock's worth of data, 1873 // the data is output a "zone" at a time. 1874 // The m_globalToBlockLocalNodeMap[zone] map is used (Ioss::Map pointer) 1875 // This map is built during the output of block connectivity, 1876 // so for cgns unstructured mesh, we need to output ElementBlock connectivity 1877 // prior to outputting nodal coordinates. 1878 std::vector<int64_t> node_offset = get_processor_zone_node_offset(); 1879 1880 const Ioss::VariableType *var_type = field.raw_storage(); 1881 size_t comp_count = var_type->component_count(); 1882 1883 double *rdata = num_to_get > 0 ? static_cast<double *>(data) : nullptr; 1884 1885 for (const auto &block : m_globalToBlockLocalNodeMap) { 1886 auto zone = block.first; 1887 // NOTE: 'block_map' has one more entry than node_count. 1888 // First entry is for something else. 'block_map' is 1889 // 1-based. 1890 const auto & block_map = block.second; 1891 std::vector<double> blk_data(block_map->size()); 1892 1893 cgsize_t start = node_offset[zone - 1] + 1; 1894 cgsize_t finish = start + block_map->size() - 1; 1895 1896 char field_suffix_separator = get_field_separator(); 1897 1898 for (size_t i = 0; i < comp_count; i++) { 1899 for (size_t j = 0; j < block_map->size(); j++) { 1900 auto global = block_map->map()[j + 1]; 1901 auto local = nodeMap.global_to_local(global) - 1; 1902 assert(local >= 0 && local < (int64_t)num_to_get); 1903 blk_data[j] = rdata[local * comp_count + i]; 1904 } 1905 std::string var_name = (comp_count > 1) ? var_type->label_name(field.get_name(), i + 1, 1906 field_suffix_separator) 1907 : field.get_name(); 1908 int cgns_field = 0; 1909 CGCHECKM(cgp_field_write(get_file_pointer(), base, zone, m_currentVertexSolutionIndex, 1910 CGNS_ENUMV(RealDouble), var_name.c_str(), &cgns_field)); 1911 1912 CGCHECKM(cgp_field_write_data(get_file_pointer(), base, zone, 1913 m_currentVertexSolutionIndex, cgns_field, &start, &finish, 1914 blk_data.data())); 1915 if (i == 0) 1916 Utils::set_field_index(field, cgns_field, CGNS_ENUMV(Vertex)); 1917 } 1918 } 1919 } 1920 else { 1921 num_to_get = Ioss::Utils::field_warning(nb, field, "output"); 1922 } 1923 return num_to_get; 1924 } 1925 put_field_internal_sub_nb(const Ioss::NodeBlock * nb,const Ioss::Field & field,void * data,size_t data_size)1926 int64_t ParallelDatabaseIO::put_field_internal_sub_nb(const Ioss::NodeBlock *nb, 1927 const Ioss::Field &field, void *data, 1928 size_t data_size) const 1929 { 1930 // Outputs field data on a NodeBlock which is a "sub" NodeBlock -- contains the nodes for a 1931 // StructuredBlock instead of for the entire model. 1932 // Currently only TRANSIENT fields are output this way. No valid reason, but that is the 1933 // current use case. 1934 1935 // Get the StructuredBlock that this NodeBlock is contained in: 1936 const Ioss::GroupingEntity *sb = nb->contained_in(); 1937 int zone = Iocgns::Utils::get_db_zone(sb); 1938 cgsize_t num_to_get = field.verify(data_size); 1939 1940 Ioss::Field::RoleType role = field.get_role(); 1941 if (role == Ioss::Field::TRANSIENT) { 1942 int base = 1; 1943 auto *rdata = static_cast<double *>(data); 1944 int cgns_field = 0; 1945 auto var_type = field.transformed_storage(); 1946 int comp_count = var_type->component_count(); 1947 char field_suffix_separator = get_field_separator(); 1948 1949 cgsize_t rmin[3] = {0, 0, 0}; 1950 cgsize_t rmax[3] = {0, 0, 0}; 1951 1952 assert(num_to_get == sb->get_property("node_count").get_int()); 1953 if (num_to_get > 0) { 1954 rmin[0] = sb->get_property("offset_i").get_int() + 1; 1955 rmin[1] = sb->get_property("offset_j").get_int() + 1; 1956 rmin[2] = sb->get_property("offset_k").get_int() + 1; 1957 1958 rmax[0] = rmin[0] + sb->get_property("ni").get_int(); 1959 rmax[1] = rmin[1] + sb->get_property("nj").get_int(); 1960 rmax[2] = rmin[2] + sb->get_property("nk").get_int(); 1961 } 1962 1963 if (comp_count == 1) { 1964 CGCHECKM(cgp_field_write(get_file_pointer(), base, zone, m_currentVertexSolutionIndex, 1965 CGNS_ENUMV(RealDouble), field.get_name().c_str(), &cgns_field)); 1966 Utils::set_field_index(field, cgns_field, CGNS_ENUMV(Vertex)); 1967 1968 CGCHECKM(cgp_field_write_data(get_file_pointer(), base, zone, m_currentVertexSolutionIndex, 1969 cgns_field, rmin, rmax, rdata)); 1970 } 1971 else { 1972 std::vector<double> cgns_data(num_to_get); 1973 for (int i = 0; i < comp_count; i++) { 1974 for (cgsize_t j = 0; j < num_to_get; j++) { 1975 cgns_data[j] = rdata[comp_count * j + i]; 1976 } 1977 std::string var_name = 1978 var_type->label_name(field.get_name(), i + 1, field_suffix_separator); 1979 1980 CGCHECKM(cgp_field_write(get_file_pointer(), base, zone, m_currentVertexSolutionIndex, 1981 CGNS_ENUMV(RealDouble), var_name.c_str(), &cgns_field)); 1982 if (i == 0) { 1983 Utils::set_field_index(field, cgns_field, CGNS_ENUMV(Vertex)); 1984 } 1985 1986 CGCHECKM(cgp_field_write_data(get_file_pointer(), base, zone, 1987 m_currentVertexSolutionIndex, cgns_field, rmin, rmax, 1988 cgns_data.data())); 1989 } 1990 } 1991 } 1992 // Ignoring all other field role types... 1993 return num_to_get; 1994 } 1995 put_field_internal(const Ioss::ElementBlock * eb,const Ioss::Field & field,void * data,size_t data_size)1996 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::ElementBlock *eb, 1997 const Ioss::Field &field, void *data, 1998 size_t data_size) const 1999 { 2000 size_t num_to_get = field.verify(data_size); 2001 2002 Ioss::Field::RoleType role = field.get_role(); 2003 2004 if (role == Ioss::Field::MESH) { 2005 // Handle the MESH fields required for a CGNS file model. 2006 // (The 'genesis' portion) 2007 if (field.get_name() == "ids") { 2008 size_t proc_offset = eb->get_property("proc_offset").get_int(); 2009 handle_element_ids(eb, data, num_to_get, proc_offset, num_to_get); 2010 } 2011 else if (field.get_name() == "connectivity") { 2012 // This blocks zone has not been defined. 2013 // Get the "node block" for this element block... 2014 size_t element_nodes = eb->topology()->number_nodes(); 2015 assert((size_t)field.raw_storage()->component_count() == element_nodes); 2016 2017 CGNSIntVector nodes; 2018 nodes.reserve(element_nodes * num_to_get); 2019 2020 if (field.get_type() == Ioss::Field::INT32) { 2021 int *idata = num_to_get > 0 ? reinterpret_cast<int *>(data) : nullptr; 2022 for (size_t i = 0; i < element_nodes * num_to_get; i++) { 2023 nodes.push_back(idata[i]); 2024 } 2025 } 2026 else { 2027 int64_t *idata = num_to_get > 0 ? reinterpret_cast<int64_t *>(data) : nullptr; 2028 for (size_t i = 0; i < element_nodes * num_to_get; i++) { 2029 nodes.push_back(idata[i]); 2030 } 2031 } 2032 Ioss::Utils::uniquify(nodes); 2033 2034 // Resolve zone-shared nodes (nodes used in this zone, but are 2035 // shared on processor boundaries). 2036 // This routine determines the mapping of each global id node 2037 // in 'nodes' to the zone-local position. 2038 // This mapping is in 'connectivity_map' and is correct for all 2039 // nodes on this processor whether they are owned or shared. 2040 // 2041 // 'resolve_zone_shared_nodes' also returns the number of nodes owned on this 2042 // processor, and the 'offset' of this processors chunk of nodes into the overall 2043 // set of nodes for the zone. Each processors chunk of nodes is contiguous 2044 // 2045 // The 'nodes' and 'connectivity_map' vectors are used later below to generate 2046 // the map of which global node data is written by this processor for this zone. 2047 CGNSIntVector connectivity_map(nodes.size()); 2048 size_t owned_node_count = 0; 2049 size_t owned_node_offset = 0; 2050 resolve_zone_shared_nodes(nodes, connectivity_map, owned_node_count, owned_node_offset); 2051 2052 // Get total count on all processors... 2053 // Note that there will be shared nodes on processor boundaries that need to be 2054 // accounted for... 2055 cgsize_t size[3] = {0, 0, 0}; 2056 size[0] = owned_node_count; 2057 size[1] = eb->entity_count(); 2058 2059 MPI_Allreduce(MPI_IN_PLACE, size, 3, cgns_mpi_type(), MPI_SUM, util().communicator()); 2060 2061 // Now, we have the node count and cell count so we can create a zone... 2062 int base = 1; 2063 int zone = 0; 2064 2065 CGCHECKM(cg_zone_write(get_file_pointer(), base, eb->name().c_str(), size, CGNS_ENUMV(Unstructured), 2066 &zone)); 2067 eb->property_update("zone", zone); 2068 eb->property_update("id", zone); 2069 eb->property_update("guid", util().generate_guid(zone)); 2070 eb->property_update("section", 1); 2071 eb->property_update("base", base); 2072 eb->property_update("zone_node_count", size[0]); 2073 eb->property_update("zone_element_count", size[1]); 2074 2075 if (eb->property_exists("assembly")) { 2076 std::string assembly = eb->get_property("assembly").get_string(); 2077 CGCHECKM(cg_goto(get_file_pointer(), base, "Zone_t", zone, "end")); 2078 CGCHECKM(cg_famname_write(assembly.c_str())); 2079 } 2080 2081 if (size[1] > 0) { 2082 CGNS_ENUMT(ElementType_t) type = Utils::map_topology_to_cgns(eb->topology()->name()); 2083 int sect = 0; 2084 CGCHECKM(cgp_section_write(get_file_pointer(), base, zone, "HexElements", type, 1, 2085 size[1], 0, §)); 2086 2087 int64_t start = 0; 2088 MPI_Exscan(&num_to_get, &start, 1, Ioss::mpi_type(start), MPI_SUM, util().communicator()); 2089 // Of the cells/elements in this zone, this proc handles 2090 // those starting at 'proc_offset+1' to 'proc_offset+num_entity' 2091 eb->property_update("proc_offset", start); 2092 2093 // Map connectivity global ids to zone-local 1-based ids. 2094 CGNSIntVector connect; 2095 connect.reserve(num_to_get * element_nodes); 2096 2097 if (field.get_type() == Ioss::Field::INT32) { 2098 int *idata = num_to_get > 0 ? reinterpret_cast<int *>(data) : nullptr; 2099 for (size_t i = 0; i < num_to_get * element_nodes; i++) { 2100 auto id = idata[i]; 2101 auto iter = std::lower_bound(nodes.cbegin(), nodes.cend(), id); 2102 assert(iter != nodes.end()); 2103 auto cur_pos = iter - nodes.cbegin(); 2104 connect.push_back(connectivity_map[cur_pos]); 2105 } 2106 } 2107 else { 2108 int64_t *idata = num_to_get > 0 ? reinterpret_cast<int64_t *>(data) : nullptr; 2109 for (size_t i = 0; i < num_to_get * element_nodes; i++) { 2110 auto id = idata[i]; 2111 auto iter = std::lower_bound(nodes.cbegin(), nodes.cend(), id); 2112 assert(iter != nodes.cend()); 2113 auto cur_pos = iter - nodes.cbegin(); 2114 connect.push_back(connectivity_map[cur_pos]); 2115 } 2116 } 2117 2118 Utils::unmap_cgns_connectivity(eb->topology(), num_to_get, connect.data()); 2119 CGCHECKM(cgp_elements_write_data(get_file_pointer(), base, zone, sect, start + 1, 2120 start + num_to_get, connect.data())); 2121 2122 int64_t eb_size = num_to_get; 2123 MPI_Allreduce(MPI_IN_PLACE, &eb_size, 1, Ioss::mpi_type(eb_size), MPI_SUM, 2124 util().communicator()); 2125 2126 m_bcOffset[zone] += eb_size; 2127 eb->property_update("section", sect); 2128 } 2129 2130 // The 'nodes' map needs to be updated to filter out any nodes 2131 // that are not owned by this processor. Currently contains both 2132 // owned and shared so we could update the connectivity... 2133 // The 'connectivity_map' value indicates whether it is owned or shared -- 2134 // if 'connectivity_map[i] > owned_node_offset, then it is owned; otherwise shared. 2135 if (!nodes.empty()) { 2136 for (size_t i = 0; i < nodes.size(); i++) { 2137 if (connectivity_map[i] <= (cgsize_t)owned_node_offset) { 2138 nodes[i] = std::numeric_limits<cgsize_t>::max(); 2139 } 2140 } 2141 connectivity_map.clear(); 2142 connectivity_map.shrink_to_fit(); 2143 2144 Ioss::Utils::uniquify(nodes); 2145 if (nodes.back() == std::numeric_limits<cgsize_t>::max()) { 2146 nodes.pop_back(); 2147 } 2148 nodes.shrink_to_fit(); 2149 } 2150 assert(nodes.size() == owned_node_count); 2151 2152 // Now we have a valid zone so can update some data structures... 2153 m_zoneOffset[zone] = m_zoneOffset[zone - 1] + size[1]; 2154 m_globalToBlockLocalNodeMap[zone] = 2155 new Ioss::Map("node", get_filename() + "::" + eb->name(), myProcessor); 2156 m_globalToBlockLocalNodeMap[zone]->map().reserve(nodes.size() + 1); 2157 m_globalToBlockLocalNodeMap[zone]->map().push_back(1); // Non one-to-one map 2158 for (auto i : nodes) { 2159 m_globalToBlockLocalNodeMap[zone]->map().push_back(i); 2160 } 2161 } 2162 else { 2163 num_to_get = Ioss::Utils::field_warning(eb, field, "output"); 2164 } 2165 } 2166 else if (role == Ioss::Field::TRANSIENT) { 2167 double * rdata = num_to_get > 0 ? static_cast<double *>(data) : nullptr; 2168 const Ioss::VariableType *var_type = field.raw_storage(); 2169 2170 int base = eb->get_property("base").get_int(); 2171 int zone = eb->get_property("zone").get_int(); 2172 2173 cgsize_t start = eb->get_property("proc_offset").get_int(); 2174 cgsize_t range_min[1] = {start + 1}; 2175 cgsize_t range_max[1] = {(cgsize_t)start + (cgsize_t)num_to_get}; 2176 2177 // get number of components, cycle through each component 2178 size_t comp_count = var_type->component_count(); 2179 if (comp_count == 1) { 2180 int cgns_field = 0; 2181 CGCHECKM(cgp_field_write(get_file_pointer(), base, zone, m_currentCellCenterSolutionIndex, 2182 CGNS_ENUMV(RealDouble), field.get_name().c_str(), &cgns_field)); 2183 CGCHECKM(cgp_field_write_data(get_file_pointer(), base, zone, 2184 m_currentCellCenterSolutionIndex, cgns_field, range_min, 2185 range_max, rdata)); 2186 Utils::set_field_index(field, cgns_field, CGNS_ENUMV(CellCenter)); 2187 } 2188 else { 2189 char field_suffix_separator = get_field_separator(); 2190 std::vector<double> cgns_data(num_to_get); 2191 for (size_t i = 0; i < comp_count; i++) { 2192 for (size_t j = 0; j < num_to_get; j++) { 2193 cgns_data[j] = rdata[comp_count * j + i]; 2194 } 2195 std::string var_name = 2196 var_type->label_name(field.get_name(), i + 1, field_suffix_separator); 2197 int cgns_field = 0; 2198 CGCHECKM(cgp_field_write(get_file_pointer(), base, zone, m_currentCellCenterSolutionIndex, 2199 CGNS_ENUMV(RealDouble), var_name.c_str(), &cgns_field)); 2200 CGCHECKM(cgp_field_write_data(get_file_pointer(), base, zone, 2201 m_currentCellCenterSolutionIndex, cgns_field, range_min, 2202 range_max, cgns_data.data())); 2203 if (i == 0) { 2204 Utils::set_field_index(field, cgns_field, CGNS_ENUMV(CellCenter)); 2205 } 2206 } 2207 } 2208 } 2209 else { 2210 num_to_get = Ioss::Utils::field_warning(eb, field, "unknown"); 2211 } 2212 return num_to_get; 2213 } 2214 put_field_internal(const Ioss::StructuredBlock * sb,const Ioss::Field & field,void * data,size_t data_size)2215 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::StructuredBlock *sb, 2216 const Ioss::Field &field, void *data, 2217 size_t data_size) const 2218 { 2219 Ioss::Field::RoleType role = field.get_role(); 2220 cgsize_t base = sb->get_property("base").get_int(); 2221 cgsize_t zone = Iocgns::Utils::get_db_zone(sb); 2222 2223 cgsize_t num_to_get = field.verify(data_size); 2224 2225 cgsize_t rmin[3] = {0, 0, 0}; 2226 cgsize_t rmax[3] = {0, 0, 0}; 2227 2228 bool cell_field = Utils::is_cell_field(field); 2229 2230 if (cell_field) { 2231 assert(num_to_get == sb->get_property("cell_count").get_int()); 2232 if (num_to_get > 0) { 2233 rmin[0] = sb->get_property("offset_i").get_int() + 1; 2234 rmin[1] = sb->get_property("offset_j").get_int() + 1; 2235 rmin[2] = sb->get_property("offset_k").get_int() + 1; 2236 2237 rmax[0] = rmin[0] + sb->get_property("ni").get_int() - 1; 2238 rmax[1] = rmin[1] + sb->get_property("nj").get_int() - 1; 2239 rmax[2] = rmin[2] + sb->get_property("nk").get_int() - 1; 2240 } 2241 } 2242 else { 2243 // cell nodal field. 2244 assert(num_to_get == sb->get_property("node_count").get_int()); 2245 if (num_to_get > 0) { 2246 rmin[0] = sb->get_property("offset_i").get_int() + 1; 2247 rmin[1] = sb->get_property("offset_j").get_int() + 1; 2248 rmin[2] = sb->get_property("offset_k").get_int() + 1; 2249 2250 rmax[0] = rmin[0] + sb->get_property("ni").get_int(); 2251 rmax[1] = rmin[1] + sb->get_property("nj").get_int(); 2252 rmax[2] = rmin[2] + sb->get_property("nk").get_int(); 2253 } 2254 } 2255 2256 assert(num_to_get == 0 || num_to_get == (rmax[0] - rmin[0] + 1) * (rmax[1] - rmin[1] + 1) * 2257 (rmax[2] - rmin[2] + 1)); 2258 double *rdata = num_to_get > 0 ? static_cast<double *>(data) : nullptr; 2259 2260 if (role == Ioss::Field::MESH) { 2261 int crd_idx = 0; 2262 if (field.get_name() == "mesh_model_coordinates_x") { 2263 CGCHECKM(cgp_coord_write(get_file_pointer(), base, zone, CGNS_ENUMV(RealDouble), "CoordinateX", 2264 &crd_idx)); 2265 CGCHECKM(cgp_coord_write_data(get_file_pointer(), base, zone, crd_idx, rmin, rmax, rdata)); 2266 } 2267 2268 else if (field.get_name() == "mesh_model_coordinates_y") { 2269 CGCHECKM(cgp_coord_write(get_file_pointer(), base, zone, CGNS_ENUMV(RealDouble), "CoordinateY", 2270 &crd_idx)); 2271 CGCHECKM(cgp_coord_write_data(get_file_pointer(), base, zone, crd_idx, rmin, rmax, rdata)); 2272 } 2273 2274 else if (field.get_name() == "mesh_model_coordinates_z") { 2275 CGCHECKM(cgp_coord_write(get_file_pointer(), base, zone, CGNS_ENUMV(RealDouble), "CoordinateZ", 2276 &crd_idx)); 2277 CGCHECKM(cgp_coord_write_data(get_file_pointer(), base, zone, crd_idx, rmin, rmax, rdata)); 2278 } 2279 2280 else if (field.get_name() == "mesh_model_coordinates") { 2281 int phys_dimension = get_region()->get_property("spatial_dimension").get_int(); 2282 2283 std::vector<double> coord(num_to_get); 2284 2285 // ======================================================================== 2286 // Repetitive code for each coordinate direction; use a lambda to consolidate... 2287 auto coord_lambda = [=, &coord](const char *ordinate, int ordinal) { 2288 // Data required by upper classes store x0, y0, z0, ... xn, 2289 // yn, zn. Data stored in cgns file is x0, ..., xn, y0, 2290 // ..., yn, z0, ..., zn so we have to allocate some scratch 2291 // memory to read in the data and then map into supplied 2292 // 'data' 2293 // Map to global coordinate position... 2294 for (cgsize_t i = 0; i < num_to_get; i++) { 2295 coord[i] = rdata[phys_dimension * i + ordinal]; 2296 } 2297 2298 int idx = 0; 2299 CGCHECKM(cgp_coord_write(get_file_pointer(), base, zone, CGNS_ENUMV(RealDouble), ordinate, &idx)); 2300 CGCHECKM( 2301 cgp_coord_write_data(get_file_pointer(), base, zone, idx, rmin, rmax, coord.data())); 2302 }; 2303 // ======================================================================== 2304 2305 coord_lambda("CoordinateX", 0); 2306 2307 if (phys_dimension >= 2) { 2308 coord_lambda("CoordinateY", 1); 2309 } 2310 2311 if (phys_dimension == 3) { 2312 coord_lambda("CoordinateZ", 2); 2313 } 2314 } 2315 else { 2316 num_to_get = Ioss::Utils::field_warning(sb, field, "output"); 2317 } 2318 } 2319 else if (role == Ioss::Field::TRANSIENT) { 2320 int cgns_field = 0; 2321 auto var_type = field.transformed_storage(); 2322 int comp_count = var_type->component_count(); 2323 char field_suffix_separator = get_field_separator(); 2324 2325 int sol_index = 0; 2326 CGNS_ENUMT(GridLocation_t) location; 2327 if (cell_field) { 2328 sol_index = m_currentCellCenterSolutionIndex; 2329 location = CGNS_ENUMV(CellCenter); 2330 } 2331 else { 2332 sol_index = m_currentVertexSolutionIndex; 2333 location = CGNS_ENUMV(Vertex); 2334 } 2335 if (comp_count == 1) { 2336 CGCHECKM(cgp_field_write(get_file_pointer(), base, zone, sol_index, CGNS_ENUMV(RealDouble), 2337 field.get_name().c_str(), &cgns_field)); 2338 Utils::set_field_index(field, cgns_field, location); 2339 2340 CGCHECKM(cgp_field_write_data(get_file_pointer(), base, zone, sol_index, cgns_field, rmin, 2341 rmax, rdata)); 2342 } 2343 else { 2344 std::vector<double> cgns_data(num_to_get); 2345 for (int i = 0; i < comp_count; i++) { 2346 for (cgsize_t j = 0; j < num_to_get; j++) { 2347 cgns_data[j] = rdata[comp_count * j + i]; 2348 } 2349 std::string var_name = 2350 var_type->label_name(field.get_name(), i + 1, field_suffix_separator); 2351 2352 CGCHECKM(cgp_field_write(get_file_pointer(), base, zone, sol_index, CGNS_ENUMV(RealDouble), 2353 var_name.c_str(), &cgns_field)); 2354 if (i == 0) { 2355 Utils::set_field_index(field, cgns_field, location); 2356 } 2357 2358 CGCHECKM(cgp_field_write_data(get_file_pointer(), base, zone, sol_index, cgns_field, rmin, 2359 rmax, cgns_data.data())); 2360 } 2361 } 2362 } 2363 else { 2364 num_to_get = Ioss::Utils::field_warning(sb, field, "output"); 2365 } 2366 return num_to_get; 2367 } 2368 put_field_internal(const Ioss::FaceBlock *,const Ioss::Field &,void *,size_t)2369 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::FaceBlock * /* nb */, 2370 const Ioss::Field & /* field */, void * /* data */, 2371 size_t /* data_size */) const 2372 { 2373 return -1; 2374 } put_field_internal(const Ioss::EdgeBlock *,const Ioss::Field &,void *,size_t)2375 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::EdgeBlock * /* nb */, 2376 const Ioss::Field & /* field */, void * /* data */, 2377 size_t /* data_size */) const 2378 { 2379 return -1; 2380 } put_field_internal(const Ioss::NodeSet *,const Ioss::Field &,void *,size_t)2381 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::NodeSet * /* ns */, 2382 const Ioss::Field & /* field */, void * /* data */, 2383 size_t /* data_size */) const 2384 { 2385 return -1; 2386 } put_field_internal(const Ioss::EdgeSet *,const Ioss::Field &,void *,size_t)2387 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::EdgeSet * /* ns */, 2388 const Ioss::Field & /* field */, void * /* data */, 2389 size_t /* data_size */) const 2390 { 2391 return -1; 2392 } put_field_internal(const Ioss::FaceSet *,const Ioss::Field &,void *,size_t)2393 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::FaceSet * /* ns */, 2394 const Ioss::Field & /* field */, void * /* data */, 2395 size_t /* data_size */) const 2396 { 2397 return -1; 2398 } put_field_internal(const Ioss::ElementSet *,const Ioss::Field &,void *,size_t)2399 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::ElementSet * /* ns */, 2400 const Ioss::Field & /* field */, void * /* data */, 2401 size_t /* data_size */) const 2402 { 2403 return -1; 2404 } put_field_internal(const Ioss::SideBlock * sb,const Ioss::Field & field,void * data,size_t data_size)2405 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::SideBlock *sb, 2406 const Ioss::Field &field, void *data, 2407 size_t data_size) const 2408 { 2409 const Ioss::EntityBlock *parent_block = sb->parent_block(); 2410 if (parent_block == nullptr) { 2411 std::ostringstream errmsg; 2412 fmt::print(errmsg, 2413 "ERROR: CGNS: SideBlock {} does not have a parent-block specified. This is " 2414 "required for CGNS output.", 2415 sb->name()); 2416 IOSS_ERROR(errmsg); 2417 } 2418 2419 int base = parent_block->get_property("base").get_int(); 2420 int zone = parent_block->get_property("zone").get_int(); 2421 ssize_t num_to_get = field.verify(data_size); 2422 2423 Ioss::Field::RoleType role = field.get_role(); 2424 2425 if (role == Ioss::Field::MESH) { 2426 // Handle the MESH fields required for a CGNS file model. 2427 // (The 'genesis' portion) 2428 if (field.get_name() == "element_side") { 2429 CGNS_ENUMT(ElementType_t) type = Utils::map_topology_to_cgns(sb->topology()->name()); 2430 int sect = 0; 2431 2432 int64_t size = num_to_get; 2433 MPI_Allreduce(MPI_IN_PLACE, &size, 1, Ioss::mpi_type(size), MPI_SUM, util().communicator()); 2434 2435 int cg_start = m_bcOffset[zone] + 1; 2436 int cg_end = m_bcOffset[zone] + size; 2437 2438 // NOTE: Currently not writing the "ElementConnectivity" data for the 2439 // boundary condition. It isn't used in the read and don't have 2440 // the data so would have to generate it. This may cause problems 2441 // with codes that use the downstream data if they base the BC off 2442 // of the nodes instead of the element/side info. 2443 // Get name from parent sideset... This is name of the ZoneBC entry 2444 auto &name = sb->owner()->name(); 2445 // This is the name of the BC_t node 2446 auto sb_name = Iocgns::Utils::decompose_sb_name(sb->name()); 2447 2448 CGNSIntVector point_range{cg_start, cg_end}; 2449 CGCHECKM(cg_boco_write(get_file_pointer(), base, zone, name.c_str(), CGNS_ENUMV(FamilySpecified), 2450 CGNS_ENUMV(PointRange), 2, point_range.data(), §)); 2451 CGCHECKM( 2452 cg_goto(get_file_pointer(), base, "Zone_t", zone, "ZoneBC_t", 1, "BC_t", sect, "end")); 2453 CGCHECKM(cg_famname_write(name.c_str())); 2454 CGCHECKM(cg_boco_gridlocation_write(get_file_pointer(), base, zone, sect, CGNS_ENUMV(FaceCenter))); 2455 2456 CGCHECKM(cgp_section_write(get_file_pointer(), base, zone, sb_name.c_str(), type, cg_start, 2457 cg_end, 0, §)); 2458 2459 sb->property_update("section", sect); 2460 2461 CGNSIntVector parent(4 * num_to_get); 2462 2463 if (field.get_type() == Ioss::Field::INT32) { 2464 int * idata = reinterpret_cast<int *>(data); 2465 size_t j = 0; 2466 for (ssize_t i = 0; i < num_to_get; i++) { 2467 parent[num_to_get * 0 + i] = elemMap.global_to_local(idata[j++]); // Element 2468 parent[num_to_get * 2 + i] = idata[j++]; 2469 } 2470 } 2471 else { 2472 auto * idata = reinterpret_cast<int64_t *>(data); 2473 size_t j = 0; 2474 for (ssize_t i = 0; i < num_to_get; i++) { 2475 parent[num_to_get * 0 + i] = elemMap.global_to_local(idata[j++]); // Element 2476 parent[num_to_get * 2 + i] = idata[j++]; 2477 } 2478 } 2479 // Adjust face numbers to CGNS convention instead of IOSS convention... 2480 Utils::map_ioss_face_to_cgns(sb->parent_element_topology(), num_to_get, parent); 2481 map_local_to_global_implicit(parent, num_to_get, m_elemGlobalImplicitMap); 2482 2483 int64_t local_face_count = num_to_get; 2484 int64_t local_face_offset = 0; 2485 MPI_Exscan(&local_face_count, &local_face_offset, 1, Ioss::mpi_type(local_face_count), 2486 MPI_SUM, util().communicator()); 2487 cg_start = m_bcOffset[zone] + local_face_offset + 1; 2488 cg_end = cg_start + local_face_count - 1; 2489 2490 auto xx = num_to_get > 0 ? parent.data() : nullptr; 2491 if (num_to_get == 0) { 2492 cg_start = cg_end = 0; 2493 } 2494 CGCHECKM(cgp_parent_data_write(get_file_pointer(), base, zone, sect, cg_start, cg_end, xx)); 2495 m_bcOffset[zone] += size; 2496 } 2497 else if (field.get_name() == "distribution_factors") { 2498 static bool warning_output = false; 2499 if (!warning_output) { 2500 if (myProcessor == 0) { 2501 fmt::print(Ioss::WARNING(), 2502 "For CGNS output, the sideset distribution factors are not output.\n"); 2503 } 2504 warning_output = true; 2505 } 2506 return 0; 2507 } 2508 else { 2509 num_to_get = Ioss::Utils::field_warning(sb, field, "output"); 2510 } 2511 } 2512 else { 2513 num_to_get = Ioss::Utils::field_warning(sb, field, "output"); 2514 } 2515 return num_to_get; 2516 } put_field_internal(const Ioss::SideSet *,const Ioss::Field &,void *,size_t)2517 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::SideSet * /* fs */, 2518 const Ioss::Field & /* field */, void * /* data */, 2519 size_t /* data_size */) const 2520 { 2521 return -1; 2522 } put_field_internal(const Ioss::CommSet *,const Ioss::Field &,void *,size_t)2523 int64_t ParallelDatabaseIO::put_field_internal(const Ioss::CommSet * /* cs */, 2524 const Ioss::Field & /* field*/, void * /*data*/, 2525 size_t /*data_size*/) const 2526 { 2527 return -1; 2528 } 2529 write_results_meta_data()2530 void ParallelDatabaseIO::write_results_meta_data() {} 2531 entity_field_support()2532 unsigned ParallelDatabaseIO::entity_field_support() const 2533 { 2534 return Ioss::NODEBLOCK | Ioss::ELEMENTBLOCK | Ioss::STRUCTUREDBLOCK | Ioss::NODESET | 2535 Ioss::SIDESET | Ioss::REGION; 2536 } 2537 handle_node_ids(void * ids,int64_t num_to_get)2538 int64_t ParallelDatabaseIO::handle_node_ids(void *ids, int64_t num_to_get) const 2539 { 2540 /*! 2541 * There are two modes we need to support in this routine: 2542 * 1. Initial definition of node map (local->global) and 2543 * nodeMap.reverse (global->local). 2544 * 2. Redefinition of node map via 'reordering' of the original 2545 * map when the nodes on this processor are the same, but their 2546 * order is changed (or count because of ghosting) 2547 * 2548 * So, there will be two maps the 'nodeMap.map' map is a 'direct lookup' 2549 * map which maps current local position to global id and the 2550 * 'nodeMap.reverse' is an associative lookup which maps the 2551 * global id to 'original local'. There is also a 2552 * 'nodeMap.reorder' which is direct lookup and maps current local 2553 * position to original local. 2554 2555 * The ids coming in are the global ids; their position is the 2556 * "local id-1" (That is, data[0] contains the global id of local 2557 * node 1 in this node block). 2558 * 2559 * int local_position = nodeMap.reverse[NodeMap[i+1]] 2560 * (the nodeMap.map and nodeMap.reverse are 1-based) 2561 * 2562 * To determine which map to update on a call to this function, we 2563 * use the following hueristics: 2564 * -- If the database state is 'STATE_MODEL:', then update the 2565 * 'nodeMap.reverse' and 'nodeMap.map' 2566 * 2567 * -- If the database state is not STATE_MODEL, then leave the 2568 * 'nodeMap.reverse' and 'nodeMap.map' alone since they correspond to the 2569 * information already written to the database. [May want to add a 2570 * STATE_REDEFINE_MODEL] 2571 * 2572 * -- In both cases, update the nodeMap.reorder 2573 * 2574 * NOTE: The mapping is done on TRANSIENT fields only; MODEL fields 2575 * should be in the original order... 2576 */ 2577 if (!nodeMap.defined()) { 2578 nodeMap.set_size(num_to_get); 2579 2580 bool in_define = (dbState == Ioss::STATE_MODEL) || (dbState == Ioss::STATE_DEFINE_MODEL); 2581 if (nodeMap.is_sequential()) { 2582 if (int_byte_size_api() == 4) { 2583 nodeMap.set_map(static_cast<int *>(ids), num_to_get, 0, in_define); 2584 } 2585 else { 2586 nodeMap.set_map(static_cast<int64_t *>(ids), num_to_get, 0, in_define); 2587 } 2588 } 2589 2590 // Only a single nodeblock and all set 2591 assert(get_region()->get_property("node_block_count").get_int() == 1); 2592 nodeMap.set_defined(true); 2593 } 2594 return num_to_get; 2595 } 2596 get_processor_zone_node_offset()2597 std::vector<int64_t> ParallelDatabaseIO::get_processor_zone_node_offset() const 2598 { 2599 size_t num_zones = m_globalToBlockLocalNodeMap.size(); 2600 std::vector<int64_t> node_count(num_zones); 2601 std::vector<int64_t> node_offset(num_zones); 2602 2603 for (const auto &block : m_globalToBlockLocalNodeMap) { 2604 auto zone = block.first; 2605 const auto &block_map = block.second; 2606 node_count[zone - 1] = block_map->size(); 2607 } 2608 MPI_Exscan(node_count.data(), node_offset.data(), num_zones, Ioss::mpi_type(node_count[0]), 2609 MPI_SUM, util().communicator()); 2610 2611 return node_offset; 2612 } 2613 2614 } // namespace Iocgns 2615