diff --git a/ebos/eclequilinitializer.hh b/ebos/eclequilinitializer.hh index 39b070684..eb5ac26a9 100644 --- a/ebos/eclequilinitializer.hh +++ b/ebos/eclequilinitializer.hh @@ -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 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 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 compressedToCartesianEquilElemIdx(numEquilElems); - std::vector 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 equilMaterialLawManager = - Opm::EclMaterialLawManager(); - equilMaterialLawManager.initFromDeck(deck, eclState, compressedToCartesianEquilElemIdx); - - Opm::EQUIL::DeckDependent::InitialStateComputer initialState(equilMaterialLawManager, + Opm::EQUIL::DeckDependent::InitialStateComputer initialState(materialLawManager, gridManager.eclState(), - equilGrid, - simulator.problem().gravity()[dimWorld - 1], - enableSwatinit); - - - std::vector 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()); - } } } diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 50ccadec0..f2af4ed41 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -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); @@ -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& 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 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 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