diff --git a/opm/simulators/flow/Transmissibility.hpp b/opm/simulators/flow/Transmissibility.hpp index a7735c522..2256b44a6 100644 --- a/opm/simulators/flow/Transmissibility.hpp +++ b/opm/simulators/flow/Transmissibility.hpp @@ -163,14 +163,20 @@ protected: void removeNonCartesianTransmissibilities_(bool removeAll); + struct FaceInfo + { + DimVector faceCenter; + int faceIdx; + unsigned elemIdx; + unsigned cartElemIdx; + }; + /// \brief Apply the Multipliers for the case PINCH(4)==TOPBOT /// /// \param pinchTop Whether PINCH(5) is TOP, otherwise ALL is assumed. void applyAllZMultipliers_(Scalar& trans, - unsigned insideFaceIdx, - unsigned outsideFaceIdx, - unsigned insideCartElemIdx, - unsigned outsideCartElemIdx, + const FaceInfo& inside, + const FaceInfo& outside, const TransMult& transMult, const std::array& cartDims); @@ -190,23 +196,15 @@ protected: template void computeFaceProperties(const Intersection& intersection, - const int, - const int, - const int, - const int, - DimVector& faceCenterInside, - DimVector& faceCenterOutside, + FaceInfo& inside, + FaceInfo& outside, DimVector& faceAreaNormal, /*isCpGrid=*/std::false_type) const; template void computeFaceProperties(const Intersection& intersection, - const int insideElemIdx, - const int insideFaceIdx, - const int outsideElemIdx, - const int outsideFaceIdx, - DimVector& faceCenterInside, - DimVector& faceCenterOutside, + FaceInfo& inside, + FaceInfo& outside, DimVector& faceAreaNormal, /*isCpGrid=*/std::true_type) const; @@ -251,16 +249,14 @@ protected: void extractDispersion_(); - void computeHalfTrans_(Scalar& halfTrans, - const DimVector& areaNormal, - int faceIdx, // in the reference element that contains the intersection - const DimVector& distance, - const DimMatrix& perm) const; + static Scalar computeHalfTrans_(const DimVector& areaNormal, + int faceIdx, // in the reference element that contains the intersection + const DimVector& distance, + const DimMatrix& perm); - void computeHalfDiffusivity_(Scalar& halfDiff, - const DimVector& areaNormal, - const DimVector& distance, - const Scalar& poro) const; + static Scalar computeHalfDiffusivity_(const DimVector& areaNormal, + const DimVector& distance, + const Scalar poro); DimVector distanceVector_(const DimVector& faceCenter, const unsigned& cellIdx) const; @@ -270,10 +266,9 @@ protected: unsigned cartElemIdx, const TransMult& transMult) const; - void applyNtg_(Scalar& trans, - unsigned faceIdx, - unsigned elemIdx, - const std::vector& ntg) const; + static void applyNtg_(Scalar& trans, + const FaceInfo& face, + const std::vector& ntg); std::vector permeability_; std::vector porosity_; diff --git a/opm/simulators/flow/Transmissibility_impl.hpp b/opm/simulators/flow/Transmissibility_impl.hpp index f2e07fb89..8d38778ec 100644 --- a/opm/simulators/flow/Transmissibility_impl.hpp +++ b/opm/simulators/flow/Transmissibility_impl.hpp @@ -177,10 +177,12 @@ update(bool global, const TransUpdateQuantities update_quantities, const bool disableNNC = eclState_.getSimulationConfig().useNONNC(); - if (map) + if (map) { extractPermeability_(map); - else + } + else { extractPermeability_(); + } // 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 @@ -225,8 +227,7 @@ update(bool global, const TransUpdateQuantities update_quantities, auto pinchTransCalcMode = eclGrid.getPinchOption(); useSmallestMultiplier = eclGrid.getMultzOption() == PinchMode::ALL; pinchOption4ALL = (pinchTransCalcMode == PinchMode::ALL); - if (pinchOption4ALL) - { + if (pinchOption4ALL) { useSmallestMultiplier = false; } } @@ -243,296 +244,236 @@ update(bool global, const TransUpdateQuantities update_quantities, centroids_cache_[elemIdx] = centroids_(elemIdx); } + auto harmonicMean = [](const Scalar half1, const Scalar half2) + { + return (std::abs(half1) < 1e-30 || std::abs(half2) < 1e-30) + ? 0.0 + : 1.0 / (1.0 / half1 + 1.0 / half2); + }; + + auto faceIdToDir = [](int insideFaceIdx) + { + switch (insideFaceIdx) { + case 0: + case 1: + return FaceDir::XPlus; + case 2: + case 3: + return FaceDir::YPlus; + break; + case 4: + case 5: + return FaceDir::ZPlus; + default: + throw std::logic_error("Could not determine a face direction"); + } + }; + + auto halfDiff = [](const DimVector& faceAreaNormal, + const unsigned, + const DimVector& distVector, + const Scalar prop) + { + return computeHalfDiffusivity_(faceAreaNormal, + distVector, + prop); + }; + // compute the transmissibilities for all intersections for (const auto& elem : elements(gridView_)) { - unsigned elemIdx = elemMapper.index(elem); + 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. + 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]); + }; - auto isIt = gridView_.ibegin(elem); - const auto& isEndIt = gridView_.iend(elem); unsigned boundaryIsIdx = 0; - for (; isIt != isEndIt; ++ isIt) { - // store intersection, this might be costly - const auto& intersection = *isIt; - + 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(); - const auto& faceCenterInside = geometry.center(); + inside.faceCenter = geometry.center(); - auto faceAreaNormal = intersection.centerUnitOuterNormal(); + faceAreaNormal = intersection.centerUnitOuterNormal(); faceAreaNormal *= geometry.volume(); - Scalar transBoundaryIs; - computeHalfTrans_(transBoundaryIs, - faceAreaNormal, - intersection.indexInInside(), - distanceVector_(faceCenterInside, - elemIdx), - permeability_[elemIdx]); + 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. - unsigned insideCartElemIdx = cartMapper_.cartesianIndex(elemIdx); - applyMultipliers_(transBoundaryIs, intersection.indexInInside(), insideCartElemIdx, transMult); - transBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = transBoundaryIs; + 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_(transBoundaryEnergyIs, - faceAreaNormal, - distanceVector_(faceCenterInside, - elemIdx), - 1.0); - thermalHalfTransBoundary_[std::make_pair(elemIdx, boundaryIsIdx)] = - transBoundaryEnergyIs; + Scalar transBoundaryEnergyIs = + computeHalfDiffusivity_(faceAreaNormal, + distanceVector_(inside.faceCenter, inside.elemIdx), + 1.0); + thermalHalfTransBoundary_.insert_or_assign(std::make_pair(inside.elemIdx, boundaryIsIdx), + transBoundaryEnergyIs); } - ++ boundaryIsIdx; + ++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; + ++boundaryIsIdx; continue; } const auto& outsideElem = intersection.outside(); - unsigned outsideElemIdx = elemMapper.index(outsideElem); + outside.elemIdx = elemMapper.index(outsideElem); - // Get the Cartesian indices of the origen cells (parent or equivalent cell on level zero), for CpGrid with LGRs. - // For genral grids and no LGRs, get the usual Cartesian Index. - unsigned insideCartElemIdx = this-> lookUpCartesianData_.template getFieldPropCartesianIdx(elemIdx); - unsigned outsideCartElemIdx = this-> lookUpCartesianData_.template getFieldPropCartesianIdx(outsideElemIdx); + // 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 insideCartElemIdx>outsideCartElemIdx does not imply elemIdx>outsideElemIdx for + // 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(insideCartElemIdx, elemIdx) > std::tie(outsideCartElemIdx, outsideElemIdx)) + if (std::tie(inside.cartElemIdx, inside.elemIdx) > std::tie(outside.cartElemIdx, outside.elemIdx)) { continue; + } // local indices of the faces of the inside and // outside elements which contain the intersection - int insideFaceIdx = intersection.indexInInside(); - int outsideFaceIdx = intersection.indexInOutside(); + inside.faceIdx = intersection.indexInInside(); + outside.faceIdx = intersection.indexInOutside(); - if (insideFaceIdx == -1) { + if (inside.faceIdx == -1) { // NNC. Set zero transmissibility, as it will be // *added to* by applyNncToGridTrans_() later. - assert(outsideFaceIdx == -1); - trans_[details::isId(elemIdx, outsideElemIdx)] = 0.0; - if (enableEnergy_ && !onlyTrans){ - thermalHalfTrans_[details::directionalIsId(elemIdx, outsideElemIdx)] = 0.0; - thermalHalfTrans_[details::directionalIsId(outsideElemIdx, elemIdx)] = 0.0; + 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); } if (updateDiffusivity && !onlyTrans) { - diffusivity_[details::isId(elemIdx, outsideElemIdx)] = 0.0; + diffusivity_.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), 0.0); } if (updateDispersivity && !onlyTrans) { - dispersivity_[details::isId(elemIdx, outsideElemIdx)] = 0.0; + dispersivity_.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), 0.0); } continue; } - DimVector faceCenterInside; - DimVector faceCenterOutside; - DimVector faceAreaNormal; - typename std::is_same::type isCpGrid; computeFaceProperties(intersection, - elemIdx, - insideFaceIdx, - outsideElemIdx, - outsideFaceIdx, - faceCenterInside, - faceCenterOutside, + inside, + outside, faceAreaNormal, isCpGrid); - Scalar halfTrans1; - Scalar halfTrans2; - - computeHalfTrans_(halfTrans1, - faceAreaNormal, - insideFaceIdx, - distanceVector_(faceCenterInside, - elemIdx), - permeability_[elemIdx]); - computeHalfTrans_(halfTrans2, - faceAreaNormal, - outsideFaceIdx, - distanceVector_(faceCenterOutside, - outsideElemIdx), - permeability_[outsideElemIdx]); - - applyNtg_(halfTrans1, insideFaceIdx, elemIdx, ntg); - applyNtg_(halfTrans2, outsideFaceIdx, outsideElemIdx, ntg); - - // convert half transmissibilities to full face - // transmissibilities using the harmonic mean - Scalar trans; - if (std::abs(halfTrans1) < 1e-30 || std::abs(halfTrans2) < 1e-30) - // avoid division by zero - trans = 0.0; - else - trans = 1.0 / (1.0/halfTrans1 + 1.0/halfTrans2); + Scalar trans = computeHalfMean(computeHalfTrans_, permeability_); // apply the full face transmissibility multipliers // for the inside ... - if(!pinchActive){ - if (insideFaceIdx > 3){// top or bottom - auto find_layer = [&cartDims](std::size_t cell){ + 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(insideCartElemIdx); - int kdown=find_layer(outsideCartElemIdx); + 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) || (insideCartElemIdx == outsideCartElemIdx)); - if(std::abs(kup -kdown) > 1){ + assert((kup != kdown) || (inside.cartElemIdx == outside.cartElemIdx)); + if (std::abs(kup -kdown) > 1) { trans = 0.0; } } } - if (useSmallestMultiplier) - { + 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, insideFaceIdx, outsideFaceIdx, insideCartElemIdx, - outsideCartElemIdx, transMult, cartDims); + applyAllZMultipliers_(trans, inside, outside, transMult, cartDims); } - else - { - applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult); + else { + applyMultipliers_(trans, inside.faceIdx, inside.cartElemIdx, transMult); // ... and outside elements - applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult); + applyMultipliers_(trans, outside.faceIdx, outside.cartElemIdx, transMult); } // apply the region multipliers (cf. the MULTREGT keyword) - FaceDir::DirEnum faceDir; - switch (insideFaceIdx) { - case 0: - case 1: - faceDir = FaceDir::XPlus; - break; + trans *= transMult.getRegionMultiplier(inside.cartElemIdx, + outside.cartElemIdx, + faceIdToDir(inside.faceIdx)); - case 2: - case 3: - faceDir = FaceDir::YPlus; - break; - - case 4: - case 5: - faceDir = FaceDir::ZPlus; - break; - - default: - throw std::logic_error("Could not determine a face direction"); - } - - trans *= transMult.getRegionMultiplier(insideCartElemIdx, - outsideCartElemIdx, - faceDir); - - trans_[details::isId(elemIdx, outsideElemIdx)] = trans; + trans_.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), trans); // update the "thermal half transmissibility" for the intersection if (enableEnergy_ && !onlyTrans) { - - Scalar halfDiffusivity1; - Scalar halfDiffusivity2; - - computeHalfDiffusivity_(halfDiffusivity1, - faceAreaNormal, - distanceVector_(faceCenterInside, - elemIdx), - 1.0); - computeHalfDiffusivity_(halfDiffusivity2, - faceAreaNormal, - distanceVector_(faceCenterOutside, - outsideElemIdx), - 1.0); - //TODO Add support for multipliers - thermalHalfTrans_[details::directionalIsId(elemIdx, outsideElemIdx)] = halfDiffusivity1; - thermalHalfTrans_[details::directionalIsId(outsideElemIdx, elemIdx)] = halfDiffusivity2; - } + 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_)); + } - Scalar halfDiffusivity1; - Scalar halfDiffusivity2; - - computeHalfDiffusivity_(halfDiffusivity1, - faceAreaNormal, - distanceVector_(faceCenterInside, - elemIdx), - porosity_[elemIdx]); - computeHalfDiffusivity_(halfDiffusivity2, - faceAreaNormal, - distanceVector_(faceCenterOutside, - outsideElemIdx), - porosity_[outsideElemIdx]); - - applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg); - applyNtg_(halfDiffusivity2, outsideFaceIdx, outsideElemIdx, ntg); - - //TODO Add support for multipliers - Scalar diffusivity; - if (std::abs(halfDiffusivity1) < 1e-30 || std::abs(halfDiffusivity2) < 1e-30) - // avoid division by zero - diffusivity = 0.0; - else - diffusivity = 1.0 / (1.0/halfDiffusivity1 + 1.0/halfDiffusivity2); - - - diffusivity_[details::isId(elemIdx, outsideElemIdx)] = diffusivity; - } - - // update the "dispersivity half transmissibility" for the intersection + // update the "dispersivity half transmissibility" for the intersection if (updateDispersivity && !onlyTrans) { - - Scalar halfDispersivity1; - Scalar halfDispersivity2; - - computeHalfDiffusivity_(halfDispersivity1, - faceAreaNormal, - distanceVector_(faceCenterInside, - elemIdx), - dispersion_[elemIdx]); - computeHalfDiffusivity_(halfDispersivity2, - faceAreaNormal, - distanceVector_(faceCenterOutside, - outsideElemIdx), - dispersion_[outsideElemIdx]); - - applyNtg_(halfDispersivity1, insideFaceIdx, elemIdx, ntg); - applyNtg_(halfDispersivity2, outsideFaceIdx, outsideElemIdx, ntg); - - //TODO Add support for multipliers - Scalar dispersivity; - if (std::abs(halfDispersivity1) < 1e-30 || std::abs(halfDispersivity2) < 1e-30) - // avoid division by zero - dispersivity = 0.0; - else - dispersivity = 1.0 / (1.0/halfDispersivity1 + 1.0/halfDispersivity2); - - - dispersivity_[details::isId(elemIdx, outsideElemIdx)] = dispersivity; - } + dispersivity_.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), + computeHalfMean(halfDiff, dispersion_)); + } } } @@ -627,19 +568,24 @@ extractPermeability_(const std::function& map) // over several processes.) const auto& fp = eclState_.fieldProps(); if (fp.has_double("PERMX")) { - const std::vector& permxData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PERMX"); + const std::vector& permxData = + this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PERMX"); std::vector permyData; - if (fp.has_double("PERMY")) - permyData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PERMY"); - else + if (fp.has_double("PERMY")){ + permyData = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PERMY"); + } + else { permyData = permxData; + } std::vector permzData; - if (fp.has_double("PERMZ")) - permzData = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PERMZ"); - else + if (fp.has_double("PERMZ")) { + permzData = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PERMZ"); + } + else { permzData = permxData; + } for (std::size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) { permeability_[dofIdx] = 0.0; @@ -650,11 +596,11 @@ extractPermeability_(const std::function& map) } // for now we don't care about non-diagonal entries - } - else + else { throw std::logic_error("Can't read the intrinsic permeability from the ecl state. " "(The PERM{X,Y,Z} keywords are missing)"); + } } template @@ -668,16 +614,18 @@ extractPorosity_() const auto& fp = eclState_.fieldProps(); if (fp.has_double("PORO")) { if constexpr (std::is_same_v) { - porosity_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PORO"); - } else { - const auto por = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PORO"); + porosity_ = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PORO"); + } + else { + const auto por = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PORO"); porosity_.resize(por.size()); std::copy(por.begin(), por.end(), porosity_.begin()); } } - else + else { throw std::logic_error("Can't read the porosityfrom the ecl state. " "(The PORO keywords are missing)"); + } } template @@ -690,9 +638,10 @@ extractDispersion_() } const auto& fp = eclState_.fieldProps(); if constexpr (std::is_same_v) { - dispersion_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"DISPERC"); - } else { - const auto disp = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"DISPERC"); + dispersion_ = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"DISPERC"); + } + else { + const auto disp = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"DISPERC"); dispersion_.resize(disp.size()); std::copy(disp.begin(), disp.end(), dispersion_.begin()); } @@ -705,7 +654,7 @@ removeNonCartesianTransmissibilities_(bool removeAll) const auto& cartDims = cartMapper_.cartesianDimensions(); for (auto&& trans: trans_) { //either remove all NNC transmissibilities or those less than the threshold (by default 1e-6 in the deck's unit system) - if (removeAll or trans.second < transmissibilityThreshold_) { + if (removeAll || trans.second < transmissibilityThreshold_) { const auto& id = trans.first; const auto& elements = details::isIdReverse(id); int gc1 = std::min(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second)); @@ -714,8 +663,9 @@ removeNonCartesianTransmissibilities_(bool removeAll) // only adjust the NNCs // When LGRs, all neighbors in the LGR are cartesian neighbours on the level grid representing the LGR. // When elements on the leaf grid view have the same parent cell, gc1 and gc2 coincide. - if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1] || gc2 - gc1 == 0) + if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1] || gc2 - gc1 == 0) { continue; + } trans.second = 0.0; } @@ -725,35 +675,32 @@ removeNonCartesianTransmissibilities_(bool removeAll) template void Transmissibility:: applyAllZMultipliers_(Scalar& trans, - unsigned insideFaceIdx, - unsigned outsideFaceIdx, - unsigned insideCartElemIdx, - unsigned outsideCartElemIdx, + const FaceInfo& inside, + const FaceInfo& outside, const TransMult& transMult, const std::array& cartDims) { - if(grid_.maxLevel()> 0) { - OPM_THROW(std::invalid_argument, "MULTZ not support with LGRS, yet."); + if (grid_.maxLevel() > 0) { + OPM_THROW(std::invalid_argument, "MULTZ not support with LGRS, yet."); } - if (insideFaceIdx > 3) { // top or or bottom - assert(insideFaceIdx==5); // as insideCartElemIdx < outsideCartElemIdx holds for the Z column + if (inside.faceIdx > 3) { // top or or bottom + assert(inside.faceIdx == 5); // as insideCartElemIdx < outsideCartElemIdx holds for the Z column // For CpGrid with LGRs, insideCartElemIdx == outsideCartElemIdx when cells on the leaf have the same parent cell on level zero. - assert(outsideCartElemIdx >= insideCartElemIdx); + assert(outside.cartElemIdx >= inside.cartElemIdx); unsigned lastCartElemIdx; - if (outsideCartElemIdx == insideCartElemIdx) { - lastCartElemIdx = outsideCartElemIdx; + if (outside.cartElemIdx == inside.cartElemIdx) { + lastCartElemIdx = outside.cartElemIdx; } else { - lastCartElemIdx = outsideCartElemIdx - cartDims[0]*cartDims[1]; + lastCartElemIdx = outside.cartElemIdx - cartDims[0]*cartDims[1]; } // Last multiplier using (Z+)*(Z-) Scalar mult = transMult.getMultiplier(lastCartElemIdx , FaceDir::ZPlus) * - transMult.getMultiplier(outsideCartElemIdx , FaceDir::ZMinus); + transMult.getMultiplier(outside.cartElemIdx , FaceDir::ZMinus); // pick the smallest multiplier using (Z+)*(Z-) while looking down // the pillar until reaching the other end of the connection - for(auto cartElemIdx = insideCartElemIdx; cartElemIdx < lastCartElemIdx;) - { + for (auto cartElemIdx = inside.cartElemIdx; cartElemIdx < lastCartElemIdx;) { auto multiplier = transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus); cartElemIdx += cartDims[0]*cartDims[1]; multiplier *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus); @@ -762,10 +709,9 @@ applyAllZMultipliers_(Scalar& trans, trans *= mult; } - else - { - applyMultipliers_(trans, insideFaceIdx, insideCartElemIdx, transMult); - applyMultipliers_(trans, outsideFaceIdx, outsideCartElemIdx, transMult); + else { + applyMultipliers_(trans, inside.faceIdx, inside.cartElemIdx, transMult); + applyMultipliers_(trans, outside.faceIdx, outside.cartElemIdx, transMult); } } @@ -781,8 +727,7 @@ updateFromEclState_(bool global) fp->tran_active("TRANY"), fp->tran_active("TRANZ")}; - if( !(is_tran[0] ||is_tran[1] || is_tran[2]) ) - { + if (!(is_tran[0] || is_tran[1] || is_tran[2])) { // Skip unneeded expensive traversals return; } @@ -792,13 +737,12 @@ updateFromEclState_(bool global) auto key = keywords.begin(); auto perform = is_tran.begin(); - for (auto it = trans.begin(); it != trans.end(); ++it, ++key, ++perform) - { - if(*perform) { - if(grid_.maxLevel()>0) { + for (auto it = trans.begin(); it != trans.end(); ++it, ++key, ++perform) { + if (*perform) { + if (grid_.maxLevel() > 0) { OPM_THROW(std::invalid_argument, "Calculations on TRANX/TRANY/TRANZ arrays are not support with LGRS, yet."); } - fp->apply_tran(*key, *it); + fp->apply_tran(*key, *it); } } @@ -814,17 +758,19 @@ createTransmissibilityArrays_(const std::array& is_tran) ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout()); auto numElem = gridView_.size(/*codim=*/0); - std::array,3> trans = - { std::vector(is_tran[0] ? numElem : 0, 0), + std::array,3> trans = { + std::vector(is_tran[0] ? numElem : 0, 0), std::vector(is_tran[1] ? numElem : 0, 0), - std::vector(is_tran[2] ? numElem : 0, 0)}; + std::vector(is_tran[2] ? numElem : 0, 0) + }; // compute the transmissibilities for all intersections for (const auto& elem : elements(gridView_)) { for (const auto& intersection : intersections(gridView_, elem)) { // store intersection, this might be costly - if (!intersection.neighbor()) + if (!intersection.neighbor()) { continue; // intersection is on the domain boundary + } // In the EclState TRANX[c1] is transmissibility in X+ // direction. we only store transmissibilities in the + @@ -841,12 +787,13 @@ createTransmissibilityArrays_(const std::array& is_tran) unsigned c2 = elemMapper.index(intersection.outside()); int gc1 = cartMapper_.cartesianIndex(c1); int gc2 = cartMapper_.cartesianIndex(c2); - if (std::tie(gc1, c1) > std::tie(gc2, c2)) + if (std::tie(gc1, c1) > std::tie(gc2, c2)) { // we only need to handle each connection once, thank you. // We do this when gc1 is smaller than the other to find the // correct place to store in parallel when ghost/overlap elements // are ordered last continue; + } auto isID = details::isId(c1, c2); @@ -855,25 +802,30 @@ createTransmissibilityArrays_(const std::array& is_tran) // 'intersection.indexInSIde()' needed to be checked to determine the direction, i.e. // add in the if/else-if 'gc2 == gc1 && intersection.indexInInside() == ... ' if ((gc2 - gc1 == 1 || (gc2 == gc1 && (intersection.indexInInside() == 0 || intersection.indexInInside() == 1))) - && cartDims[0] > 1) { - if (is_tran[0]) + && cartDims[0] > 1) + { + if (is_tran[0]) { // set simulator internal transmissibilities to values from inputTranx trans[0][c1] = trans_[isID]; + } } else if ((gc2 - gc1 == cartDims[0] || (gc2 == gc1 && (intersection.indexInInside() == 2 || intersection.indexInInside() == 3))) - && cartDims[1] > 1) { - if (is_tran[1]) + && cartDims[1] > 1) + { + if (is_tran[1]) { // set simulator internal transmissibilities to values from inputTrany trans[1][c1] = trans_[isID]; + } } else if (gc2 - gc1 == cartDims[0]*cartDims[1] || - (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) { - if (is_tran[2]) + (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) + { + if (is_tran[2]) { // set simulator internal transmissibilities to values from inputTranz trans[2][c1] = trans_[isID]; + } } - - //else.. We don't support modification of NNC at the moment. + // else.. We don't support modification of NNC at the moment. } } @@ -891,8 +843,9 @@ resetTransmissibilityFromArrays_(const std::array& is_tran, // compute the transmissibilities for all intersections for (const auto& elem : elements(gridView_)) { for (const auto& intersection : intersections(gridView_, elem)) { - if (!intersection.neighbor()) + if (!intersection.neighbor()) { continue; // intersection is on the domain boundary + } // In the EclState TRANX[c1] is transmissibility in X+ // direction. we only store transmissibilities in the + @@ -909,12 +862,13 @@ resetTransmissibilityFromArrays_(const std::array& is_tran, unsigned c2 = elemMapper.index(intersection.outside()); int gc1 = cartMapper_.cartesianIndex(c1); int gc2 = cartMapper_.cartesianIndex(c2); - if (std::tie(gc1, c1) > std::tie(gc2, c2)) + if (std::tie(gc1, c1) > std::tie(gc2, c2)) { // we only need to handle each connection once, thank you. // We do this when gc1 is smaller than the other to find the // correct place to read in parallel when ghost/overlap elements // are ordered last continue; + } auto isID = details::isId(c1, c2); @@ -923,25 +877,31 @@ resetTransmissibilityFromArrays_(const std::array& is_tran, // 'intersection.indexInSIde()' needed to be checked to determine the direction, i.e. // add in the if/else-if 'gc2 == gc1 && intersection.indexInInside() == ... ' if ((gc2 - gc1 == 1 || (gc2 == gc1 && (intersection.indexInInside() == 0 || intersection.indexInInside() == 1))) - && cartDims[0] > 1) { - if (is_tran[0]) + && cartDims[0] > 1) + { + if (is_tran[0]) { // set simulator internal transmissibilities to values from inputTranx trans_[isID] = trans[0][c1]; + } } else if ((gc2 - gc1 == cartDims[0] || (gc2 == gc1 && (intersection.indexInInside() == 2|| intersection.indexInInside() == 3))) - && cartDims[1] > 1) { - if (is_tran[1]) + && cartDims[1] > 1) + { + if (is_tran[1]) { // set simulator internal transmissibilities to values from inputTrany trans_[isID] = trans[1][c1]; + } } else if (gc2 - gc1 == cartDims[0]*cartDims[1] || - (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) { - if (is_tran[2]) + (gc2 == gc1 && (intersection.indexInInside() == 4 || intersection.indexInInside() == 5))) + { + if (is_tran[2]) { // set simulator internal transmissibilities to values from inputTranz trans_[isID] = trans[2][c1]; + } } - //else.. We don't support modification of NNC at the moment. + // else.. We don't support modification of NNC at the moment. } } } @@ -950,19 +910,14 @@ template void Transmissibility:: computeFaceProperties(const Intersection& intersection, - const int, - const int, - const int, - const int, - DimVector& faceCenterInside, - DimVector& faceCenterOutside, + FaceInfo& inside, + FaceInfo& outside, DimVector& faceAreaNormal, /*isCpGrid=*/std::false_type) const { // default implementation for DUNE grids const auto& geometry = intersection.geometry(); - faceCenterInside = geometry.center(); - faceCenterOutside = faceCenterInside; + outside.faceCenter = inside.faceCenter = geometry.center(); faceAreaNormal = intersection.centerUnitOuterNormal(); faceAreaNormal *= geometry.volume(); @@ -972,43 +927,41 @@ template void Transmissibility:: computeFaceProperties(const Intersection& intersection, - const int insideElemIdx, - const int insideFaceIdx, - const int outsideElemIdx, - const int outsideFaceIdx, - DimVector& faceCenterInside, - DimVector& faceCenterOutside, + FaceInfo& inside, + FaceInfo& outside, DimVector& faceAreaNormal, /*isCpGrid=*/std::true_type) const { int faceIdx = intersection.id(); - if(grid_.maxLevel() == 0) { - faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection); - faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection); + if (grid_.maxLevel() == 0) { + inside.faceCenter = grid_.faceCenterEcl(inside.elemIdx, inside.faceIdx, intersection); + outside.faceCenter = grid_.faceCenterEcl(outside.elemIdx, outside.faceIdx, intersection); faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx); } else { if ((intersection.inside().level() != intersection.outside().level())) { - // For CpGrid with LGRs, intersection laying on the boundary of an LGR, having two neighboring cells: // one coarse neighboring cell and one refined neighboring cell, we get the corresponding parent // intersection (from level 0), and use the center of the parent intersection for the coarse // neighboring cell. // Get parent intersection and its geometry - const auto& parentIntersection = grid_.getParentIntersectionFromLgrBoundaryFace(intersection); + const auto& parentIntersection = + grid_.getParentIntersectionFromLgrBoundaryFace(intersection); const auto& parentIntersectionGeometry = parentIntersection.geometry(); // For the coarse neighboring cell, take the center of the parent intersection. - // For the refined neighboring cell, take the 'usual' center. - faceCenterInside = (intersection.inside().level() == 0) ? parentIntersectionGeometry.center() : - grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection); - faceCenterOutside = (intersection.outside().level() == 0) ? parentIntersectionGeometry.center() : - grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection); + // For the refined neighboring cell, take the 'usual' center. + inside.faceCenter = (intersection.inside().level() == 0) + ? parentIntersectionGeometry.center() + : grid_.faceCenterEcl(inside.elemIdx, inside.faceIdx, intersection); + outside.faceCenter = (intersection.outside().level() == 0) + ? parentIntersectionGeometry.center() + : grid_.faceCenterEcl(outside.elemIdx, outside.faceIdx, intersection); // For some computations, it seems to be benefitial to replace the actual area of the refined face, by - // the area of its parent face. + // the area of its parent face. // faceAreaNormal = parentIntersection.centerUnitOuterNormal(); // faceAreaNormal *= parentIntersectionGeometry.volume(); @@ -1018,10 +971,10 @@ computeFaceProperties(const Intersection& intersection, } else { assert(intersection.inside().level() == intersection.outside().level()); - - faceCenterInside = grid_.faceCenterEcl(insideElemIdx, insideFaceIdx, intersection); - faceCenterOutside = grid_.faceCenterEcl(outsideElemIdx, outsideFaceIdx, intersection); - + + inside.faceCenter = grid_.faceCenterEcl(inside.elemIdx, inside.faceIdx, intersection); + outside.faceCenter = grid_.faceCenterEcl(outside.elemIdx, outside.faceIdx, intersection); + // When the CpGrid has LGRs, we compute the face area normal differently. if (intersection.inside().level() > 0) { // remove intersection.inside().level() > 0 faceAreaNormal = intersection.centerUnitOuterNormal(); @@ -1033,6 +986,7 @@ computeFaceProperties(const Intersection& intersection, } } } + template void Transmissibility:: @@ -1049,12 +1003,14 @@ applyPinchNncToGridTrans_(const std::unordered_map& cartesianTo int low = (lowIt == cartesianToCompressed.end())? -1 : lowIt->second; int high = (highIt == cartesianToCompressed.end())? -1 : highIt->second; - if (low > high) + if (low > high) { std::swap(low, high); + } - if (low == -1 && high == -1) + if (low == -1 && high == -1) { // Silently discard as it is not between active cells continue; + } if (low == -1 || high == -1) { // We can end up here if one of the cells is overlap/ghost, because those @@ -1090,12 +1046,14 @@ applyNncToGridTrans_(const std::unordered_map& cartesianToCompr int low = (lowIt == cartesianToCompressed.end())? -1 : lowIt->second; int high = (highIt == cartesianToCompressed.end())? -1 : highIt->second; - if (low > high) + if (low > high) { std::swap(low, high); + } - if (low == -1 && high == -1) + if (low == -1 && high == -1) { // Silently discard as it is not between active cells continue; + } if (low == -1 || high == -1) { // Discard the NNC if it is between active cell and inactive cell @@ -1106,14 +1064,11 @@ applyNncToGridTrans_(const std::unordered_map& cartesianToCompr continue; } - { - auto candidate = trans_.find(details::isId(low, high)); - if (candidate != trans_.end()) { - // NNC is represented by the grid and might be a neighboring connection - // In this case the transmissibilty is added to the value already - // set or computed. - candidate->second += nncEntry.trans; - } + if (auto candidate = trans_.find(details::isId(low, high)); candidate != trans_.end()) { + // NNC is represented by the grid and might be a neighboring connection + // In this case the transmissibilty is added to the value already + // set or computed. + candidate->second += nncEntry.trans; } // if (enableEnergy_) { // auto candidate = thermalHalfTrans_.find(details::directionalIsId(low, high)); @@ -1177,11 +1132,13 @@ applyEditNncToGridTransHelper_(const std::unordered_map& global const std::function& getLocation, const std::function& apply) { - if (nncs.empty()) + if (nncs.empty()) { return; + } const auto& cartDims = cartMapper_.cartesianDimensions(); - auto format_ijk = [&cartDims](std::size_t cell) -> std::string { + auto format_ijk = [&cartDims](std::size_t cell) -> std::string + { auto i = cell % cartDims[0]; cell /= cartDims[0]; auto j = cell % cartDims[1]; auto k = cell / cartDims[1]; @@ -1189,13 +1146,14 @@ applyEditNncToGridTransHelper_(const std::unordered_map& global return fmt::format("({},{},{})", i + 1,j + 1,k + 1); }; - auto print_warning = [&format_ijk, &getLocation, &keyword] (const NNCdata& nnc) { - const auto& location = getLocation( nnc ); - auto warning = fmt::format("Problem with {} keyword\n" - "In {} line {} \n" - "No NNC defined for connection {} -> {}", keyword, location.filename, - location.lineno, format_ijk(nnc.cell1), format_ijk(nnc.cell2)); - OpmLog::warning(keyword, warning); + auto print_warning = [&format_ijk, &getLocation, &keyword] (const NNCdata& nnc) + { + const auto& location = getLocation( nnc ); + auto warning = fmt::format("Problem with {} keyword\n" + "In {} line {} \n" + "No NNC defined for connection {} -> {}", keyword, location.filename, + location.lineno, format_ijk(nnc.cell1), format_ijk(nnc.cell2)); + OpmLog::warning(keyword, warning); }; // editNnc is supposed to only reference non-neighboring connections and not @@ -1212,7 +1170,7 @@ applyEditNncToGridTransHelper_(const std::unordered_map& global if (lowIt == globalToLocal.end() || highIt == globalToLocal.end()) { // Prevent warnings for NNCs stored on other processes in parallel (both cells inactive) - if ( lowIt != highIt && warnEditNNC_) { + if (lowIt != highIt && warnEditNNC_) { print_warning(*nnc); warning_count++; } @@ -1222,8 +1180,9 @@ applyEditNncToGridTransHelper_(const std::unordered_map& global auto low = lowIt->second, high = highIt->second; - if (low > high) + if (low > high) { std::swap(low, high); + } auto candidate = trans_.find(details::isId(low, high)); if (candidate == trans_.end() && warnEditNNC_) { @@ -1233,7 +1192,7 @@ applyEditNncToGridTransHelper_(const std::unordered_map& global } else { // NNC exists - while (nnc!= end && c1==nnc->cell1 && c2==nnc->cell2) { + while (nnc != end && c1 == nnc->cell1 && c2 == nnc->cell2) { apply(candidate->second, nnc->trans); ++nnc; } @@ -1272,7 +1231,7 @@ applyNncMultreg_(const std::unordered_map& cartesianToCompresse // those act as regular multipliers and have already been fully // accounted for in the multiplier part of the main loop of update() and // the applyEditNncToGridTrans_() member function. - for (const auto& nncList : { &NNC::input, &NNC::editr }) { + for (const auto& nncList : {&NNC::input, &NNC::editr}) { for (const auto& nncEntry : (inputNNC.*nncList)()) { const auto c1 = nncEntry.cell1; const auto c2 = nncEntry.cell2; @@ -1297,31 +1256,35 @@ applyNncMultreg_(const std::unordered_map& cartesianToCompresse } template -void Transmissibility:: -computeHalfTrans_(Scalar& halfTrans, - const DimVector& areaNormal, +Scalar +Transmissibility:: +computeHalfTrans_(const DimVector& areaNormal, int faceIdx, // in the reference element that contains the intersection const DimVector& distance, - const DimMatrix& perm) const + const DimMatrix& perm) { assert(faceIdx >= 0); - unsigned dimIdx = faceIdx/2; + unsigned dimIdx = faceIdx / 2; assert(dimIdx < dimWorld); - halfTrans = perm[dimIdx][dimIdx]; + Scalar halfTrans = perm[dimIdx][dimIdx]; halfTrans *= std::abs(Dune::dot(areaNormal, distance)); halfTrans /= distance.two_norm2(); + + return halfTrans; } template -void Transmissibility:: -computeHalfDiffusivity_(Scalar& halfDiff, - const DimVector& areaNormal, +Scalar +Transmissibility:: +computeHalfDiffusivity_(const DimVector& areaNormal, const DimVector& distance, - const Scalar& poro) const + const Scalar poro) { - halfDiff = poro; + Scalar halfDiff = poro; halfDiff *= std::abs(Dune::dot(areaNormal, distance)); halfDiff /= distance.two_norm2(); + + return halfDiff; } template @@ -1333,8 +1296,9 @@ distanceVector_(const DimVector& faceCenter, const auto& cellCenter = centroids_cache_.empty() ? centroids_(cellIdx) : centroids_cache_[cellIdx]; DimVector x = faceCenter; - for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) { x[dimIdx] -= cellCenter[dimIdx]; + } return x; } @@ -1346,59 +1310,25 @@ applyMultipliers_(Scalar& trans, unsigned cartElemIdx, const TransMult& transMult) const { - // apply multiplyer for the transmissibility of the face. (the + // apply multiplier for the transmissibility of the face. (the // face index is the index of the reference-element face which // contains the intersection of interest.) - switch (faceIdx) { - case 0: // left - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XMinus); - break; - case 1: // right - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::XPlus); - break; - - case 2: // front - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YMinus); - break; - case 3: // back - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::YPlus); - break; - - case 4: // bottom - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus); - break; - case 5: // top - trans *= transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus); - break; - } + trans *= transMult.getMultiplier(cartElemIdx, + FaceDir::FromIntersectionIndex(faceIdx)); } template void Transmissibility:: applyNtg_(Scalar& trans, - unsigned faceIdx, - unsigned elemIdx, - const std::vector& ntg) const + const FaceInfo& face, + const std::vector& ntg) { - // apply multiplyer for the transmissibility of the face. (the + // apply multiplier for the transmissibility of the face. (the // face index is the index of the reference-element face which // contains the intersection of interest.) - switch (faceIdx) { - case 0: // left - trans *= ntg[elemIdx]; - break; - case 1: // right - trans *= ntg[elemIdx]; - break; - - case 2: // front - trans *= ntg[elemIdx]; - break; - case 3: // back - trans *= ntg[elemIdx]; - break; - - // NTG does not apply to top and bottom faces + // NTG does not apply to top and bottom faces + if (face.faceIdx >= 0 && face.faceIdx <= 3) { + trans *= ntg[face.elemIdx]; } }