From 7a7d6d868dc3658154697b3fc937e956f32ac8b6 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 18 Jan 2017 16:46:53 +0100 Subject: [PATCH] ebos: fix the interactions between SWATINIT and threshold pressures 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. --- ebos/eclequilinitializer.hh | 34 ++++++--------- ebos/eclproblem.hh | 85 +++++++++++++++++++++++++++++-------- 2 files changed, 81 insertions(+), 38 deletions(-) diff --git a/ebos/eclequilinitializer.hh b/ebos/eclequilinitializer.hh index 107569978..de0441731 100644 --- a/ebos/eclequilinitializer.hh +++ b/ebos/eclequilinitializer.hh @@ -31,6 +31,7 @@ #include #include +#include // 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 + template 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 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(); } /*! diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index ec65f6522..b1c8f25b0 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -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(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 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(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 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 pffDofData_; - - std::vector maxPcnw_; }; } // namespace Ewoms