From 427741fe840abc7872fe4b01c71f162609717c40 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 13 Feb 2019 12:54:52 +0100 Subject: [PATCH 1/3] ebos: Fix restart from ECL files Some time loop stuff was missing in the doobly-doo, the init() method of the well model was not called and there was the slightly deeper issue that the initial solutions where not calculated on restarts which breaks everything that relies on them. (at the moment, that's everything which is related to non-trivial boundary contitions.) --- ebos/eclproblem.hh | 80 +++++++++++++++++++++++++++++++++------------- ebos/eclwriter.hh | 10 +++++- 2 files changed, 67 insertions(+), 23 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index ab624a6b3..fceaea6f6 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -578,28 +578,18 @@ public: const auto& initconfig = eclState.getInitConfig(); const auto& timeMap = simulator.vanguard().schedule().getTimeMap(); - if(initconfig.restartRequested()) { - // Set the start time of the simulation - simulator.setStartTime( timeMap.getStartTime(/*timeStepIdx=*/initconfig.getRestartStep()) ); - simulator.setEpisodeIndex(initconfig.getRestartStep()); - simulator.setEpisodeLength(0.0); - simulator.setTimeStepSize(0.0); - readEclRestartSolution_(); - } - else { - readInitialCondition_(); - // Set the start time of the simulation - simulator.setStartTime( timeMap.getStartTime(/*timeStepIdx=*/0) ); + readInitialCondition_(); + // Set the start time of the simulation + simulator.setStartTime(timeMap.getStartTime(/*timeStepIdx=*/0)); - // We want the episode index to be the same as the report step index to make - // things simpler, so we have to set the episode index to -1 because it is - // incremented inside beginEpisode(). The size of the initial time step and - // length of the initial episode is set to zero for the same reason. - simulator.setEpisodeIndex(-1); - simulator.setEpisodeLength(0.0); - simulator.setTimeStepSize(0.0); - } + // We want the episode index to be the same as the report step index to make + // things simpler, so we have to set the episode index to -1 because it is + // incremented inside beginEpisode(). The size of the initial time step and + // length of the initial episode is set to zero for the same reason. + simulator.setEpisodeIndex(-1); + simulator.setEpisodeLength(0.0); + simulator.setTimeStepSize(0.0); updatePffDofData_(); @@ -1350,10 +1340,12 @@ public: */ void initialSolutionApplied() { + const auto& eclState = this->simulator().vanguard().eclState(); + // initialize the wells. Note that this needs to be done after initializing the // intrinsic permeabilities and the after applying the initial solution because // the well model uses these... - wellModel_.init(this->simulator().vanguard().eclState(), this->simulator().vanguard().schedule()); + wellModel_.init(eclState, this->simulator().vanguard().schedule()); // let the object for threshold pressures initialize itself. this is done only at // this point, because determining the threshold pressures may require to access @@ -1366,6 +1358,10 @@ public: updateCompositionChangeLimits_(); aquiferModel_.initialSolutionApplied(); + + const auto& initconfig = eclState.getInitConfig(); + if (initconfig.restartRequested()) + readEclRestartSolution_(); } /*! @@ -1837,7 +1833,19 @@ private: void readEclRestartSolution_() { - eclWriter_->restartBegin(); + // Set the start time of the simulation + const auto& eclState = this->simulator().vanguard().eclState(); + const auto& timeMap = this->simulator().vanguard().schedule().getTimeMap(); + const auto& initconfig = eclState.getInitConfig(); + int episodeIdx = initconfig.getRestartStep() - 1; + + this->simulator().setStartTime(timeMap.getStartTime(/*timeStepIdx=*/0)); + this->simulator().setTime(timeMap.getTimePassedUntil(episodeIdx)); + this->simulator().setEpisodeIndex(episodeIdx); + this->simulator().setEpisodeLength(timeMap.getTimeStepLength(episodeIdx)); + this->simulator().setTimeStepSize(eclWriter_->restartTimeStepSize()); + + eclWriter_->beginRestart(); size_t numElems = this->model().numGridDof(); initialFluidStates_.resize(numElems); @@ -1855,6 +1863,9 @@ private: polymerMoleWeight_.resize(numElems, 0.0); } + // this is a hack to preserve the initial fluid states + auto tmpInitialFs = initialFluidStates_; + for (size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) { auto& elemFluidState = initialFluidStates_[elemIdx]; elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx)); @@ -1887,6 +1898,31 @@ private: if (tracerModel().numTracers() > 0) std::cout << "Warning: Restart is not implemented for the tracer model, it will initialize with initial tracer concentration" << std::endl; + // assign the restart solution to the current solution. note that we still need + // to compute real initial solution after this because the initial fluid states + // need to be correct for stuff like boundary conditions. + auto& sol = this->model().solution(/*timeIdx=*/0); + const auto& gridView = this->gridView(); + ElementContext elemCtx(this->simulator()); + auto elemIt = gridView.template begin(); + const auto& elemEndIt = gridView.template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const auto& elem = *elemIt; + if (elem.partitionType() != Dune::InteriorEntity) + continue; + + elemCtx.updatePrimaryStencil(elem); + int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + initial(sol[elemIdx], elemCtx, /*spaceIdx=*/0, /*timeIdx=*/0); + } + + // this is a hack to preserve the initial fluid states + initialFluidStates_ = tmpInitialFs; + // make sure that the stuff which needs to be done at the beginning of an episode + // is run. + this->beginEpisode(/*isOnRestart=*/true); + + eclWriter_->endRestart(); } void processRestartSaturations_(InitialFluidState& elemFluidState) diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index 7e89d4420..a43ab0d7c 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -238,7 +238,7 @@ public: } } - void restartBegin() + void beginRestart() { bool enableHysteresis = simulator_.problem().materialLawManager()->enableHysteresis(); bool enableSwatinit = simulator_.vanguard().eclState().get3DProperties().hasDeckDoubleGridProperty("SWATINIT"); @@ -277,11 +277,18 @@ public: const auto& thpresValues = restartValues.getExtra("THRESHPR"); thpres.setFromRestart(thpresValues); } + restartTimeStepSize_ = restartValues.getExtra("OPMEXTRA")[0]; } + void endRestart() + {} + const EclOutputBlackOilModule& eclOutputModule() const { return eclOutputModule_; } + Scalar restartTimeStepSize() const + { return restartTimeStepSize_; } + private: static bool enableEclOutput_() @@ -481,6 +488,7 @@ private: std::unique_ptr eclIO_; Grid globalGrid_; std::unique_ptr taskletRunner_; + Scalar restartTimeStepSize_; }; From 65d44a055f69df8b9f13f521062e6bc7d1da111c Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 13 Feb 2019 12:54:52 +0100 Subject: [PATCH 2/3] EclProblem: fix a few minor style issues --- ebos/eclproblem.hh | 95 ++++++++++++++++++++-------------------------- 1 file changed, 41 insertions(+), 54 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index fceaea6f6..f162f9a55 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -357,7 +357,6 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) }; enum { enablePolymerMolarWeight = GET_PROP_VALUE(TypeTag, EnablePolymerMW) }; - enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; enum { enableThermalFluxBoundaries = GET_PROP_VALUE(TypeTag, EnableThermalFluxBoundaries) }; @@ -737,12 +736,12 @@ public: const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx); if (drsdtActive_()) // DRSDT is enabled - for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx ) + for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx) maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*this->simulator().timeStepSize(); if (drvdtActive_()) // DRVDT is enabled - for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx ) + for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*this->simulator().timeStepSize(); wellModel_.beginTimeStep(); @@ -958,7 +957,6 @@ public: Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const { return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); } - const EclThresholdPressure& thresholdPressure() const { return thresholdPressures_; } @@ -1064,7 +1062,6 @@ public: return thermalLawManager_->solidEnergyLawParams(globalSpaceIdx); } - /*! * \copydoc FvBaseMultiPhaseProblem::thermalConductionParams */ @@ -1113,7 +1110,6 @@ public: return polymerConcentration_[elemIdx]; } - /*! * \brief Returns the polymer molecule weight for a given cell index */ @@ -1191,7 +1187,7 @@ public: { return plmixnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); } /*! - * \brief Returns the index the relevant PLMIXNUM ( for polymer module) region given a cell index + * \brief Returns the index the relevant PLMIXNUM (for polymer module) region given a cell index */ unsigned plmixnumRegionIndex(unsigned elemIdx) const { @@ -1219,8 +1215,6 @@ public: return maxPolymerAdsorption_[elemIdx]; } - - /*! * \copydoc FvBaseProblem::name */ @@ -1266,17 +1260,15 @@ public: } if (hasFreeBoundaryConditions()) { - unsigned indexInInside = context.intersection(spaceIdx).indexInInside(); unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx); unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx); switch (indexInInside) { case 0: - { if (freebcXMinus_[globalDofIdx]) values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]); + break; - } case 1: if (freebcX_[globalDofIdx]) values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]); @@ -1400,7 +1392,7 @@ public: { int pvtRegionIdx = pvtRegionIndex(globalDofIdx); if (!drsdtActive_() || maxDRs_[pvtRegionIdx] < 0.0) - return std::numeric_limits::max()/2; + return std::numeric_limits::max()/2.0; // this is a bit hacky because it assumes that a time discretization with only // two time indices is used. @@ -1418,7 +1410,7 @@ public: { int pvtRegionIdx = pvtRegionIndex(globalDofIdx); if (!drvdtActive_() || maxDRv_[pvtRegionIdx] < 0.0) - return std::numeric_limits::max()/2; + return std::numeric_limits::max()/2.0; // this is a bit hacky because it assumes that a time discretization with only // two time indices is used. @@ -1470,7 +1462,7 @@ public: { return wellModel_; } // temporary solution to facilitate output of initial state from flow - const InitialFluidState& initialFluidState(unsigned globalDofIdx ) const + const InitialFluidState& initialFluidState(unsigned globalDofIdx) const { return initialFluidStates_[globalDofIdx]; } const Opm::EclipseIO& eclIO() const @@ -1478,7 +1470,7 @@ public: bool vapparsActive() const { - int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0 ); + int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0); const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx); return (oilVaporizationControl.getType() == Opm::OilVaporizationEnum::VAPPARS); } @@ -1489,31 +1481,30 @@ public: private: bool drsdtActive_() const { - int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0 ); + int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0); const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx); return (oilVaporizationControl.drsdtActive()); - } + bool drvdtActive_() const { - int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0 ); + int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0); const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx); return (oilVaporizationControl.drvdtActive()); } - Scalar cellCenterDepth( const Element& element ) const + + Scalar cellCenterDepth(const Element& element) const { - typedef typename Element :: Geometry Geometry; - static constexpr int zCoord = Element :: dimension - 1; + typedef typename Element::Geometry Geometry; + static constexpr int zCoord = Element::dimension - 1; Scalar zz = 0.0; const Geometry geometry = element.geometry(); - const int corners = geometry.corners(); - for (int i=0; isimulator().episodeIndex(), 0 ); + int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0); const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx); if (oilVaporizationControl.drsdtActive()) { @@ -1559,7 +1550,7 @@ private: const auto& iq = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& fs = iq.fluidState(); - typedef typename std::decay::type FluidState; + typedef typename std::decay::type FluidState; int pvtRegionIdx = pvtRegionIndex(compressedDofIdx); if (oilVaporizationControl.getOption(pvtRegionIdx) || fs.saturation(gasPhaseIdx) > freeGasMinSaturation_) @@ -1589,7 +1580,7 @@ private: const auto& iq = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& fs = iq.fluidState(); - typedef typename std::decay::type FluidState; + typedef typename std::decay::type FluidState; lastRv_[compressedDofIdx] = Opm::BlackOil::template getRv_simulator().vanguard().schedule(); const auto& eclState = this->simulator().vanguard().eclState(); const auto& timeMap = this->simulator().vanguard().schedule().getTimeMap(); const auto& initconfig = eclState.getInitConfig(); @@ -1850,10 +1842,10 @@ private: size_t numElems = this->model().numGridDof(); initialFluidStates_.resize(numElems); if (enableSolvent) - solventSaturation_.resize(numElems,0.0); + solventSaturation_.resize(numElems, 0.0); if (enablePolymer) - polymerConcentration_.resize(numElems,0.0); + polymerConcentration_.resize(numElems, 0.0); if (enablePolymerMolarWeight) { const std::string msg {"Support of the RESTART for polymer molecular weight " @@ -1887,12 +1879,12 @@ private: const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx); if (drsdtActive_()) // DRSDT is enabled - for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx ) + for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx) maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*this->simulator().timeStepSize(); if (drvdtActive_()) // DRVDT is enabled - for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx ) + for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*this->simulator().timeStepSize(); if (tracerModel().numTracers() > 0) @@ -1930,17 +1922,17 @@ private: // each phase needs to be above certain value to be claimed to be existing // this is used to recover some RESTART running with the defaulted single-precision format const Scalar smallSaturationTolerance = 1.e-6; - Scalar sumSaturation = 0.; + Scalar sumSaturation = 0.0; for (size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (FluidSystem::phaseIsActive(phaseIdx)) { if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance) - elemFluidState.setSaturation(phaseIdx, 0.); + elemFluidState.setSaturation(phaseIdx, 0.0); sumSaturation += elemFluidState.saturation(phaseIdx); } } - assert(sumSaturation > 0.); + assert(sumSaturation > 0.0); for (size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { if (FluidSystem::phaseIsActive(phaseIdx)) { @@ -2055,7 +2047,7 @@ private: // this assumes that capillary pressures only depend on the phase saturations // and possibly on temperature. (this is always the case for ECL problems.) - Dune::FieldVector< Scalar, numPhases > pc( 0 ); + Dune::FieldVector pc(0.0); const auto& matParams = materialLawParams(dofIdx); MaterialLaw::capillaryPressures(pc, matParams, dofFluidState); Opm::Valgrind::CheckDefined(oilPressure); @@ -2102,7 +2094,7 @@ private: if (enableSolvent) { const std::vector& solventSaturationData = eclState.get3DProperties().getDoubleGridProperty("SSOL").getData(); - solventSaturation_.resize(numDof,0.0); + solventSaturation_.resize(numDof, 0.0); for (size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { size_t cartesianDofIdx = vanguard.cartesianIndex(dofIdx); assert(0 <= cartesianDofIdx); @@ -2113,7 +2105,7 @@ private: if (enablePolymer) { const std::vector& polyConcentrationData = eclState.get3DProperties().getDoubleGridProperty("SPOLY").getData(); - polymerConcentration_.resize(numDof,0.0); + polymerConcentration_.resize(numDof, 0.0); for (size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { size_t cartesianDofIdx = vanguard.cartesianIndex(dofIdx); assert(0 <= cartesianDofIdx); @@ -2124,7 +2116,7 @@ private: if (enablePolymerMolarWeight) { const std::vector& polyMoleWeightData = eclState.get3DProperties().getDoubleGridProperty("SPOLYMW").getData(); - polymerMoleWeight_.resize(numDof,0.0); + polymerMoleWeight_.resize(numDof, 0.0); for (size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { const size_t cartesianDofIdx = vanguard.cartesianIndex(dofIdx); assert(0 <= cartesianDofIdx); @@ -2132,11 +2124,8 @@ private: polymerMoleWeight_[dofIdx] = polyMoleWeightData[cartesianDofIdx]; } } - } - - // update the hysteresis parameters of the material laws for the whole grid bool updateHysteresis_() { @@ -2291,15 +2280,15 @@ private: void readBoundaryConditions_() { hasFreeBoundaryConditions_ = false; - readBoundaryConditionKeyword_("FREEBCX", freebcX_ ); - readBoundaryConditionKeyword_("FREEBCX-", freebcXMinus_ ); - readBoundaryConditionKeyword_("FREEBCY", freebcY_ ); - readBoundaryConditionKeyword_("FREEBCY-", freebcYMinus_ ); - readBoundaryConditionKeyword_("FREEBCZ", freebcZ_ ); - readBoundaryConditionKeyword_("FREEBCZ-", freebcZMinus_ ); + readBoundaryConditionKeyword_("FREEBCX", freebcX_); + readBoundaryConditionKeyword_("FREEBCX-", freebcXMinus_); + readBoundaryConditionKeyword_("FREEBCY", freebcY_); + readBoundaryConditionKeyword_("FREEBCY-", freebcYMinus_); + readBoundaryConditionKeyword_("FREEBCZ", freebcZ_); + readBoundaryConditionKeyword_("FREEBCZ-", freebcZMinus_); } - void readBoundaryConditionKeyword_(const std::string& name, std::vector& compressedData ) + void readBoundaryConditionKeyword_(const std::string& name, std::vector& compressedData) { const auto& eclProps = this->simulator().vanguard().eclState().get3DProperties(); const auto& vanguard = this->simulator().vanguard(); @@ -2368,8 +2357,6 @@ private: std::vector freebcYMinus_; std::vector freebcZ_; std::vector freebcZMinus_; - - }; template From 5b032a6a280568215aefcf1c2f426ab4cd74a3f6 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 19 Feb 2019 09:07:12 +0100 Subject: [PATCH 3/3] fix the issues found by [at]tosa82 in his review in particular the missing synchronization after restarts was very nasty to find. thanks a ton for pointing this out! also, IIRC changing DR[SV]DT in the schedule section has been working properly for a while, so the comment which stated the opposite is removed as well. --- ebos/eclproblem.hh | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index f162f9a55..74fe4acb4 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -552,11 +552,6 @@ public: if (!deck.hasKeyword("NOGRAV") && EWOMS_GET_PARAM(TypeTag, bool, EnableGravity)) this->gravity_[dim - 1] = 9.80665; - // this is actually not fully correct: the latest occurence of VAPPARS and DRSDT - // or DRVDT up to the current time step in the schedule section counts, presence - // of VAPPARS alone is not sufficient to disable DR[SV]DT. TODO: implment support - // for this in opm-parser's Schedule object" - // deal with DRSDT const auto& eclState = simulator.vanguard().eclState(); unsigned ntpvt = eclState.runspec().tabdims().getNumPVTTables(); @@ -575,7 +570,6 @@ public: readThermalParameters_(); transmissibilities_.finishInit(); - const auto& initconfig = eclState.getInitConfig(); const auto& timeMap = simulator.vanguard().schedule().getTimeMap(); readInitialCondition_(); @@ -1827,7 +1821,7 @@ private: // Set the start time of the simulation const auto& schedule = this->simulator().vanguard().schedule(); const auto& eclState = this->simulator().vanguard().eclState(); - const auto& timeMap = this->simulator().vanguard().schedule().getTimeMap(); + const auto& timeMap = schedule.getTimeMap(); const auto& initconfig = eclState.getInitConfig(); int episodeIdx = initconfig.getRestartStep() - 1; @@ -1908,6 +1902,13 @@ private: initial(sol[elemIdx], elemCtx, /*spaceIdx=*/0, /*timeIdx=*/0); } + // make sure that the ghost and overlap entities exhibit the correct + // solution. alternatively, this could be done in the loop above by also + // considering non-interior elements. Since the initial() method might not work + // 100% correctly for such elements, let's play safe and explicitly synchronize + // using message passing. + this->model().syncOverlap(); + // this is a hack to preserve the initial fluid states initialFluidStates_ = tmpInitialFs; // make sure that the stuff which needs to be done at the beginning of an episode