OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

eclthresholdpressure: split in typetag dependent and typetag-independent parts

Browse Source 
this allows using explicit template instantation to only
compile this code per grid, not per simulator object
This commit is contained in:
Arne Morten Kvarving 2021-05-05 12:13:25 +02:00
parent af68511a63
commit f5fa602abb
4 changed files with 408 additions and 246 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -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

265
ebos/eclgenericthresholdpressure.cc Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <config.h>
#include <ebos/eclgenericthresholdpressure.hh>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/Eqldims.hpp>
#include <opm/parser/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
#include <opm/parser/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
#include <dune/grid/common/mcmgmapper.hh>
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/polyhedralgrid.hh>
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif
#include <algorithm>
#include <cassert>
#include <stdexcept>
namespace Opm {
template<class Grid, class GridView, class ElementMapper, class Scalar>
EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
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<class Grid, class GridView, class ElementMapper,class Scalar>
Scalar EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
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<int>(cartElemFaultIdx_.size()) > cartElem1Idx);
assert(0 <= cartElem2Idx && static_cast<int>(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<class Grid, class GridView, class ElementMapper, class Scalar>
void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
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<class Grid, class GridView, class ElementMapper, class Scalar>
void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
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</*codim=*/ 0>();
const auto& elemEndIt = gridView_.template end</*codim=*/ 0>();
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<class Grid, class GridView, class ElementMapper, class Scalar>
void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
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::CpGrid,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
double>;
#else
template class EclGenericThresholdPressure<Dune::CpGrid,
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,Dune::Impl::MCMGFailLayout>,
double>;
#endif
template class EclGenericThresholdPressure<Dune::PolyhedralGrid<3,3,double>,
Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>, Dune::Impl::MCMGFailLayout>,
double>;
} // namespace Opm

111
ebos/eclgenericthresholdpressure.hh Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <opm/grid/common/CartesianIndexMapper.hpp>
#include <vector>
namespace Opm {
class Deck;
class DeckKeyword;
class EclipseState;
template<class Grid, class GridView, class ElementMapper, class Scalar>
class EclGenericThresholdPressure {
public:
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
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<Scalar>& 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<Scalar>& 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<Scalar> thpresDefault_;
std::vector<Scalar> thpres_;
unsigned numEquilRegions_;
std::vector<unsigned char> elemEquilRegion_;
// threshold pressure accross faults. EXPERIMENTAL!
std::vector<Scalar> thpresftValues_;
std::vector<int> cartElemFaultIdx_;
bool enableThresholdPressure_;
bool enableExperiments_;
};
} // namespace Opm
#endif

277
ebos/eclthresholdpressure.hh
View File

@ -29,23 +29,15 @@
#define EWOMS_ECL_THRESHOLD_PRESSURE_HH
#include <opm/models/utils/propertysystem.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/common/multiphasebaseproperties.hh>
#include <ebos/eclgenericthresholdpressure.hh>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/Eqldims.hpp>
#include <opm/parser/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
#include <opm/parser/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
#include <dune/grid/common/gridenums.hh>
#include <dune/common/version.hh>
#include <array>
#include <algorithm>
#include <vector>
#include <unordered_map>
namespace Opm {
@ -60,8 +52,15 @@ namespace Opm {
* than the threshold pressure, it is reduced by the threshold pressure.
*/
template <class TypeTag>
class EclThresholdPressure
class EclThresholdPressure : public EclGenericThresholdPressure<GetPropType<TypeTag, Properties::Grid>,
GetPropType<TypeTag, Properties::GridView>,
GetPropType<TypeTag, Properties::ElementMapper>,
GetPropType<TypeTag, Properties::Scalar>>
{
using BaseType = EclGenericThresholdPressure<GetPropType<TypeTag, Properties::Grid>,
GetPropType<TypeTag, Properties::GridView>,
GetPropType<TypeTag, Properties::ElementMapper>,
GetPropType<TypeTag, Properties::Scalar>>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
@ -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<int>(cartElemFaultIdx_.size()) > cartElem1Idx);
assert(0 <= cartElem2Idx && static_cast<int>(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<Scalar>& 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<Scalar>& 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</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
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<Scalar> thpresDefault_;
std::vector<Scalar> thpres_;
unsigned numEquilRegions_;
std::vector<unsigned char> elemEquilRegion_;
// threshold pressure accross faults. EXPERIMENTAL!
std::vector<Scalar> thpresftValues_;
std::vector<int> cartElemFaultIdx_;
bool enableThresholdPressure_;
};
} // namespace Opm