From f5fa602abb65e860fc50afb3dc6b786001ec3321 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 5 May 2021 12:13:25 +0200 Subject: [PATCH] eclthresholdpressure: split in typetag dependent and typetag-independent parts this allows using explicit template instantation to only compile this code per grid, not per simulator object --- CMakeLists_files.cmake | 1 + ebos/eclgenericthresholdpressure.cc | 265 ++++++++++++++++++++++++++ ebos/eclgenericthresholdpressure.hh | 111 +++++++++++ ebos/eclthresholdpressure.hh | 277 ++++------------------------ 4 files changed, 408 insertions(+), 246 deletions(-) create mode 100644 ebos/eclgenericthresholdpressure.cc create mode 100644 ebos/eclgenericthresholdpressure.hh diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index a34a1e9b5..23f2671f2 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -25,6 +25,7 @@ list (APPEND MAIN_SOURCE_FILES ebos/collecttoiorank.cc ebos/eclgenericcpgridvanguard.cc + ebos/eclgenericthresholdpressure.cc ebos/eclgenericvanguard.cc ebos/ecltransmissibility.cc opm/core/props/phaseUsageFromDeck.cpp diff --git a/ebos/eclgenericthresholdpressure.cc b/ebos/eclgenericthresholdpressure.cc new file mode 100644 index 000000000..66a7bbfbc --- /dev/null +++ b/ebos/eclgenericthresholdpressure.cc @@ -0,0 +1,265 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#if HAVE_DUNE_FEM +#include +#include +#include +#endif + +#include +#include +#include + +namespace Opm { + +template +EclGenericThresholdPressure:: +EclGenericThresholdPressure(const CartesianIndexMapper& cartMapper, + const GridView& gridView, + const ElementMapper& elementMapper, + const EclipseState& eclState, + const Deck& deck, + bool enableExperiments) + : cartMapper_(cartMapper) + , gridView_(gridView) + , elementMapper_(elementMapper) + , eclState_(eclState) + , deck_(deck) + , enableExperiments_(enableExperiments) +{ +} + +template +Scalar EclGenericThresholdPressure:: +thresholdPressure(int elem1Idx, int elem2Idx) const +{ + if (!enableThresholdPressure_) + return 0.0; + + if (enableExperiments_) { + // threshold pressure accross faults + if (!thpresftValues_.empty()) { + int cartElem1Idx = cartMapper_.cartesianIndex(elem1Idx); + int cartElem2Idx = cartMapper_.cartesianIndex(elem2Idx); + + assert(0 <= cartElem1Idx && static_cast(cartElemFaultIdx_.size()) > cartElem1Idx); + assert(0 <= cartElem2Idx && static_cast(cartElemFaultIdx_.size()) > cartElem2Idx); + + int fault1Idx = cartElemFaultIdx_[cartElem1Idx]; + int fault2Idx = cartElemFaultIdx_[cartElem2Idx]; + if (fault1Idx != -1 && fault1Idx == fault2Idx) + // inside a fault there's no threshold pressure, even accross EQUIL + // regions. + return 0.0; + if (fault1Idx != fault2Idx) { + // TODO: which value if a cell is part of multiple faults? we take + // the maximum here. + Scalar val1 = (fault1Idx >= 0) ? thpresftValues_[fault1Idx] : 0.0; + Scalar val2 = (fault2Idx >= 0) ? thpresftValues_[fault2Idx] : 0.0; + return std::max(val1, val2); + } + } + } + + // threshold pressure accross EQUIL regions + unsigned short equilRegion1Idx = elemEquilRegion_[elem1Idx]; + unsigned short equilRegion2Idx = elemEquilRegion_[elem2Idx]; + + if (equilRegion1Idx == equilRegion2Idx) + return 0.0; + + return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx]; +} + +template +void EclGenericThresholdPressure:: +finishInit() +{ + unsigned numElements = gridView_.size(/*codim=*/0); + + const auto& simConfig = eclState_.getSimulationConfig(); + + enableThresholdPressure_ = simConfig.useThresholdPressure(); + if (!enableThresholdPressure_) + return; + + numEquilRegions_ = eclState_.getTableManager().getEqldims().getNumEquilRegions(); + if (numEquilRegions_ > 0xff) { + // make sure that the index of an equilibration region can be stored in a + // single byte + throw std::runtime_error("The maximum number of supported equilibration regions is 255!"); + } + + // internalize the data specified using the EQLNUM keyword + const auto& fp = eclState_.fieldProps(); + const auto& equilRegionData = fp.get_int("EQLNUM"); + elemEquilRegion_.resize(numElements, 0); + for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { + elemEquilRegion_[elemIdx] = equilRegionData[elemIdx] - 1; + } + + /* + If this is a restart run the ThresholdPressure object will be active, + but it will *not* be properly initialized with numerical values. The + values must instead come from the THPRES vector in the restart file. + */ + if (simConfig.getThresholdPressure().restart()) + return; + + // allocate the array which specifies the threshold pressures + thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0); + thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0); +} + +template +void EclGenericThresholdPressure:: +applyExplicitThresholdPressures_() +{ + const SimulationConfig& simConfig = eclState_.getSimulationConfig(); + const auto& thpres = simConfig.getThresholdPressure(); + + // set the threshold pressures for all EQUIL region boundaries which have a + // intersection in the grid + auto elemIt = gridView_.template begin(); + const auto& elemEndIt = gridView_.template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const auto& elem = *elemIt; + if (elem.partitionType() != Dune::InteriorEntity) + continue; + + auto isIt = gridView_.ibegin(elem); + const auto& isEndIt = gridView_.iend(elem); + for (; isIt != isEndIt; ++ isIt) { + // store intersection, this might be costly + const auto& intersection = *isIt; + + if (intersection.boundary()) + continue; // ignore boundary intersections for now (TODO?) + else if (!intersection.neighbor()) //processor boundary but not domain boundary + continue; + + const auto& inside = intersection.inside(); + const auto& outside = intersection.outside(); + + unsigned insideElemIdx = elementMapper_.index(inside); + unsigned outsideElemIdx = elementMapper_.index(outside); + + unsigned equilRegionInside = elemEquilRegion_[insideElemIdx]; + unsigned equilRegionOutside = elemEquilRegion_[outsideElemIdx]; + if (thpres.hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1)) { + Scalar pth = 0.0; + if (thpres.hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) { + // threshold pressure explicitly specified + pth = thpres.getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1); + } + else { + // take the threshold pressure from the initial condition + unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside; + pth = thpresDefault_[offset]; + } + + unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside; + unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside; + + thpres_[offset1] = pth; + thpres_[offset2] = pth; + } + } + } + + if (enableExperiments_) { + // apply threshold pressures accross faults (experimental!) + if (deck_.hasKeyword("THPRESFT")) + extractThpresft_(deck_.getKeyword("THPRESFT")); + } +} + +template +void EclGenericThresholdPressure:: +extractThpresft_(const DeckKeyword& thpresftKeyword) +{ + // retrieve the faults collection. + const FaultCollection& faults = eclState_.getFaults(); + + // extract the multipliers from the deck keyword + int numFaults = faults.size(); + int numCartesianElem = eclState_.getInputGrid().getCartesianSize(); + thpresftValues_.resize(numFaults, -1.0); + cartElemFaultIdx_.resize(numCartesianElem, -1); + for (size_t recordIdx = 0; recordIdx < thpresftKeyword.size(); ++ recordIdx) { + const DeckRecord& record = thpresftKeyword.getRecord(recordIdx); + + const std::string& faultName = record.getItem("FAULT_NAME").getTrimmedString(0); + Scalar thpresValue = record.getItem("VALUE").getSIDouble(0); + + for (size_t faultIdx = 0; faultIdx < faults.size(); faultIdx++) { + auto& fault = faults.getFault(faultIdx); + if (fault.getName() != faultName) + continue; + + thpresftValues_[faultIdx] = thpresValue; + for (const FaultFace& face: fault) + // "face" is a misnomer because the object describes a set of cell + // indices, but we go with the conventions of the parser here... + for (size_t cartElemIdx: face) + cartElemFaultIdx_[cartElemIdx] = faultIdx; + } + } +} + +#if HAVE_DUNE_FEM +template class EclGenericThresholdPressure>>, + Dune::MultipleCodimMultipleGeomTypeMapper>>, + double>; +#else +template class EclGenericThresholdPressure>, + Dune::MultipleCodimMultipleGeomTypeMapper>,Dune::Impl::MCMGFailLayout>, + double>; +#endif + +template class EclGenericThresholdPressure, + Dune::GridView>, + Dune::MultipleCodimMultipleGeomTypeMapper>, Dune::Impl::MCMGFailLayout>, + double>; + +} // namespace Opm diff --git a/ebos/eclgenericthresholdpressure.hh b/ebos/eclgenericthresholdpressure.hh new file mode 100644 index 000000000..bbc4779b2 --- /dev/null +++ b/ebos/eclgenericthresholdpressure.hh @@ -0,0 +1,111 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::EclThresholdPressure + */ +#ifndef EWOMS_ECL_GENERIC_THRESHOLD_PRESSURE_HH +#define EWOMS_ECL_GENERIC_THRESHOLD_PRESSURE_HH + +#include + +#include + +namespace Opm { + +class Deck; +class DeckKeyword; +class EclipseState; + +template +class EclGenericThresholdPressure { +public: + using CartesianIndexMapper = Dune::CartesianIndexMapper; + + EclGenericThresholdPressure(const CartesianIndexMapper& cartMapper, + const GridView& gridView, + const ElementMapper& elementMapper, + const EclipseState& eclState, + const Deck& deck, + bool enableExperiments); + + /*! + * \brief Returns the theshold pressure [Pa] for the intersection between two elements. + * + * This is tailor made for the E100 threshold pressure mechanism and it is thus quite + * a hack: First of all threshold pressures in general are unphysical, and second, + * they should be different for the fluid phase but are not. Anyway, this seems to be + * E100's way of doing things, so we do it the same way. + */ + Scalar thresholdPressure(int elem1Idx, int elem2Idx) const; + + /*! + * \brief Return the raw array with the threshold pressures + * + * This is used for the restart capability. + */ + const std::vector& data() const + { return thpres_; } + + /*! + * \brief Set the threshold pressures from a raw array + * + * This is used for the restart capability. + */ + void setFromRestart(const std::vector& values) + { thpres_ = values; } + +protected: + /*! + * \brief Actually compute the threshold pressures over a face as a pre-compute step. + */ + void finishInit(); + + // internalize the threshold pressures which where explicitly specified via the + // THPRES keyword. + void applyExplicitThresholdPressures_(); + + void extractThpresft_(const DeckKeyword& thpresftKeyword); + + const CartesianIndexMapper& cartMapper_; + const GridView& gridView_; + const ElementMapper& elementMapper_; + const EclipseState& eclState_; + const Deck& deck_; + std::vector thpresDefault_; + std::vector thpres_; + unsigned numEquilRegions_; + std::vector elemEquilRegion_; + + // threshold pressure accross faults. EXPERIMENTAL! + std::vector thpresftValues_; + std::vector cartElemFaultIdx_; + + bool enableThresholdPressure_; + bool enableExperiments_; +}; + +} // namespace Opm + +#endif diff --git a/ebos/eclthresholdpressure.hh b/ebos/eclthresholdpressure.hh index d20214c50..a3ef95e2f 100644 --- a/ebos/eclthresholdpressure.hh +++ b/ebos/eclthresholdpressure.hh @@ -29,23 +29,15 @@ #define EWOMS_ECL_THRESHOLD_PRESSURE_HH #include +#include +#include +#include #include #include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include +#include #include -#include namespace Opm { @@ -60,8 +52,15 @@ namespace Opm { * than the threshold pressure, it is reduced by the threshold pressure. */ template -class EclThresholdPressure +class EclThresholdPressure : public EclGenericThresholdPressure, + GetPropType, + GetPropType, + GetPropType> { + using BaseType = EclGenericThresholdPressure, + GetPropType, + GetPropType, + GetPropType>; using Simulator = GetPropType; using Scalar = GetPropType; using Evaluation = GetPropType; @@ -73,9 +72,14 @@ class EclThresholdPressure public: EclThresholdPressure(const Simulator& simulator) - : simulator_(simulator) + : BaseType(simulator.vanguard().cartesianIndexMapper(), + simulator.vanguard().gridView(), + simulator.model().elementMapper(), + simulator.vanguard().eclState(), + simulator.vanguard().deck(), + enableExperiments) + , simulator_(simulator) { - enableThresholdPressure_ = false; } /*! @@ -83,120 +87,13 @@ public: */ void finishInit() { - const auto& gridView = simulator_.gridView(); - - unsigned numElements = gridView.size(/*codim=*/0); - - // this code assumes that the DOFs are the elements. (i.e., an - // ECFV spatial discretization with TPFA). if you try to use - // it with something else, you're currently out of luck, - // sorry! - assert(simulator_.model().numGridDof() == numElements); - - const auto& vanguard = simulator_.vanguard(); - const auto& eclState = vanguard.eclState(); - const auto& simConfig = eclState.getSimulationConfig(); - - enableThresholdPressure_ = simConfig.useThresholdPressure(); - if (!enableThresholdPressure_) - return; - - numEquilRegions_ = eclState.getTableManager().getEqldims().getNumEquilRegions(); - if (numEquilRegions_ > 0xff) { - // make sure that the index of an equilibration region can be stored in a - // single byte - throw std::runtime_error("The maximum number of supported equilibration regions is 255!"); + this->BaseType::finishInit(); + if (this->enableThresholdPressure_ && !this->thpresDefault_.empty()) { + this->computeDefaultThresholdPressures_(); + this->applyExplicitThresholdPressures_(); } - - // internalize the data specified using the EQLNUM keyword - const auto& fp = eclState.fieldProps(); - const auto& equilRegionData = fp.get_int("EQLNUM"); - elemEquilRegion_.resize(numElements, 0); - for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { - elemEquilRegion_[elemIdx] = equilRegionData[elemIdx] - 1; - } - - /* - If this is a restart run the ThresholdPressure object will be active, - but it will *not* be properly initialized with numerical values. The - values must instead come from the THPRES vector in the restart file. - */ - if (simConfig.getThresholdPressure().restart()) - return; - - // allocate the array which specifies the threshold pressures - thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0); - thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0); - - computeDefaultThresholdPressures_(); - applyExplicitThresholdPressures_(); } - /*! - * \brief Returns the theshold pressure [Pa] for the intersection between two elements. - * - * This is tailor made for the E100 threshold pressure mechanism and it is thus quite - * a hack: First of all threshold pressures in general are unphysical, and second, - * they should be different for the fluid phase but are not. Anyway, this seems to be - * E100's way of doing things, so we do it the same way. - */ - Scalar thresholdPressure(int elem1Idx, int elem2Idx) const - { - if (!enableThresholdPressure_) - return 0.0; - - if (enableExperiments) { - // threshold pressure accross faults - if (!thpresftValues_.empty()) { - const auto& vanguard = simulator_.vanguard(); - int cartElem1Idx = vanguard.cartesianIndex(elem1Idx); - int cartElem2Idx = vanguard.cartesianIndex(elem2Idx); - - assert(0 <= cartElem1Idx && static_cast(cartElemFaultIdx_.size()) > cartElem1Idx); - assert(0 <= cartElem2Idx && static_cast(cartElemFaultIdx_.size()) > cartElem2Idx); - - int fault1Idx = cartElemFaultIdx_[cartElem1Idx]; - int fault2Idx = cartElemFaultIdx_[cartElem2Idx]; - if (fault1Idx != -1 && fault1Idx == fault2Idx) - // inside a fault there's no threshold pressure, even accross EQUIL - // regions. - return 0.0; - if (fault1Idx != fault2Idx) { - // TODO: which value if a cell is part of multiple faults? we take - // the maximum here. - Scalar val1 = (fault1Idx >= 0) ? thpresftValues_[fault1Idx] : 0.0; - Scalar val2 = (fault2Idx >= 0) ? thpresftValues_[fault2Idx] : 0.0; - return std::max(val1, val2); - } - } - } - - // threshold pressure accross EQUIL regions - unsigned short equilRegion1Idx = elemEquilRegion_[elem1Idx]; - unsigned short equilRegion2Idx = elemEquilRegion_[elem2Idx]; - - if (equilRegion1Idx == equilRegion2Idx) - return 0.0; - - return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx]; - } - - /*! - * \brief Return the raw array with the threshold pressures - * - * This is used for the restart capability. - */ - const std::vector& data() const - { return thpres_; } - - /*! - * \brief Set the threshold pressures from a raw array - * - * This is used for the restart capability. - */ - void setFromRestart(const std::vector& values) - { thpres_ = values; } - private: // compute the defaults of the threshold pressures using the initial condition void computeDefaultThresholdPressures_() @@ -228,8 +125,8 @@ private: unsigned insideElemIdx = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0); unsigned outsideElemIdx = elemCtx.globalSpaceIndex(j, /*timeIdx=*/0); - unsigned equilRegionInside = elemEquilRegion_[insideElemIdx]; - unsigned equilRegionOutside = elemEquilRegion_[outsideElemIdx]; + unsigned equilRegionInside = this->elemEquilRegion_[insideElemIdx]; + unsigned equilRegionOutside = this->elemEquilRegion_[outsideElemIdx]; if (equilRegionInside == equilRegionOutside) // the current face is not at the boundary between EQUIL regions! @@ -255,133 +152,21 @@ private: } } - int offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside; - int offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside; + int offset1 = equilRegionInside*this->numEquilRegions_ + equilRegionOutside; + int offset2 = equilRegionOutside*this->numEquilRegions_ + equilRegionInside; - thpresDefault_[offset1] = std::max(thpresDefault_[offset1], pth); - thpresDefault_[offset2] = std::max(thpresDefault_[offset2], pth); + this->thpresDefault_[offset1] = std::max(this->thpresDefault_[offset1], pth); + this->thpresDefault_[offset2] = std::max(this->thpresDefault_[offset2], pth); } } // make sure that the threshold pressures is consistent for parallel // runs. (i.e. take the maximum of all processes) - for (unsigned i = 0; i < thpresDefault_.size(); ++i) - thpresDefault_[i] = gridView.comm().max(thpresDefault_[i]); - } - - // internalize the threshold pressures which where explicitly specified via the - // THPRES keyword. - void applyExplicitThresholdPressures_() - { - const auto& vanguard = simulator_.vanguard(); - const auto& gridView = vanguard.gridView(); - const auto& elementMapper = simulator_.model().elementMapper(); - const auto& eclState = simulator_.vanguard().eclState(); - const SimulationConfig& simConfig = eclState.getSimulationConfig(); - const auto& thpres = simConfig.getThresholdPressure(); - - // set the threshold pressures for all EQUIL region boundaries which have a - // intersection in the grid - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - for (; elemIt != elemEndIt; ++elemIt) { - const auto& elem = *elemIt; - if (elem.partitionType() != Dune::InteriorEntity) - continue; - - auto isIt = gridView.ibegin(elem); - const auto& isEndIt = gridView.iend(elem); - for (; isIt != isEndIt; ++ isIt) { - // store intersection, this might be costly - const auto& intersection = *isIt; - - if (intersection.boundary()) - continue; // ignore boundary intersections for now (TODO?) - else if (!intersection.neighbor()) //processor boundary but not domain boundary - continue; - - const auto& inside = intersection.inside(); - const auto& outside = intersection.outside(); - - unsigned insideElemIdx = elementMapper.index(inside); - unsigned outsideElemIdx = elementMapper.index(outside); - - unsigned equilRegionInside = elemEquilRegion_[insideElemIdx]; - unsigned equilRegionOutside = elemEquilRegion_[outsideElemIdx]; - if (thpres.hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1)) { - Scalar pth = 0.0; - if (thpres.hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) { - // threshold pressure explicitly specified - pth = thpres.getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1); - } - else { - // take the threshold pressure from the initial condition - unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside; - pth = thpresDefault_[offset]; - } - - unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside; - unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside; - - thpres_[offset1] = pth; - thpres_[offset2] = pth; - } - } - } - - if (enableExperiments) { - // apply threshold pressures accross faults (experimental!) - const auto& deck = simulator_.vanguard().deck(); - if (deck.hasKeyword("THPRESFT")) - extractThpresft_(deck.getKeyword("THPRESFT")); - } - - } - - void extractThpresft_(const DeckKeyword& thpresftKeyword) - { - // retrieve the faults collection. - const EclipseState& eclState = simulator_.vanguard().eclState(); - const FaultCollection& faults = eclState.getFaults(); - - // extract the multipliers from the deck keyword - int numFaults = faults.size(); - int numCartesianElem = eclState.getInputGrid().getCartesianSize(); - thpresftValues_.resize(numFaults, -1.0); - cartElemFaultIdx_.resize(numCartesianElem, -1); - for (size_t recordIdx = 0; recordIdx < thpresftKeyword.size(); ++ recordIdx) { - const DeckRecord& record = thpresftKeyword.getRecord(recordIdx); - - const std::string& faultName = record.getItem("FAULT_NAME").getTrimmedString(0); - Scalar thpresValue = record.getItem("VALUE").getSIDouble(0); - - for (size_t faultIdx = 0; faultIdx < faults.size(); faultIdx++) { - auto& fault = faults.getFault(faultIdx); - if (fault.getName() != faultName) - continue; - - thpresftValues_[faultIdx] = thpresValue; - for (const FaultFace& face: fault) - // "face" is a misnomer because the object describes a set of cell - // indices, but we go with the conventions of the parser here... - for (size_t cartElemIdx: face) - cartElemFaultIdx_[cartElemIdx] = faultIdx; - } - } + for (unsigned i = 0; i < this->thpresDefault_.size(); ++i) + this->thpresDefault_[i] = gridView.comm().max(this->thpresDefault_[i]); } const Simulator& simulator_; - - std::vector thpresDefault_; - std::vector thpres_; - unsigned numEquilRegions_; - std::vector elemEquilRegion_; - - // threshold pressure accross faults. EXPERIMENTAL! - std::vector thpresftValues_; - std::vector cartElemFaultIdx_; - - bool enableThresholdPressure_; }; } // namespace Opm