diff --git a/opm/simulators/flow/Transmissibility_impl.hpp b/opm/simulators/flow/Transmissibility_impl.hpp index 60c1925c6..67141cc6e 100644 --- a/opm/simulators/flow/Transmissibility_impl.hpp +++ b/opm/simulators/flow/Transmissibility_impl.hpp @@ -23,21 +23,27 @@ #ifndef OPM_TRANSMISSIBILITY_IMPL_HPP #define OPM_TRANSMISSIBILITY_IMPL_HPP +#ifndef OPM_TRANSMISSIBILITY_HPP +#include +#include +#endif + #include #include -#include - #include +#include + +#include +#include + #include #include #include #include #include -#include - -#include +#include #include #include @@ -53,6 +59,8 @@ #include #include +#include + namespace Opm { namespace details { @@ -184,19 +192,25 @@ update(bool global, const TransUpdateQuantities update_quantities, extractPermeability_(); } + const int num_threads = ThreadManager::maxThreads(); + // reserving some space in the hashmap upfront saves quite a bit of time because // resizes are costly for hashmaps and there would be quite a few of them if we // would not have a rough idea of how large the final map will be (the rough idea // is a conforming Cartesian grid). trans_.clear(); - trans_.reserve(numElements*3*1.05); + if (num_threads == 1) { + trans_.reserve(numElements*3*1.05); + } transBoundary_.clear(); // if energy is enabled, let's do the same for the "thermal half transmissibilities" if (enableEnergy_ && !onlyTrans) { thermalHalfTrans_.clear(); - thermalHalfTrans_.reserve(numElements*6*1.05); + if (num_threads == 1) { + thermalHalfTrans_.reserve(numElements*6*1.05); + } thermalHalfTransBoundary_.clear(); } @@ -204,14 +218,18 @@ update(bool global, const TransUpdateQuantities update_quantities, // if diffusion is enabled, let's do the same for the "diffusivity" if (updateDiffusivity && !onlyTrans) { diffusivity_.clear(); - diffusivity_.reserve(numElements*3*1.05); + if (num_threads == 1) { + diffusivity_.reserve(numElements*3*1.05); + } extractPorosity_(); } // if dispersion is enabled, let's do the same for the "dispersivity" if (updateDispersivity && !onlyTrans) { dispersivity_.clear(); - dispersivity_.reserve(numElements*3*1.05); + if (num_threads == 1) { + dispersivity_.reserve(numElements*3*1.05); + } extractDispersion_(); } @@ -279,206 +297,252 @@ update(bool global, const TransUpdateQuantities update_quantities, prop); }; - // compute the transmissibilities for all intersections - for (const auto& elem : elements(gridView_)) { - FaceInfo inside; - FaceInfo outside; - DimVector faceAreaNormal; + ThreadSafeMapBuilder transBoundary(transBoundary_, num_threads, + MapBuilderInsertionMode::Insert_Or_Assign); + ThreadSafeMapBuilder transMap(trans_, num_threads, + MapBuilderInsertionMode::Insert_Or_Assign); + ThreadSafeMapBuilder thermalHalfTransBoundary(thermalHalfTransBoundary_, num_threads, + MapBuilderInsertionMode::Insert_Or_Assign); + ThreadSafeMapBuilder thermalHalfTrans(thermalHalfTrans_, num_threads, + MapBuilderInsertionMode::Insert_Or_Assign); + ThreadSafeMapBuilder diffusivity(diffusivity_, num_threads, + MapBuilderInsertionMode::Insert_Or_Assign); + ThreadSafeMapBuilder dispersivity(dispersivity_, num_threads, + MapBuilderInsertionMode::Insert_Or_Assign); - inside.elemIdx = elemMapper.index(elem); - // Get the Cartesian index of the origin cells (parent or equivalent cell on level zero), - // for CpGrid with LGRs. For general grids and no LGRs, get the usual Cartesian Index. - inside.cartElemIdx = this->lookUpCartesianData_. - template getFieldPropCartesianIdx(inside.elemIdx); - - auto computeHalf = [this, &faceAreaNormal, &inside, &outside] - (const auto& halfComputer, - const auto& prop1, const auto& prop2) -> std::array - { - return { - halfComputer(faceAreaNormal, - inside.faceIdx, - distanceVector_(inside.faceCenter, inside.elemIdx), - prop1), - halfComputer(faceAreaNormal, - outside.faceIdx, - distanceVector_(outside.faceCenter, outside.elemIdx), - prop2) - }; - }; - - auto computeHalfMean = [&inside, &outside, &computeHalf, &ntg, &harmonicMean] - (const auto& halfComputer, const auto& prop) - { - auto half = computeHalf(halfComputer, prop[inside.elemIdx], prop[outside.elemIdx]); - applyNtg_(half[0], inside, ntg); - applyNtg_(half[1], outside, ntg); - - //TODO Add support for multipliers - return harmonicMean(half[0], half[1]); - }; - - unsigned boundaryIsIdx = 0; - for (const auto& intersection : intersections(gridView_, elem)) { - // deal with grid boundaries - if (intersection.boundary()) { - // compute the transmissibilty for the boundary intersection - const auto& geometry = intersection.geometry(); - inside.faceCenter = geometry.center(); - - faceAreaNormal = intersection.centerUnitOuterNormal(); - faceAreaNormal *= geometry.volume(); - - Scalar transBoundaryIs = - computeHalfTrans_(faceAreaNormal, - intersection.indexInInside(), - distanceVector_(inside.faceCenter, inside.elemIdx), - permeability_[inside.elemIdx]); - - // normally there would be two half-transmissibilities that would be - // averaged. on the grid boundary there only is the half - // transmissibility of the interior element. - applyMultipliers_(transBoundaryIs, intersection.indexInInside(), inside.cartElemIdx, transMult); - transBoundary_.insert_or_assign(std::make_pair(inside.elemIdx, boundaryIsIdx), transBoundaryIs); - - // for boundary intersections we also need to compute the thermal - // half transmissibilities - if (enableEnergy_ && !onlyTrans) { - Scalar transBoundaryEnergyIs = - computeHalfDiffusivity_(faceAreaNormal, - distanceVector_(inside.faceCenter, inside.elemIdx), - 1.0); - thermalHalfTransBoundary_.insert_or_assign(std::make_pair(inside.elemIdx, boundaryIsIdx), - transBoundaryEnergyIs); - } - - ++boundaryIsIdx; - continue; - } - - if (!intersection.neighbor()) { - // elements can be on process boundaries, i.e. they are not on the - // domain boundary yet they don't have neighbors. - ++boundaryIsIdx; - continue; - } - - const auto& outsideElem = intersection.outside(); - outside.elemIdx = elemMapper.index(outsideElem); +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (const auto& chunk : ElementChunks(gridView_, num_threads)) { + for (const auto& elem : chunk) { + FaceInfo inside; + FaceInfo outside; + DimVector faceAreaNormal; + inside.elemIdx = elemMapper.index(elem); // Get the Cartesian index of the origin cells (parent or equivalent cell on level zero), // for CpGrid with LGRs. For general grids and no LGRs, get the usual Cartesian Index. - outside.cartElemIdx = this->lookUpCartesianData_. - template getFieldPropCartesianIdx(outside.elemIdx); + inside.cartElemIdx = this->lookUpCartesianData_. + template getFieldPropCartesianIdx(inside.elemIdx); - // we only need to calculate a face's transmissibility - // once... - // In a parallel run inside.cartElemIdx > outside.cartElemIdx does not imply inside.elemIdx > outside.elemIdx for - // ghost cells and we need to use the cartesian index as this will be used when applying Z multipliers - // To cover the case where both cells are part of an LGR and as a consequence might have - // the same cartesian index, we tie their Cartesian indices and the ones on the leaf grid view. - if (std::tie(inside.cartElemIdx, inside.elemIdx) > std::tie(outside.cartElemIdx, outside.elemIdx)) { - continue; - } + auto computeHalf = [this, &faceAreaNormal, &inside, &outside] + (const auto& halfComputer, + const auto& prop1, const auto& prop2) -> std::array + { + return { + halfComputer(faceAreaNormal, + inside.faceIdx, + distanceVector_(inside.faceCenter, inside.elemIdx), + prop1), + halfComputer(faceAreaNormal, + outside.faceIdx, + distanceVector_(outside.faceCenter, outside.elemIdx), + prop2) + }; + }; - // local indices of the faces of the inside and - // outside elements which contain the intersection - inside.faceIdx = intersection.indexInInside(); - outside.faceIdx = intersection.indexInOutside(); + auto computeHalfMean = [&inside, &outside, &computeHalf, &ntg, &harmonicMean] + (const auto& halfComputer, const auto& prop) + { + auto half = computeHalf(halfComputer, prop[inside.elemIdx], prop[outside.elemIdx]); + applyNtg_(half[0], inside, ntg); + applyNtg_(half[1], outside, ntg); - if (inside.faceIdx == -1) { - // NNC. Set zero transmissibility, as it will be - // *added to* by applyNncToGridTrans_() later. - assert(outside.faceIdx == -1); - trans_.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), 0.0); - if (enableEnergy_ && !onlyTrans) { - thermalHalfTrans_.insert_or_assign(details::directionalIsId(inside.elemIdx, outside.elemIdx), 0.0); - thermalHalfTrans_.insert_or_assign(details::directionalIsId(outside.elemIdx, inside.elemIdx), 0.0); + //TODO Add support for multipliers + return harmonicMean(half[0], half[1]); + }; + + unsigned boundaryIsIdx = 0; + for (const auto& intersection : intersections(gridView_, elem)) { + // deal with grid boundaries + if (intersection.boundary()) { + // compute the transmissibilty for the boundary intersection + const auto& geometry = intersection.geometry(); + inside.faceCenter = geometry.center(); + + faceAreaNormal = intersection.centerUnitOuterNormal(); + faceAreaNormal *= geometry.volume(); + + Scalar transBoundaryIs = + computeHalfTrans_(faceAreaNormal, + intersection.indexInInside(), + distanceVector_(inside.faceCenter, inside.elemIdx), + permeability_[inside.elemIdx]); + + // normally there would be two half-transmissibilities that would be + // averaged. on the grid boundary there only is the half + // transmissibility of the interior element. + applyMultipliers_(transBoundaryIs, intersection.indexInInside(), inside.cartElemIdx, transMult); + transBoundary.insert_or_assign(std::make_pair(inside.elemIdx, boundaryIsIdx), transBoundaryIs); + + // for boundary intersections we also need to compute the thermal + // half transmissibilities + if (enableEnergy_ && !onlyTrans) { + Scalar transBoundaryEnergyIs = + computeHalfDiffusivity_(faceAreaNormal, + distanceVector_(inside.faceCenter, inside.elemIdx), + 1.0); + thermalHalfTransBoundary.insert_or_assign(std::make_pair(inside.elemIdx, boundaryIsIdx), + transBoundaryEnergyIs); + } + + ++boundaryIsIdx; + continue; } - if (updateDiffusivity && !onlyTrans) { - diffusivity_.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), 0.0); + if (!intersection.neighbor()) { + // elements can be on process boundaries, i.e. they are not on the + // domain boundary yet they don't have neighbors. + ++boundaryIsIdx; + continue; } - if (updateDispersivity && !onlyTrans) { - dispersivity_.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), 0.0); + + const auto& outsideElem = intersection.outside(); + outside.elemIdx = elemMapper.index(outsideElem); + + // Get the Cartesian index of the origin cells (parent or equivalent cell on level zero), + // for CpGrid with LGRs. For general grids and no LGRs, get the usual Cartesian Index. + outside.cartElemIdx = this->lookUpCartesianData_. + template getFieldPropCartesianIdx(outside.elemIdx); + + // we only need to calculate a face's transmissibility + // once... + // In a parallel run inside.cartElemIdx > outside.cartElemIdx does not imply inside.elemIdx > outside.elemIdx for + // ghost cells and we need to use the cartesian index as this will be used when applying Z multipliers + // To cover the case where both cells are part of an LGR and as a consequence might have + // the same cartesian index, we tie their Cartesian indices and the ones on the leaf grid view. + if (std::tie(inside.cartElemIdx, inside.elemIdx) > std::tie(outside.cartElemIdx, outside.elemIdx)) { + continue; } - continue; - } - typename std::is_same::type isCpGrid; - computeFaceProperties(intersection, - inside, - outside, - faceAreaNormal, - isCpGrid); + // local indices of the faces of the inside and + // outside elements which contain the intersection + inside.faceIdx = intersection.indexInInside(); + outside.faceIdx = intersection.indexInOutside(); - Scalar trans = computeHalfMean(computeHalfTrans_, permeability_); + if (inside.faceIdx == -1) { + // NNC. Set zero transmissibility, as it will be + // *added to* by applyNncToGridTrans_() later. + assert(outside.faceIdx == -1); + transMap.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), 0.0); + if (enableEnergy_ && !onlyTrans) { + thermalHalfTrans.insert_or_assign(details::directionalIsId(inside.elemIdx, outside.elemIdx), 0.0); + thermalHalfTrans.insert_or_assign(details::directionalIsId(outside.elemIdx, inside.elemIdx), 0.0); + } - // apply the full face transmissibility multipliers - // for the inside ... - if (!pinchActive) { - if (inside.faceIdx > 3) { // top or bottom - auto find_layer = [&cartDims](std::size_t cell) { - cell /= cartDims[0]; - auto k = cell / cartDims[1]; - return k; - }; - int kup = find_layer(inside.cartElemIdx); - int kdown = find_layer(outside.cartElemIdx); - // When a grid is a CpGrid with LGRs, insideCartElemIdx coincides with outsideCartElemIdx - // for cells on the leaf with the same parent cell on level zero. - assert((kup != kdown) || (inside.cartElemIdx == outside.cartElemIdx)); - if (std::abs(kup -kdown) > 1) { - trans = 0.0; + if (updateDiffusivity && !onlyTrans) { + diffusivity.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), 0.0); + } + if (updateDispersivity && !onlyTrans) { + dispersivity.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), 0.0); + } + continue; + } + + typename std::is_same::type isCpGrid; + computeFaceProperties(intersection, + inside, + outside, + faceAreaNormal, + isCpGrid); + + Scalar trans = computeHalfMean(computeHalfTrans_, permeability_); + + // apply the full face transmissibility multipliers + // for the inside ... + if (!pinchActive) { + if (inside.faceIdx > 3) { // top or bottom + auto find_layer = [&cartDims](std::size_t cell) { + cell /= cartDims[0]; + auto k = cell / cartDims[1]; + return k; + }; + int kup = find_layer(inside.cartElemIdx); + int kdown = find_layer(outside.cartElemIdx); + // When a grid is a CpGrid with LGRs, insideCartElemIdx coincides with outsideCartElemIdx + // for cells on the leaf with the same parent cell on level zero. + assert((kup != kdown) || (inside.cartElemIdx == outside.cartElemIdx)); + if (std::abs(kup -kdown) > 1) { + trans = 0.0; + } } } - } - if (useSmallestMultiplier) { - // PINCH(4) == TOPBOT is assumed here as we set useSmallestMultipliers - // to false if PINCH(4) == ALL holds - // In contrast to the name this will also apply - applyAllZMultipliers_(trans, inside, outside, transMult, cartDims); - } - else { - applyMultipliers_(trans, inside.faceIdx, inside.cartElemIdx, transMult); - // ... and outside elements - applyMultipliers_(trans, outside.faceIdx, outside.cartElemIdx, transMult); - } + if (useSmallestMultiplier) { + // PINCH(4) == TOPBOT is assumed here as we set useSmallestMultipliers + // to false if PINCH(4) == ALL holds + // In contrast to the name this will also apply + applyAllZMultipliers_(trans, inside, outside, transMult, cartDims); + } + else { + applyMultipliers_(trans, inside.faceIdx, inside.cartElemIdx, transMult); + // ... and outside elements + applyMultipliers_(trans, outside.faceIdx, outside.cartElemIdx, transMult); + } - // apply the region multipliers (cf. the MULTREGT keyword) - trans *= transMult.getRegionMultiplier(inside.cartElemIdx, - outside.cartElemIdx, - faceIdToDir(inside.faceIdx)); + // apply the region multipliers (cf. the MULTREGT keyword) + trans *= transMult.getRegionMultiplier(inside.cartElemIdx, + outside.cartElemIdx, + faceIdToDir(inside.faceIdx)); - trans_.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), trans); + transMap.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), trans); - // update the "thermal half transmissibility" for the intersection - if (enableEnergy_ && !onlyTrans) { - const auto half = computeHalf(halfDiff, 1.0, 1.0); - // TODO Add support for multipliers - thermalHalfTrans_.insert_or_assign(details::directionalIsId(inside.elemIdx, outside.elemIdx), - half[0]); - thermalHalfTrans_.insert_or_assign(details::directionalIsId(outside.elemIdx, inside.elemIdx), - half[1]); - } + // update the "thermal half transmissibility" for the intersection + if (enableEnergy_ && !onlyTrans) { + const auto half = computeHalf(halfDiff, 1.0, 1.0); + // TODO Add support for multipliers + thermalHalfTrans.insert_or_assign(details::directionalIsId(inside.elemIdx, outside.elemIdx), + half[0]); + thermalHalfTrans.insert_or_assign(details::directionalIsId(outside.elemIdx, inside.elemIdx), + half[1]); + } - // update the "diffusive half transmissibility" for the intersection - if (updateDiffusivity && !onlyTrans) { - diffusivity_.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), - computeHalfMean(halfDiff, porosity_)); - } + // update the "diffusive half transmissibility" for the intersection + if (updateDiffusivity && !onlyTrans) { + diffusivity.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), + computeHalfMean(halfDiff, porosity_)); + } - // update the "dispersivity half transmissibility" for the intersection - if (updateDispersivity && !onlyTrans) { - dispersivity_.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), - computeHalfMean(halfDiff, dispersion_)); + // update the "dispersivity half transmissibility" for the intersection + if (updateDispersivity && !onlyTrans) { + dispersivity.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), + computeHalfMean(halfDiff, dispersion_)); + } } } } - centroids_cache_.clear(); +#ifdef _OPENMP +#pragma omp parallel sections +#endif + { +#ifdef _OPENMP +#pragma omp section +#endif + transMap.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + transBoundary.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + thermalHalfTransBoundary.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + thermalHalfTrans.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + diffusivity.finalize(); +#ifdef _OPENMP +#pragma omp section +#endif + dispersivity.finalize(); + } + // Potentially overwrite and/or modify transmissibilities based on input from deck this->updateFromEclState_(global);