| 24 | |
| 25 | namespace seissol::dr::initializer { |
| 26 | void BaseDRInitializer::initializeFault(const seissol::initializer::DynamicRupture* const dynRup, |
| 27 | seissol::initializer::LTSTree* const dynRupTree) { |
| 28 | const int rank = seissol::MPI::mpi.rank(); |
| 29 | logInfo(rank) << "Initializing Fault, using a quadrature rule with " |
| 30 | << misc::NumBoundaryGaussPoints << " points."; |
| 31 | seissol::initializer::FaultParameterDB faultParameterDB; |
| 32 | for (auto& layer : dynRupTree->leaves(Ghost)) { |
| 33 | // parameters to be read from fault parameters yaml file |
| 34 | std::unordered_map<std::string, real*> parameterToStorageMap; |
| 35 | |
| 36 | // read initial stress and nucleation stress |
| 37 | auto addStressesToStorageMap = [¶meterToStorageMap, &layer, this]( |
| 38 | StressTensor& initialStress, bool readNucleation) { |
| 39 | // return pointer to first element |
| 40 | auto getRawData = [](StressTensor::VectorOfArraysT& vectorOfArrays) { |
| 41 | return vectorOfArrays.data()->data(); |
| 42 | }; |
| 43 | // fault can be either initialized by traction or by cartesian stress |
| 44 | // this method reads either the nucleation stress or the initial stress |
| 45 | auto [identifiers, parametrization] = this->stressIdentifiers(readNucleation); |
| 46 | const bool isFaultParameterizedByTraction = parametrization == Parametrization::Traction; |
| 47 | if (isFaultParameterizedByTraction) { |
| 48 | // only read traction in normal, strike and dip direction |
| 49 | parameterToStorageMap.insert({identifiers[0], getRawData(initialStress.xx)}); |
| 50 | parameterToStorageMap.insert({identifiers[1], getRawData(initialStress.xy)}); |
| 51 | parameterToStorageMap.insert({identifiers[2], getRawData(initialStress.xz)}); |
| 52 | // set the rest to zero |
| 53 | for (unsigned ltsFace = 0; ltsFace < layer.getNumberOfCells(); ++ltsFace) { |
| 54 | for (unsigned pointIndex = 0; pointIndex < init::QInterpolated::Stop[0]; ++pointIndex) { |
| 55 | initialStress.yy[ltsFace][pointIndex] = 0.0; |
| 56 | initialStress.zz[ltsFace][pointIndex] = 0.0; |
| 57 | initialStress.yz[ltsFace][pointIndex] = 0.0; |
| 58 | } |
| 59 | } |
| 60 | } else { // read all stress components from the parameter file |
| 61 | parameterToStorageMap.insert({identifiers[0], getRawData(initialStress.xx)}); |
| 62 | parameterToStorageMap.insert({identifiers[1], getRawData(initialStress.yy)}); |
| 63 | parameterToStorageMap.insert({identifiers[2], getRawData(initialStress.zz)}); |
| 64 | parameterToStorageMap.insert({identifiers[3], getRawData(initialStress.xy)}); |
| 65 | parameterToStorageMap.insert({identifiers[4], getRawData(initialStress.yz)}); |
| 66 | parameterToStorageMap.insert({identifiers[5], getRawData(initialStress.xz)}); |
| 67 | } |
| 68 | #ifdef USE_POROELASTIC |
| 69 | if (isFaultParameterizedByTraction) { |
| 70 | parameterToStorageMap.insert({identifiers[3], getRawData(initialStress.p)}); |
| 71 | } else { |
| 72 | parameterToStorageMap.insert({identifiers[6], getRawData(initialStress.p)}); |
| 73 | } |
| 74 | #else |
| 75 | for (unsigned ltsFace = 0; ltsFace < layer.getNumberOfCells(); ++ltsFace) { |
| 76 | for (unsigned pointIndex = 0; pointIndex < init::QInterpolated::Stop[0]; ++pointIndex) { |
| 77 | initialStress.p[ltsFace][pointIndex] = 0.0; |
| 78 | } |
| 79 | } |
| 80 | #endif |
| 81 | |
| 82 | return isFaultParameterizedByTraction; |
| 83 | }; |
no test coverage detected