| 465 | // TODO Add support for passive netCDF |
| 466 | #if defined(USE_NETCDF) && !defined(NETCDF_PASSIVE) |
| 467 | auto loadSourcesFromNRF(const char* fileName, |
| 468 | const seissol::geometry::MeshReader& mesh, |
| 469 | seissol::initializer::LTSTree* ltsTree, |
| 470 | seissol::initializer::LTS* lts, |
| 471 | seissol::initializer::Lut* ltsLut, |
| 472 | seissol::memory::Memkind memkind) |
| 473 | -> std::unordered_map<LayerType, std::vector<seissol::kernels::PointSourceClusterPair>> { |
| 474 | const int rank = seissol::MPI::mpi.rank(); |
| 475 | |
| 476 | logInfo(rank) << "Reading" << fileName; |
| 477 | NRF nrf; |
| 478 | readNRF(fileName, nrf); |
| 479 | |
| 480 | auto contained = std::vector<short>(nrf.size()); |
| 481 | auto meshIds = std::vector<unsigned>(nrf.size()); |
| 482 | |
| 483 | logInfo(rank) << "Finding meshIds for point sources..."; |
| 484 | initializer::findMeshIds(nrf.centres.data(), mesh, nrf.size(), contained.data(), meshIds.data()); |
| 485 | |
| 486 | #ifdef USE_MPI |
| 487 | logInfo(rank) << "Cleaning possible double occurring point sources for MPI..."; |
| 488 | initializer::cleanDoubles(contained.data(), nrf.size()); |
| 489 | #endif |
| 490 | |
| 491 | auto originalIndex = std::vector<unsigned>(nrf.size()); |
| 492 | unsigned numSources = 0; |
| 493 | for (unsigned source = 0; source < nrf.size(); ++source) { |
| 494 | originalIndex[numSources] = source; |
| 495 | meshIds[numSources] = meshIds[source]; |
| 496 | numSources += contained[source]; |
| 497 | } |
| 498 | |
| 499 | // Checking that all sources are within the domain |
| 500 | unsigned globalnumSources = numSources; |
| 501 | #ifdef USE_MPI |
| 502 | MPI_Reduce(&numSources, &globalnumSources, 1, MPI_UNSIGNED, MPI_SUM, 0, seissol::MPI::mpi.comm()); |
| 503 | #endif |
| 504 | |
| 505 | if (rank == 0) { |
| 506 | const int numSourceOutside = nrf.size() - globalnumSources; |
| 507 | if (numSourceOutside > 0) { |
| 508 | logError() << nrf.size() - globalnumSources << " point sources are outside the domain."; |
| 509 | } |
| 510 | } |
| 511 | |
| 512 | logInfo(rank) << "Mapping point sources to LTS cells..."; |
| 513 | |
| 514 | auto layeredClusterMapping = |
| 515 | mapPointSourcesToClusters(meshIds.data(), numSources, ltsTree, lts, ltsLut, memkind); |
| 516 | std::unordered_map<LayerType, std::vector<seissol::kernels::PointSourceClusterPair>> |
| 517 | layeredSourceClusters; |
| 518 | |
| 519 | for (auto layer : {Interior, Copy}) { |
| 520 | auto& sourceCluster = layeredSourceClusters[layer]; |
| 521 | sourceCluster.resize(ltsTree->numChildren()); |
| 522 | auto& clusterMappings = layeredClusterMapping[layer]; |
| 523 | for (unsigned cluster = 0; cluster < ltsTree->numChildren(); ++cluster) { |
| 524 | auto numberOfSources = clusterMappings[cluster].sources.size(); |
no test coverage detected