OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-27 03:30:17 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1774 from andlaus/support_THPRESFT

Browse Source 
Support THPRESFT
This commit is contained in:
Tor Harald Sandve 2019-04-03 09:20:53 +02:00 committed by GitHub
commit ae6edf6a06
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 84 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
ebos/eclproblem.hh
View File

@ -1691,6 +1691,7 @@ private:
void checkDeckCompatibility_() const
{
const auto& deck = this->simulator().vanguard().deck();
bool beVerbose = this->simulator().gridView().comm().rank() == 0;
if (enableApiTracking)
throw std::logic_error("API tracking is not yet implemented but requested at compile time.");
@ -1716,17 +1717,17 @@ private:
else if (!enableEnergy && deckEnergyEnabled)
throw std::runtime_error("The deck enables the TEMP or the THERMAL option, but the simulator is not compiled to support either.");
if (deckEnergyEnabled && deck.hasKeyword("TEMP"))
std::cout << "WARNING: The deck requests the TEMP option, i.e., treating energy "
if (deckEnergyEnabled && deck.hasKeyword("TEMP") && beVerbose)
std::cerr << "WARNING: The deck requests the TEMP option, i.e., treating energy "
<< "conservation as a post processing step. This is currently unsupported, "
<< "i.e., energy conservation is always handled fully implicitly.";
<< "i.e., energy conservation is always handled fully implicitly." << std::endl;
int numDeckPhases = FluidSystem::numActivePhases();
if (numDeckPhases < Indices::numPhases)
std::cout << "WARNING: The number of active phases specified by the deck ("
if (numDeckPhases < Indices::numPhases && beVerbose)
std::cerr << "WARNING: The number of active phases specified by the deck ("
<< numDeckPhases << ") is smaller than the number of compiled-in phases ("
<< Indices::numPhases << "). This usually results in a significant "
<< "performance degradation compared to using a specialized simulator.";
<< "performance degradation compared to using a specialized simulator." << std::endl;
else if (numDeckPhases < Indices::numPhases)
throw std::runtime_error("The deck enables "+std::to_string(numDeckPhases)+" phases "
"while this simulator can only handle "+
@ -1739,7 +1740,7 @@ private:
throw std::runtime_error("The deck enables gas, but this simulator cannot handle it.");
if (FluidSystem::phaseIsActive(waterPhaseIdx) && !Indices::waterEnabled)
throw std::runtime_error("The deck enables water, but this simulator cannot handle it.");
// the opposite cases should be fine (albeit a bit slower than possible)
// the opposite cases should be fine (albeit a bit slower than what's possible)
}
bool drsdtActive_() const

79
ebos/eclthresholdpressure.hh
View File

@ -56,6 +56,7 @@ NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(Evaluation);
NEW_PROP_TAG(ElementContext);
NEW_PROP_TAG(FluidSystem);
NEW_PROP_TAG(EnableExperiments);
END_PROPERTIES
@ -80,6 +81,7 @@ class EclThresholdPressure
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
enum { enableExperiments = GET_PROP_VALUE(TypeTag, EnableExperiments) };
enum { numPhases = FluidSystem::numPhases };
public:
@ -158,13 +160,40 @@ public:
* 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 elemIdx1, int elemIdx2) const
Scalar thresholdPressure(int elem1Idx, int elem2Idx) const
{
if (!enableThresholdPressure_)
return 0.0;
unsigned short equilRegion1Idx = elemEquilRegion_[elemIdx1];
unsigned short equilRegion2Idx = elemEquilRegion_[elemIdx2];
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 && cartElemFaultIdx_.size() > cartElem1Idx);
assert(0 <= cartElem2Idx && 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;
@ -257,6 +286,7 @@ private:
const auto& gridView = vanguard.gridView();
const auto& elementMapper = simulator_.model().elementMapper();
const auto& eclState = simulator_.vanguard().eclState();
const auto& deck = simulator_.vanguard().deck();
const Opm::SimulationConfig& simConfig = eclState.getSimulationConfig();
const auto& thpres = simConfig.getThresholdPressure();
@ -307,6 +337,45 @@ private:
}
}
}
if (enableExperiments) {
// apply threshold pressures accross faults (experimental!)
if (deck.hasKeyword("THPRESFT"))
extractThpresft_(deck.getKeyword("THPRESFT"));
}
}
void extractThpresft_(const Opm::DeckKeyword& thpresftKeyword)
{
// retrieve the faults collection.
const Opm::EclipseState& eclState = simulator_.vanguard().eclState();
const Opm::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 Opm::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 Opm::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;
}
}
}
const Simulator& simulator_;
@ -316,6 +385,10 @@ private:
unsigned numEquilRegions_;
std::vector<unsigned char> elemEquilRegion_;
// threshold pressure accross faults. EXPERIMENTAL!
std::vector<Scalar> thpresftValues_;
std::vector<int> cartElemFaultIdx_;
bool enableThresholdPressure_;
};