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

ebos: fix the interactions between SWATINIT and threshold pressures

Browse Source 
hopefully this makes standalone `ebos` arrive at the same initial
condition as `flow_ebos` if both, SWATINIT and threshold pressures are
enabled. we need to calculate the initial condition twice either
threshold pressures and SWATINIT are enabled. (`ebos` and `flow_ebos`
diverged after OPM/opm-core#1129.)

the proposed patch is a kludge IMO, but in the light that in my
opinion, SWATINIT and threshold pressures are both physically not
justified and given the fact that SWATINIT must not be considered for
the threshold pressues should be considered to be a bug of the
reference simulator, I think the patch is okay.
This commit is contained in:
Andreas Lauser 2017-01-18 16:46:53 +01:00
parent 41d258f1c8
commit 7a7d6d868d
2 changed files with 81 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
ebos/eclequilinitializer.hh
View File

@ -31,6 +31,7 @@
#include <ewoms/common/propertysystem.hh>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
// the ordering of these includes matters. do not touch it if you're not prepared to deal
// with some trouble!
@ -85,9 +86,10 @@ class EclEquilInitializer
enum { dimWorld = GridView::dimensionworld };
public:
template <class MaxPcnwVector>
template <class EclMaterialLawManager>
EclEquilInitializer(const Simulator& simulator,
MaxPcnwVector& maxPcnw)
EclMaterialLawManager& materialLawManager,
bool enableSwatinit)
: simulator_(simulator)
{
const auto& gridManager = simulator.gridManager();
@ -137,9 +139,6 @@ public:
/*numFaces=*/0, // we don't care here
opmPhaseUsage.num_phases);
// tell the initializers whether SWATINIT should be applied or not.
bool enableSwatInit = GET_PROP_VALUE(TypeTag, EnableSwatinit);
// do the actual computation.
Opm::initStateEquil(equilGrid,
opmBlackoilProps,
@ -147,7 +146,7 @@ public:
gridManager.eclState(),
simulator.problem().gravity()[dimWorld - 1],
opmBlackoilState,
enableSwatInit);
enableSwatinit);
std::vector<int> localToEquilIndex( numElems, -1 );
for( unsigned int elemIdx = 0; elemIdx < numElems; ++elemIdx )
@ -157,10 +156,6 @@ public:
localToEquilIndex[ elemIdx ] = equilCartesianToCompressed[ cartesianIndex ];
}
// resize the array of maximum water-oil capillary pressures (needed to apply the
// SWATINIT keyword).
maxPcnw.resize(numElems);
// copy the result into the array of initial fluid states
initialFluidStates_.resize(numCartesianElems);
for (unsigned int elemIdx = 0; elemIdx < numElems; ++elemIdx) {
@ -257,17 +252,16 @@ public:
fluidState.setMoleFraction(gasPhaseIdx, gasCompIdx, 1 - xgO);
}
// store the maximum oil-water capillary pressure for later (required by the
// SWATINIT keyword)
const auto& equilScalingPoints =
equilMaterialLawManager->oilWaterScaledEpsPointsDrainage(equilElemIdx);
maxPcnw[elemIdx] = equilScalingPoints.maxPcnw();
// deal with the changed pressure scaling due to SWATINIT if SWATINIT is
// requested to be applied. this is quite hacky but hey it works!
if (enableSwatinit) {
const auto& equilScalingPoints =
equilMaterialLawManager->oilWaterScaledEpsPointsDrainage(equilElemIdx);
auto& scalingPoints =
materialLawManager.oilWaterScaledEpsPointsDrainage(elemIdx);
scalingPoints.setMaxPcnw(equilScalingPoints.maxPcnw());
}
}
// immediately forget about the SWATINIT stuff if the SWATINIT keyword is not
// present in the deck
if (!deck.hasKeyword("SWATINIT"))
maxPcnw.clear();
}
/*!

85
ebos/eclproblem.hh
View File

@ -817,6 +817,28 @@ public:
// this point, because determining the threshold pressures may require to access
// the initial solution.
thresholdPressures_.finishInit();
// apply SWATINIT if requested by the programmer. this is only necessary if
// SWATINIT has not yet been considered at the first time the EquilInitializer
// was used, i.e., only if threshold pressures are enabled in addition to
// SWATINIT.
const auto& deck = this->simulator().gridManager().deck();
const auto& eclState = this->simulator().gridManager().eclState();
int numEquilRegions =
deck.getKeyword("EQLDIMS").getRecord(0).getItem("NTEQUL").template get<int>(0);
bool useThpres = deck.hasKeyword("THPRES") && numEquilRegions > 1;
bool useSwatinit =
GET_PROP_VALUE(TypeTag, EnableSwatinit) &&
eclState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT");
if (useThpres && useSwatinit)
applySwatinit();
// release the memory of the EQUIL grid since it's no longer needed after this point
this->simulator().gridManager().releaseEquilGrid();
// the fluid states specifying the initial condition are also no longer required.
initialFluidStates_.clear();
}
/*!
@ -860,16 +882,39 @@ public:
*/
void applySwatinit()
{
if (maxPcnw_.empty())
return; // SWATINIT not applicable
const auto& deck = this->simulator().gridManager().deck();
const auto& eclState = this->simulator().gridManager().eclState();
const auto& deckProps = eclState.get3DProperties();
unsigned numElems = this->simulator().gridView().size(/*codim=*/0);
for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
auto& scalingPoints = materialLawManager_->oilWaterScaledEpsPointsDrainage(elemIdx);
scalingPoints.setMaxPcnw(maxPcnw_[elemIdx]);
if (!deckProps.hasDeckDoubleGridProperty("SWATINIT"))
return; // SWATINIT is not in the deck
else if (!deck.hasKeyword("EQUIL"))
// SWATINIT only applies if the initial solution is specified using the EQUIL
// keyword.
return;
// SWATINIT applies. We have to do a complete re-initialization here. this is a
// kludge, but to calculate the threshold pressures the initial solution without
// SWATINIT is required!
typedef Ewoms::EclEquilInitializer<TypeTag> EquilInitializer;
EquilInitializer equilInitializer(this->simulator(),
*materialLawManager_,
/*enableSwatinit=*/true);
auto& model = this->model();
size_t numElems = model.numGridDof();
for (size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
const auto& fs = equilInitializer.initialFluidState(elemIdx);
auto& priVars0 = model.solution(/*timeIdx=*/0)[elemIdx];
auto& priVars1 = model.solution(/*timeIdx=*/1)[elemIdx];
priVars1.assignNaive(fs);
priVars0 = priVars1;
}
maxPcnw_.clear();
// the solution (most likely) has changed, so we need to invalidate the cache for
// intensive quantities
this->simulator().model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
}
private:
@ -1050,20 +1095,26 @@ private:
readExplicitInitialCondition_();
else
readEquilInitialCondition_();
// release the memory of the EQUIL grid since it's no longer needed after this point
this->simulator().gridManager().releaseEquilGrid();
// apply SWATINIT if requested by the problem and if it is enabled by the deck
if (GET_PROP_VALUE(TypeTag, EnableSwatinit))
applySwatinit();
}
void readEquilInitialCondition_()
{
// initial condition corresponds to hydrostatic conditions
const auto& deck = this->simulator().gridManager().deck();
const auto& eclState = this->simulator().gridManager().eclState();
int numEquilRegions =
deck.getKeyword("EQLDIMS").getRecord(0).getItem("NTEQUL").template get<int>(0);
bool useThpres = deck.hasKeyword("THPRES") && numEquilRegions > 1;
bool useSwatinit =
GET_PROP_VALUE(TypeTag, EnableSwatinit) &&
eclState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT");
// initial condition corresponds to hydrostatic conditions. The SWATINIT keyword
// can be considered here directly if threshold pressures are disabled.
typedef Ewoms::EclEquilInitializer<TypeTag> EquilInitializer;
EquilInitializer equilInitializer(this->simulator(), maxPcnw_);
EquilInitializer equilInitializer(this->simulator(),
*materialLawManager_,
/*enableSwatinit=*/useThpres && useSwatinit);
// since the EquilInitializer provides fluid states that are consistent with the
// black-oil model, we can use naive instead of mass conservative determination
@ -1365,8 +1416,6 @@ private:
EclSummaryWriter summaryWriter_;
PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
std::vector<Scalar> maxPcnw_;
};
} // namespace Ewoms