OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-23 01:36:25 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #241 from totto82/removeInitDupl

Browse Source 
Prepare for using the initial solution from ebos directly
This commit is contained in:
Atgeirr Flø Rasmussen 2017-12-05 08:44:29 +01:00 committed by GitHub
commit 6a5c6ea227
2 changed files with 17 additions and 134 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
ebos/eclequilinitializer.hh
View File

@ -48,7 +48,6 @@ NEW_PROP_TAG(FluidSystem);
NEW_PROP_TAG(GridView);
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(MaterialLaw);
NEW_PROP_TAG(EnableSwatinit);
}
/*!
@ -87,54 +86,19 @@ class EclEquilInitializer
public:
template <class EclMaterialLawManager>
EclEquilInitializer(const Simulator& simulator,
EclMaterialLawManager& materialLawManager,
bool enableSwatinit)
EclMaterialLawManager& materialLawManager)
: simulator_(simulator)
{
const auto& gridManager = simulator.gridManager();
const auto& deck = gridManager.deck();
const auto& eclState = gridManager.eclState();
const auto& equilGrid = gridManager.equilGrid();
unsigned numElems = gridManager.grid().size(0);
unsigned numEquilElems = gridManager.equilGrid().size(0);
unsigned numCartesianElems = gridManager.cartesianSize();
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef Opm::ThreePhaseMaterialTraits<double,
/*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> EquilTraits;
// create a separate instance of the material law manager just because opm-core
// only supports double as the type for scalars (but ebos may use float or quad)
std::vector<int> compressedToCartesianEquilElemIdx(numEquilElems);
std::vector<int> equilCartesianToCompressed( gridManager.equilCartesianSize(), -1 );
for (unsigned equilElemIdx = 0; equilElemIdx < numEquilElems; ++equilElemIdx)
{
unsigned int equilCartesianIdx = gridManager.equilCartesianIndex(equilElemIdx);
compressedToCartesianEquilElemIdx[equilElemIdx] = equilCartesianIdx;
equilCartesianToCompressed[ equilCartesianIdx ] = equilElemIdx;
}
Opm::EclMaterialLawManager<EquilTraits> equilMaterialLawManager =
Opm::EclMaterialLawManager<EquilTraits>();
equilMaterialLawManager.initFromDeck(deck, eclState, compressedToCartesianEquilElemIdx);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> initialState(equilMaterialLawManager,
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> initialState(materialLawManager,
gridManager.eclState(),
equilGrid,
simulator.problem().gravity()[dimWorld - 1],
enableSwatinit);
std::vector<int> localToEquilIndex( numElems, -1 );
for( unsigned int elemIdx = 0; elemIdx < numElems; ++elemIdx )
{
const int cartesianIndex = gridManager.cartesianIndex( elemIdx );
assert( equilCartesianToCompressed[ cartesianIndex ] >= 0 );
localToEquilIndex[ elemIdx ] = equilCartesianToCompressed[ cartesianIndex ];
}
gridManager.grid(),
simulator.problem().gravity()[dimWorld - 1]);
// copy the result into the array of initial fluid states
initialFluidStates_.resize(numCartesianElems);
@ -142,8 +106,6 @@ public:
unsigned cartesianElemIdx = gridManager.cartesianIndex(elemIdx);
auto& fluidState = initialFluidStates_[cartesianElemIdx];
const unsigned int equilElemIdx = localToEquilIndex[ elemIdx ];
// get the PVT region index of the current element
unsigned regionIdx = simulator_.problem().pvtRegionIndex(elemIdx);
@ -153,7 +115,7 @@ public:
if (!FluidSystem::phaseIsActive(phaseIdx))
S = 0.0;
else {
S = initialState.saturation()[phaseIdx][equilElemIdx];
S = initialState.saturation()[phaseIdx][elemIdx];
}
fluidState.setSaturation(phaseIdx, S);
}
@ -164,7 +126,7 @@ public:
// set the phase pressures.
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
fluidState.setPressure(phaseIdx, initialState.press()[phaseIdx][equilElemIdx]);
fluidState.setPressure(phaseIdx, initialState.press()[phaseIdx][elemIdx]);
// reset the phase compositions
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
@ -178,7 +140,7 @@ public:
if (FluidSystem::enableDissolvedGas()) {
// for gas and oil we have to translate surface volumes to mole fractions
// before we can set the composition in the fluid state
Scalar Rs = initialState.rs()[equilElemIdx];
Scalar Rs = initialState.rs()[elemIdx];
Scalar RsSat = FluidSystem::saturatedDissolutionFactor(fluidState, oilPhaseIdx, regionIdx);
if (Rs > RsSat)
@ -194,7 +156,7 @@ public:
// retrieve the surface volume of vaporized gas
if (FluidSystem::enableVaporizedOil()) {
Scalar Rv = initialState.rv()[equilElemIdx];
Scalar Rv = initialState.rv()[elemIdx];
Scalar RvSat = FluidSystem::saturatedDissolutionFactor(fluidState, gasPhaseIdx, regionIdx);
if (Rv > RvSat)
@ -207,16 +169,6 @@ public:
fluidState.setMoleFraction(gasPhaseIdx, oilCompIdx, xgO);
fluidState.setMoleFraction(gasPhaseIdx, gasCompIdx, 1 - xgO);
}
// 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());
}
}
}

87
ebos/eclproblem.hh
View File

@ -119,9 +119,6 @@ NEW_PROP_TAG(DisableWells);
// print statements in debug mode.
NEW_PROP_TAG(EnableDebuggingChecks);
// If this property is set to false, the SWATINIT keyword will not be handled by ebos.
NEW_PROP_TAG(EnableSwatinit);
// Set the problem property
SET_TYPE_PROP(EclBaseProblem, Problem, Ewoms::EclProblem<TypeTag>);
@ -238,8 +235,6 @@ SET_BOOL_PROP(EclBaseProblem, DisableWells, false);
// By default, we enable the debugging checks if we're compiled in debug mode
SET_BOOL_PROP(EclBaseProblem, EnableDebuggingChecks, true);
// ebos handles the SWATINIT keyword by default
SET_BOOL_PROP(EclBaseProblem, EnableSwatinit, true);
} // namespace Properties
/*!
@ -556,6 +551,10 @@ public:
// write the summary information after each time step
summaryWriter_.write(wellManager_);
}
// we no longer need the initial soluiton
if (this->simulator().episodeIndex() == 0)
initialFluidStates_.clear();
}
/*!
@ -972,26 +971,8 @@ public:
// 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 = eclState.getTableManager().getEqldims().getNumEquilRegions();
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();
}
/*!
@ -1026,48 +1007,9 @@ public:
const EclWellManager<TypeTag>& wellManager() const
{ return wellManager_; }
/*!
* \brief Apply the necessary measures mandated by the SWATINIT keyword the the
* material law parameters.
*
* If SWATINIT is not used or if the method has already been called, this method is a
* no-op.
*/
void applySwatinit()
{
const auto& deck = this->simulator().gridManager().deck();
const auto& eclState = this->simulator().gridManager().eclState();
const auto& deckProps = eclState.get3DProperties();
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;
}
// the solution (most likely) has changed, so we need to invalidate the cache for
// intensive quantities
this->simulator().model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
// temporary solution to facilitate output of initial state from flow
const ScalarFluidState& initialFluidState(unsigned globalDofIdx ) const {
return initialFluidStates_[globalDofIdx];
}
private:
@ -1275,20 +1217,9 @@ private:
void readEquilInitialCondition_()
{
const auto& deck = this->simulator().gridManager().deck();
const auto& eclState = this->simulator().gridManager().eclState();
int numEquilRegions = eclState.getTableManager().getEqldims().getNumEquilRegions();
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.
// initial condition corresponds to hydrostatic conditions.
typedef Ewoms::EclEquilInitializer<TypeTag> EquilInitializer;
EquilInitializer equilInitializer(this->simulator(),
*materialLawManager_,
/*enableSwatinit=*/useThpres && useSwatinit);
EquilInitializer equilInitializer(this->simulator(), *materialLawManager_);
// since the EquilInitializer provides fluid states that are consistent with the
// black-oil model, we can use naive instead of mass conservative determination