| 77 | namespace seissol::initializer { |
| 78 | |
| 79 | void projectInitialField(const std::vector<std::unique_ptr<physics::InitialField>>& iniFields, |
| 80 | const GlobalData& globalData, |
| 81 | const seissol::geometry::MeshReader& meshReader, |
| 82 | seissol::initializer::MemoryManager& memoryManager, |
| 83 | LTS const& lts, |
| 84 | const Lut& ltsLut) { |
| 85 | const auto& vertices = meshReader.getVertices(); |
| 86 | const auto& elements = meshReader.getElements(); |
| 87 | |
| 88 | constexpr auto QuadPolyDegree = ConvergenceOrder + 1; |
| 89 | constexpr auto NumQuadPoints = QuadPolyDegree * QuadPolyDegree * QuadPolyDegree; |
| 90 | |
| 91 | double quadraturePoints[NumQuadPoints][3]; |
| 92 | double quadratureWeights[NumQuadPoints]; |
| 93 | seissol::quadrature::TetrahedronQuadrature(quadraturePoints, quadratureWeights, QuadPolyDegree); |
| 94 | |
| 95 | #if defined(_OPENMP) && !NVHPC_AVOID_OMP |
| 96 | #pragma omp parallel |
| 97 | { |
| 98 | #endif |
| 99 | alignas(Alignment) real iniCondData[tensor::iniCond::size()] = {}; |
| 100 | auto iniCond = init::iniCond::view::create(iniCondData); |
| 101 | |
| 102 | std::vector<std::array<double, 3>> quadraturePointsXyz; |
| 103 | quadraturePointsXyz.resize(NumQuadPoints); |
| 104 | |
| 105 | kernel::projectIniCond krnl; |
| 106 | krnl.projectQP = globalData.projectQPMatrix; |
| 107 | krnl.iniCond = iniCondData; |
| 108 | kernels::set_selectAneFull(krnl, kernels::get_static_ptr_Values<init::selectAneFull>()); |
| 109 | kernels::set_selectElaFull(krnl, kernels::get_static_ptr_Values<init::selectElaFull>()); |
| 110 | |
| 111 | #if defined(_OPENMP) && !NVHPC_AVOID_OMP |
| 112 | #pragma omp for schedule(static) |
| 113 | #endif |
| 114 | for (unsigned int meshId = 0; meshId < elements.size(); ++meshId) { |
| 115 | const double* elementCoords[4]; |
| 116 | for (size_t v = 0; v < 4; ++v) { |
| 117 | elementCoords[v] = vertices[elements[meshId].vertices[v]].coords; |
| 118 | } |
| 119 | for (size_t i = 0; i < NumQuadPoints; ++i) { |
| 120 | seissol::transformations::tetrahedronReferenceToGlobal(elementCoords[0], |
| 121 | elementCoords[1], |
| 122 | elementCoords[2], |
| 123 | elementCoords[3], |
| 124 | quadraturePoints[i], |
| 125 | quadraturePointsXyz[i].data()); |
| 126 | } |
| 127 | |
| 128 | const CellMaterialData& material = ltsLut.lookup(lts.material, meshId); |
| 129 | #ifdef MULTIPLE_SIMULATIONS |
| 130 | for (int s = 0; s < MULTIPLE_SIMULATIONS; ++s) { |
| 131 | auto sub = iniCond.subtensor(s, yateto::slice<>(), yateto::slice<>()); |
| 132 | iniFields[s % iniFields.size()]->evaluate(0.0, quadraturePointsXyz, material, sub); |
| 133 | } |
| 134 | #else |
| 135 | iniFields[0]->evaluate(0.0, quadraturePointsXyz, material, iniCond); |
| 136 | #endif |
no test coverage detected