Function loadSourcesFromFSRM

src/SourceTerm/Manager.cpp:356–463

}

auto loadSourcesFromFSRM(const char* fileName,
                         const seissol::geometry::MeshReader& mesh,
                         seissol::initializer::LTSTree* ltsTree,
                         seissol::initializer::LTS* lts,
                         seissol::initializer::Lut* ltsLut,
                         seissol::memory::Memkind memkind)
    -> std::unordered_map<LayerType, std::vector<seissol::kernels::PointSourceClusterPair>> {
  // until further rewrite, we'll leave most of the raw pointers/arrays in here.

  const int rank = seissol::MPI::mpi.rank();

  seissol::sourceterm::FSRMSource fsrm;
  fsrm.read(std::string(fileName));

  logInfo(rank) << "Finding meshIds for point sources...";

  auto contained = std::vector<short>(fsrm.numberOfSources);
  auto meshIds = std::vector<unsigned>(fsrm.numberOfSources);

  initializer::findMeshIds(
      fsrm.centers.data(), mesh, fsrm.numberOfSources, contained.data(), meshIds.data());

#ifdef USE_MPI
  logInfo(rank) << "Cleaning possible double occurring point sources for MPI...";
  initializer::cleanDoubles(contained.data(), fsrm.numberOfSources);
#endif

  auto originalIndex = std::vector<unsigned>(fsrm.numberOfSources);
  unsigned numSources = 0;
  for (unsigned source = 0; source < fsrm.numberOfSources; ++source) {
    originalIndex[numSources] = source;
    meshIds[numSources] = meshIds[source];
    numSources += contained[source];
  }

  logInfo(rank) << "Mapping point sources to LTS cells...";
  auto layeredClusterMapping =
      mapPointSourcesToClusters(meshIds.data(), numSources, ltsTree, lts, ltsLut, memkind);
  std::unordered_map<LayerType, std::vector<seissol::kernels::PointSourceClusterPair>>
      layeredSourceClusters;

  for (auto layer : {Interior, Copy}) {
    auto& sourceCluster = layeredSourceClusters[layer];
    sourceCluster.resize(ltsTree->numChildren());
    auto& clusterMappings = layeredClusterMapping[layer];
    for (unsigned cluster = 0; cluster < ltsTree->numChildren(); ++cluster) {
      auto numberOfSources = clusterMappings[cluster].sources.size();
      auto sources = PointSources{memkind};
      sources.mode = PointSourceMode::Fsrm;
      sources.numberOfSources = numberOfSources;
      sources.mInvJInvPhisAtSources.resize(numberOfSources);
      sources.tensor.resize(numberOfSources);
      sources.onsetTime.resize(numberOfSources);
      sources.samplingInterval.resize(numberOfSources);
      sources.sampleOffsets[0].resize(numberOfSources + 1);
      sources.sampleOffsets[0][0] = 0;
      sources.sample[0].resize(fsrm.numberOfSamples * numberOfSources);