From cdb173f6d1f56f0873ed20417a3e27e89c9c9631 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 15 Nov 2017 11:45:02 +0100 Subject: [PATCH 01/38] eclproblem option to avoid restart spesific deserialize added comment on probable wrong behavoir of deserialize which also is restart spesific Added grid name in restart --- ebos/eclproblem.hh | 9 ++++++--- ewoms/common/simulator.hh | 12 +++++++++++- ewoms/disc/ecfv/ecfvdiscretization.hh | 2 +- ewoms/io/restart.hh | 7 ++++--- 4 files changed, 22 insertions(+), 8 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index ce409297d8..d63f73342d 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -622,11 +622,14 @@ public: * \param res The deserializer object */ template - void deserialize(Restarter& res) + void deserialize(Restarter& res, bool isOnRestart) { // reload the current episode/report step from the deck - beginEpisode(/*isOnRestart=*/true); - + //beginEpisode(/*isOnRestart=*/true); + // should probably be removed + if(isOnRestart){ + beginEpisode(/*isOnRestart=*/true); + } // deserialize the wells wellModel_.deserialize(res); } diff --git a/ewoms/common/simulator.hh b/ewoms/common/simulator.hh index 57bc9e9238..7691cae68c 100644 --- a/ewoms/common/simulator.hh +++ b/ewoms/common/simulator.hh @@ -550,7 +550,7 @@ public: if (verbose_) std::cout << "Deserialize from file '" << res.fileName() << "'\n" << std::flush; this->deserialize(res); - problem_->deserialize(res); + problem_->deserialize(res, true);// to things for restart model_->deserialize(res); res.deserializeEnd(); if (verbose_) @@ -808,6 +808,16 @@ public: res.serializeEnd(); } + void deserializeAll(Scalar t) + { + typedef Ewoms::Restart Restarter; + Restarter res; + res.deserializeBegin(*this, t); + this->deserialize(res); + problem_->deserialize(res, false); + model_->deserialize(res); + res.deserializeEnd(); + } /*! * \brief Write the time manager's state to a restart file. * diff --git a/ewoms/disc/ecfv/ecfvdiscretization.hh b/ewoms/disc/ecfv/ecfvdiscretization.hh index 9ae5805e5b..913f88a814 100644 --- a/ewoms/disc/ecfv/ecfvdiscretization.hh +++ b/ewoms/disc/ecfv/ecfvdiscretization.hh @@ -202,7 +202,7 @@ public: void deserialize(Restarter& res) { res.template deserializeEntities(asImp_(), this->gridView_); - this->solution(/*timeIdx=*/1) = this->solution(/*timeIdx=*/0); + this->solution(/*timeIdx=*/1) = this->solution(/*timeIdx=*/0);// may wrong and only valid for restart } private: diff --git a/ewoms/io/restart.hh b/ewoms/io/restart.hh index fc0633c829..171f769f91 100644 --- a/ewoms/io/restart.hh +++ b/ewoms/io/restart.hh @@ -46,7 +46,7 @@ class Restart template static const std::string magicRestartCookie_(const GridView& gridView) { - static const std::string gridName = "blubb"; // gridView.grid().name(); + static const std::string gridName = gridView.grid().name(); static const int dim = GridView::dimension; int numVertices = gridView.size(dim); @@ -105,7 +105,6 @@ public: simulator.problem().outputDir(), simulator.problem().name(), simulator.time()); - // open output file and write magic cookie outStream_.open(fileName_.c_str()); outStream_.precision(20); @@ -164,10 +163,12 @@ public: void serializeEnd() { outStream_.close(); } - /*! + + /*! * \brief Start reading a restart file at a certain simulated * time. */ + template void deserializeBegin(Simulator& simulator, Scalar t) { From 231fe5e223bc19a9a0a1f80d02916a08db7fca36 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 16 Nov 2017 21:24:20 +0100 Subject: [PATCH 02/38] Tried to fix serialization of simulator by adding time step to it commented probably unwanted code --- ewoms/common/simulator.hh | 6 ++++++ ewoms/models/blackoil/blackoilmodel.hh | 1 + 2 files changed, 7 insertions(+) diff --git a/ewoms/common/simulator.hh b/ewoms/common/simulator.hh index 7691cae68c..1c19468ad6 100644 --- a/ewoms/common/simulator.hh +++ b/ewoms/common/simulator.hh @@ -350,11 +350,15 @@ public: */ Scalar timeStepSize() const { + return timeStepSize_; + + /* This should be handeled by ghe time stepper Scalar maximumTimeStepSize = std::min(episodeMaxTimeStepSize(), std::max( Scalar(0), endTime() - this->time())); return std::min(timeStepSize_, maximumTimeStepSize); + */ } /*! @@ -835,6 +839,7 @@ public: << episodeLength_ << " " << startTime_ << " " << time_ << " " + << timeStepSize_ << " " << timeStepIdx_ << " "; restarter.serializeSectionEnd(); } @@ -856,6 +861,7 @@ public: >> episodeLength_ >> startTime_ >> time_ + >> timeStepSize_ >> timeStepIdx_; restarter.deserializeSectionEnd(); } diff --git a/ewoms/models/blackoil/blackoilmodel.hh b/ewoms/models/blackoil/blackoilmodel.hh index fc8c6c22d7..7bc36e1821 100644 --- a/ewoms/models/blackoil/blackoilmodel.hh +++ b/ewoms/models/blackoil/blackoilmodel.hh @@ -503,6 +503,7 @@ public: } } + // this is restart spesific behavoir and to avoid this one have to hack outside for adjiont this->solution(/*timeIdx=*/1) = this->solution(/*timeIdx=*/0); } From 36bfda86110b71691f38defbb084526520fe397f Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 16 Nov 2017 22:02:06 +0100 Subject: [PATCH 03/38] disable storageCache --- ebos/eclproblem.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index d63f73342d..01c601c89b 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -262,7 +262,7 @@ SET_STRING_PROP(EclBaseProblem, OutputDir, "."); SET_BOOL_PROP(EclBaseProblem, EnableIntensiveQuantityCache, true); // the cache for the storage term can also be used and also yields a decent speedup -SET_BOOL_PROP(EclBaseProblem, EnableStorageCache, true); +SET_BOOL_PROP(EclBaseProblem, EnableStorageCache, false); // Use the "velocity module" which uses the Eclipse "NEWTRAN" transmissibilities SET_TYPE_PROP(EclBaseProblem, FluxModule, Ewoms::EclTransFluxModule); From 0d2f024b76475fabfbd51223356b45ebbef186d9 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 1 Feb 2018 16:36:55 +0100 Subject: [PATCH 04/38] use storage cach on --- ebos/eclproblem.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 01c601c89b..d63f73342d 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -262,7 +262,7 @@ SET_STRING_PROP(EclBaseProblem, OutputDir, "."); SET_BOOL_PROP(EclBaseProblem, EnableIntensiveQuantityCache, true); // the cache for the storage term can also be used and also yields a decent speedup -SET_BOOL_PROP(EclBaseProblem, EnableStorageCache, false); +SET_BOOL_PROP(EclBaseProblem, EnableStorageCache, true); // Use the "velocity module" which uses the Eclipse "NEWTRAN" transmissibilities SET_TYPE_PROP(EclBaseProblem, FluxModule, Ewoms::EclTransFluxModule); From 4fa97f7a3eb6ca78c27ff04995e4172bce610804 Mon Sep 17 00:00:00 2001 From: hnil Date: Fri, 2 Feb 2018 01:18:33 +0100 Subject: [PATCH 05/38] tried to add timefocus --- ewoms/disc/common/fvbaseelementcontext.hh | 10 ++++++++-- ewoms/disc/common/fvbaseprimaryvariables.hh | 4 ++-- .../blackoil/blackoilintensivequantities.hh | 20 +++++++++---------- 3 files changed, 20 insertions(+), 14 deletions(-) diff --git a/ewoms/disc/common/fvbaseelementcontext.hh b/ewoms/disc/common/fvbaseelementcontext.hh index eea92af2e0..ddc7a218e7 100644 --- a/ewoms/disc/common/fvbaseelementcontext.hh +++ b/ewoms/disc/common/fvbaseelementcontext.hh @@ -100,6 +100,7 @@ public: enableStorageCache_ = EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache); stashedDofIdx_ = -1; focusDofIdx_ = -1; + focusTimeIdx_ = 0; } static void *operator new(size_t size) { @@ -261,7 +262,8 @@ public: */ void setFocusDofIndex(unsigned dofIdx) { focusDofIdx_ = dofIdx; } - + void setFocusTimeIndex(unsigned timeIdx) + { focusTimeIdx_ = timeIdx; } /*! * \brief Returns the degree of freedom on which the simulator is currently "focused" on * @@ -270,6 +272,9 @@ public: unsigned focusDofIndex() const { return focusDofIdx_; } + unsigned focusTimeIndex() const + { return focusTimeIdx_; } + /*! * \brief Return a reference to the simulator. */ @@ -583,7 +588,7 @@ protected: #endif dofVars_[dofIdx].priVars[timeIdx] = priVars; - dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/asImp_(), dofIdx, timeIdx); + dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/asImp_(), dofIdx, timeIdx, focusTimeIdx_); } IntensiveQuantities intensiveQuantitiesStashed_; @@ -601,6 +606,7 @@ protected: int stashedDofIdx_; int focusDofIdx_; + int focusTimeIdx_; bool enableStorageCache_; }; diff --git a/ewoms/disc/common/fvbaseprimaryvariables.hh b/ewoms/disc/common/fvbaseprimaryvariables.hh index 8d87d4a0e4..01962ea4a2 100644 --- a/ewoms/disc/common/fvbaseprimaryvariables.hh +++ b/ewoms/disc/common/fvbaseprimaryvariables.hh @@ -87,9 +87,9 @@ public: * it represents the a constant f = x_i. (the difference is that in the first case, * the derivative w.r.t. x_i is 1, while it is 0 in the second case. */ - Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx) const + Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx,unsigned focustimeidx) const { - if (timeIdx == 0) + if (timeIdx == focustimeidx) return Toolbox::createVariable((*this)[varIdx], varIdx); else return Toolbox::createConstant((*this)[varIdx]); diff --git a/ewoms/models/blackoil/blackoilintensivequantities.hh b/ewoms/models/blackoil/blackoilintensivequantities.hh index 2f7edc96b4..7678856dc3 100644 --- a/ewoms/models/blackoil/blackoilintensivequantities.hh +++ b/ewoms/models/blackoil/blackoilintensivequantities.hh @@ -107,9 +107,9 @@ public: /*! * \copydoc IntensiveQuantities::update */ - void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx) + void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focustimeidx) { - ParentType::update(elemCtx, dofIdx, timeIdx); + ParentType::update(elemCtx, dofIdx, timeIdx);// this only gives extrusion factor const auto& problem = elemCtx.problem(); const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx); @@ -123,21 +123,21 @@ public: // extract the water and the gas saturations for convenience Evaluation Sw = 0.0; if (waterEnabled) - Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx); + Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx, focustimeidx); Evaluation Sg = 0.0; if (compositionSwitchEnabled) { if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) // -> threephase case - Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx); + Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focustimeidx); else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) { // -> gas-water case Sg = 1.0 - Sw; // deal with solvent if (enableSolvent) - Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx); + Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx, focustimeidx); } else { @@ -154,7 +154,7 @@ public: // deal with solvent if (enableSolvent) - So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx); + So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx, focustimeidx); fluidState_.setSaturation(waterPhaseIdx, Sw); fluidState_.setSaturation(gasPhaseIdx, Sg); @@ -169,13 +169,13 @@ public: //oil is the reference phase for pressure if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) { - const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx); + const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focustimeidx); for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx])); } else { - const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx); + const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focustimeidx); for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx])); } @@ -224,7 +224,7 @@ public: Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(globalSpaceIdx); // oil phase, we can directly set the composition of the oil phase - const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx); + const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focustimeidx); fluidState_.setRs(Opm::min(RsMax, Rs)); if (FluidSystem::enableVaporizedOil()) { @@ -245,7 +245,7 @@ public: else { assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv); - const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx); + const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focustimeidx); fluidState_.setRv(Rv); if (FluidSystem::enableDissolvedGas()) { From 19bd8e8094ea9a4d72a23dca4c032284c924fa1a Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 6 Feb 2018 11:38:58 +0100 Subject: [PATCH 06/38] changes to make focus time work in particular for storage term, need retinging of consepts of derivatives... --- ewoms/disc/common/fvbaselinearizer.hh | 19 ++++++++++--------- ewoms/disc/common/fvbaselocalresidual.hh | 8 +++++++- 2 files changed, 17 insertions(+), 10 deletions(-) diff --git a/ewoms/disc/common/fvbaselinearizer.hh b/ewoms/disc/common/fvbaselinearizer.hh index 50b69cd55c..5aff556588 100644 --- a/ewoms/disc/common/fvbaselinearizer.hh +++ b/ewoms/disc/common/fvbaselinearizer.hh @@ -163,10 +163,10 @@ public: * * This means the spatial domain plus all auxiliary equations. */ - void linearize() + void linearize(unsigned focustimeIndex) { - linearizeDomain(); - linearizeAuxiliaryEquations(); + linearizeDomain(focustimeIndex); + linearizeAuxiliaryEquations(focustimeIndex); } /*! @@ -179,7 +179,7 @@ public: * The current state of affairs (esp. the previous and the current solutions) is * represented by the model object. */ - void linearizeDomain() + void linearizeDomain(unsigned focustimeIndex) { // we defer the initialization of the Jacobian matrix until here because the // auxiliary modules usually assume the problem, model and grid to be fully @@ -189,7 +189,7 @@ public: int succeeded; try { - linearize_(); + linearize_(focustimeIndex); succeeded = 1; } catch (const std::exception& e) @@ -225,7 +225,7 @@ public: * \brief Linearize the part of the non-linear system of equations that is associated * with the spatial domain. */ - void linearizeAuxiliaryEquations() + void linearizeAuxiliaryEquations(OPM_UNUSED unsigned focustimeIndex) { auto& model = model_(); const auto& comm = simulator_().gridView().comm(); @@ -434,7 +434,7 @@ private: } // linearize the whole system - void linearize_() + void linearize_(unsigned focustimeindex) { resetSystem_(); @@ -474,7 +474,7 @@ private: if (!linearizeNonLocalElements && elem.partitionType() != Dune::InteriorEntity) continue; - linearizeElement_(elem); + linearizeElement_(elem, focustimeindex); } } @@ -482,7 +482,7 @@ private: } // linearize an element in the interior of the process' grid partition - void linearizeElement_(const Element& elem) + void linearizeElement_(const Element& elem, unsigned focustimeindex) { unsigned threadId = ThreadManager::threadId(); @@ -490,6 +490,7 @@ private: auto& localLinearizer = model_().localLinearizer(threadId); // the actual work of linearization is done by the local linearizer class + elementCtx->setFocusTimeIndex(focustimeindex); localLinearizer.linearize(*elementCtx, elem); // update the right hand side and the Jacobian matrix diff --git a/ewoms/disc/common/fvbaselocalresidual.hh b/ewoms/disc/common/fvbaselocalresidual.hh index b6a36bf79c..ac763816d3 100644 --- a/ewoms/disc/common/fvbaselocalresidual.hh +++ b/ewoms/disc/common/fvbaselocalresidual.hh @@ -491,7 +491,8 @@ protected: const ElementContext& elemCtx) const { EvalVector tmp; - EqVector tmp2; + //EqVector tmp2; + EvalVector tmp2; RateVector sourceRate; tmp = 0.0; @@ -566,6 +567,11 @@ protected: // Use the implicit Euler time discretization for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { + if(elemCtx.focusTimeIndex()==0){ + tmp2[eqIdx] .clearDerivatives();// remove derivatives for the time index + }else{ + tmp[eqIdx] .clearDerivatives(); + } tmp[eqIdx] -= tmp2[eqIdx]; tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize(); From d2adee248f886bfc38a88b38b2ec580f949bf028 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 14 Feb 2018 14:42:08 +0100 Subject: [PATCH 07/38] added numAdjoint in ecl def, and fixed storage term to work also with adjoint --- ebos/eclproblem.hh | 4 ++ ewoms/disc/common/fvbaselocalresidual.hh | 65 ++++++++++++++++++++---- 2 files changed, 58 insertions(+), 11 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index d63f73342d..d2e8e669f5 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -237,6 +237,10 @@ SET_SCALAR_PROP(EclBaseProblem, NewtonMaxError, 0.1); // step succeeds at more than 14 Newton iteration is rather small SET_INT_PROP(EclBaseProblem, NewtonMaxIterations, 14); +// set no adjoint variables ad default +NEW_PROP_TAG(numAdjoint); +SET_INT_PROP(EclBaseProblem, numAdjoint, 0); + // also, reduce the target for the "optimum" number of Newton iterations to 6. Note that // this is only relevant if the time step is reduced from the report step size for some // reason. (because ebos first tries to do a report step using a single time step.) diff --git a/ewoms/disc/common/fvbaselocalresidual.hh b/ewoms/disc/common/fvbaselocalresidual.hh index ac763816d3..7471ac1cb0 100644 --- a/ewoms/disc/common/fvbaselocalresidual.hh +++ b/ewoms/disc/common/fvbaselocalresidual.hh @@ -491,8 +491,8 @@ protected: const ElementContext& elemCtx) const { EvalVector tmp; - //EqVector tmp2; - EvalVector tmp2; + EqVector tmp2; + EvalVector tmp2_der;// this is added to support calculation of storage term derivative focusTime!=0 for need for adjoint: always used if cachedStororage term si false RateVector sourceRate; tmp = 0.0; @@ -565,19 +565,62 @@ protected: Opm::Valgrind::CheckDefined(tmp2); } - // Use the implicit Euler time discretization - for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { - if(elemCtx.focusTimeIndex()==0){ - tmp2[eqIdx] .clearDerivatives();// remove derivatives for the time index - }else{ - tmp[eqIdx] .clearDerivatives(); + if (elemCtx.enableStorageCache()) { + const auto& model = elemCtx.model(); + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + if (model.newtonMethod().numIterations() == 0 && + !elemCtx.haveStashedIntensiveQuantities()) + { + // if the storage term is cached and we're in the first iteration of + // the time step, update the cache of the storage term (this assumes + // that the initial guess for the solution at the end of the time + // step is the same as the solution at the beginning of the time + // step. This is usually true, but some fancy preprocessing scheme + // might invalidate that assumption.) + for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) + tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]); + Opm::Valgrind::CheckDefined(tmp2); + + model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2); + } + else { + // if the storage term is cached and we're not looking at the first + // iteration of the time step, we take the cached data. + tmp2 = model.cachedStorage(globalDofIdx, /*timeIdx=*/1); + Opm::Valgrind::CheckDefined(tmp2); } - tmp[eqIdx] -= tmp2[eqIdx]; - tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize(); + } else { + // if the mass storage at the beginning of the time step is not cached, + // we re-calculate it from scratch. + // if storage cache not enables we always use derivatives + tmp2_der = 0.0; + asImp_().computeStorage(tmp2_der, elemCtx, dofIdx, /*timeIdx=*/1); + Opm::Valgrind::CheckDefined(tmp2_der); + } - residual[dofIdx][eqIdx] += tmp[eqIdx]; + if( !elemCtx.enableStorageCache() ){ + // Use the implicit Euler time discretization + for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { + if(elemCtx.focusTimeIndex()==0){ + tmp2_der[eqIdx] .clearDerivatives();// remove derivatives for the time index + }else{ + tmp[eqIdx] .clearDerivatives(); + } + tmp[eqIdx] -= tmp2_der[eqIdx]; + tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize(); + residual[dofIdx][eqIdx] += tmp[eqIdx]; + } + }else{ + for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { + // cached storage should only be used in forward mode + assert(elemCtx.focusTimeIndex()==0); + tmp[eqIdx] -= tmp2[eqIdx]; + tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize(); + residual[dofIdx][eqIdx] += tmp[eqIdx]; + } } + Opm::Valgrind::CheckDefined(residual[dofIdx]); // deal with the source term From fb5c54f10ac37afaedf81baf0cba68a9cf990124 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 14 Feb 2018 16:59:33 +0100 Subject: [PATCH 08/38] fixed time focus for solvent and polymer --- .../blackoil/blackoilintensivequantities.hh | 11 +++++---- .../models/blackoil/blackoilpolymermodules.hh | 9 ++++--- .../models/blackoil/blackoilsolventmodules.hh | 24 ++++++++++++------- 3 files changed, 27 insertions(+), 17 deletions(-) diff --git a/ewoms/models/blackoil/blackoilintensivequantities.hh b/ewoms/models/blackoil/blackoilintensivequantities.hh index 7678856dc3..3e1ee67e0d 100644 --- a/ewoms/models/blackoil/blackoilintensivequantities.hh +++ b/ewoms/models/blackoil/blackoilintensivequantities.hh @@ -160,7 +160,7 @@ public: fluidState_.setSaturation(gasPhaseIdx, Sg); fluidState_.setSaturation(oilPhaseIdx, So); - asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx); + asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx, focustimeidx); // now we compute all phase pressures Evaluation pC[numPhases]; @@ -186,7 +186,7 @@ public: Opm::Valgrind::CheckDefined(mobility_); // update the Saturation functions for the blackoil solvent module. - asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx); + asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx, focustimeidx); Scalar SoMax = elemCtx.problem().maxOilSaturation(globalSpaceIdx); @@ -327,12 +327,13 @@ public: porosity_ *= 1.0 + x + 0.5*x*x; } - asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx); - asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx); - asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache); + asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx, focustimeidx);// strictly do not need focus time since no makeEvaluation is done + asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx, focustimeidx); + asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache, focustimeidx); // update the quantities which are required by the chosen // velocity model + // TODO: use the focus time here FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx); #ifndef NDEBUG diff --git a/ewoms/models/blackoil/blackoilpolymermodules.hh b/ewoms/models/blackoil/blackoilpolymermodules.hh index 8614a843bd..84114b9217 100644 --- a/ewoms/models/blackoil/blackoilpolymermodules.hh +++ b/ewoms/models/blackoil/blackoilpolymermodules.hh @@ -861,10 +861,11 @@ public: */ void polymerPropertiesUpdate_(const ElementContext& elemCtx, unsigned dofIdx, - unsigned timeIdx) + unsigned timeIdx, + unsigned focustime) { const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx); - polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx); + polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx, focustime); const Scalar cmax = PolymerModule::plymaxMaxConcentration(elemCtx, dofIdx, timeIdx); // permeability reduction due to polymer @@ -952,7 +953,9 @@ class BlackOilPolymerIntensiveQuantities public: void polymerPropertiesUpdate_(const ElementContext& elemCtx OPM_UNUSED, unsigned scvIdx OPM_UNUSED, - unsigned timeIdx OPM_UNUSED) + unsigned timeIdx OPM_UNUSED, + unsigned focustimeIdx OPM_UNUSED) + { } const Evaluation& polymerConcentration() const diff --git a/ewoms/models/blackoil/blackoilsolventmodules.hh b/ewoms/models/blackoil/blackoilsolventmodules.hh index 2392ba4e2e..f7d7a72adb 100644 --- a/ewoms/models/blackoil/blackoilsolventmodules.hh +++ b/ewoms/models/blackoil/blackoilsolventmodules.hh @@ -899,11 +899,12 @@ public: */ void solventPreSatFuncUpdate_(const ElementContext& elemCtx, unsigned dofIdx, - unsigned timeIdx) + unsigned timeIdx, + unsigned focustimeIdx) { const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx); auto& fs = asImp_().fluidState_; - solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx); + solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx,focustimeIdx); hydrocarbonSaturation_ = fs.saturation(gasPhaseIdx); // apply a cut-off. Don't waste calculations if no solvent @@ -924,7 +925,8 @@ public: */ void solventPostSatFuncUpdate_(const ElementContext& elemCtx, unsigned dofIdx, - unsigned timeIdx) + unsigned timeIdx, + unsigned focustimeidx) { // revert the gas "saturation" of the fluid state back to the saturation of the // hydrocarbon gas. @@ -953,9 +955,9 @@ public: //oil is the reference phase for pressure if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) - pgMisc = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx); + pgMisc = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focustimeidx); else { - const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx); + const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focustimeidx); pgMisc = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]); } @@ -1039,7 +1041,8 @@ public: */ void solventPvtUpdate_(const ElementContext& elemCtx, unsigned scvIdx, - unsigned timeIdx) + unsigned timeIdx, + unsigned focustimeidx) { const auto& iq = asImp_(); const auto& fs = iq.fluidState(); @@ -1268,17 +1271,20 @@ class BlackOilSolventIntensiveQuantities public: void solventPreSatFuncUpdate_(const ElementContext& elemCtx OPM_UNUSED, unsigned scvIdx OPM_UNUSED, - unsigned timeIdx OPM_UNUSED) + unsigned timeIdx OPM_UNUSED, + unsigned focustimeidx OPM_UNUSED) { } void solventPostSatFuncUpdate_(const ElementContext& elemCtx OPM_UNUSED, unsigned scvIdx OPM_UNUSED, - unsigned timeIdx OPM_UNUSED) + unsigned timeIdx OPM_UNUSED, + unsigned focustimeidx OPM_UNUSED ) { } void solventPvtUpdate_(const ElementContext& elemCtx OPM_UNUSED, unsigned scvIdx OPM_UNUSED, - unsigned timeIdx OPM_UNUSED) + unsigned timeIdx OPM_UNUSED, + unsigned focustimeIdx OPM_UNUSED) { } const Evaluation& solventSaturation() const From a6cafd5eb9141ef6c91eb8699b02d5eecbb3fb91 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 14 Feb 2018 18:38:18 +0100 Subject: [PATCH 09/38] fix after timefocus added --- ewoms/nonlinear/newtonmethod.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ewoms/nonlinear/newtonmethod.hh b/ewoms/nonlinear/newtonmethod.hh index 98d1fd6721..5210592694 100644 --- a/ewoms/nonlinear/newtonmethod.hh +++ b/ewoms/nonlinear/newtonmethod.hh @@ -643,10 +643,10 @@ protected: * spatial domain. */ void linearizeDomain_() - { model().linearizer().linearizeDomain(); } + { model().linearizer().linearizeDomain(/*focus time*/ 0); } void linearizeAuxiliaryEquations_() - { model().linearizer().linearizeAuxiliaryEquations(); } + { model().linearizer().linearizeAuxiliaryEquations(/*focus time*/ 0); } void preSolve_(const SolutionVector& currentSolution OPM_UNUSED, const GlobalEqVector& currentResidual) From f3167a7672477e1a50be1455a73bbe26e7eb7b28 Mon Sep 17 00:00:00 2001 From: hnil Date: Sun, 18 Feb 2018 17:36:52 +0100 Subject: [PATCH 10/38] introduced directory for ebos serialization, NB probably break all ewoms not using eclproblem.hh --- ewoms/io/restart.hh | 35 +++++++++++++++++++++++------------ 1 file changed, 23 insertions(+), 12 deletions(-) diff --git a/ewoms/io/restart.hh b/ewoms/io/restart.hh index 171f769f91..82cee57281 100644 --- a/ewoms/io/restart.hh +++ b/ewoms/io/restart.hh @@ -31,6 +31,7 @@ #include #include #include +# include namespace Ewoms { @@ -70,21 +71,32 @@ class Restart * \brief Return the restart file name. */ template - static const std::string restartFileName_(const GridView& gridView, - const std::string& outputDir, - const std::string& simName, - Scalar t) + static const std::string restartFileName_(const Simulator& simulator, Scalar t) { - std::string dir = outputDir; + // Directory. + std::string dir = simulator.problem().outputDir(); if (dir == ".") dir = ""; else if (!dir.empty() && dir.back() != '/') dir += "/"; - + namespace fs = boost::filesystem; + fs::path output_dir(dir); + fs::path subdir("ebos_restart"); + output_dir = output_dir / subdir; + if(!(fs::exists(output_dir))){ + fs::create_directory(output_dir); + } + + // Filename. int rank = gridView.comm().rank(); + std::string simName = simulator.problem().name(); std::ostringstream oss; - oss << dir << simName << "_time=" << t << "_rank=" << rank << ".ers"; - return oss.str(); + oss << simName << "_time=" << t << "_rank=" << rank << ".ers"; + fs::path output_file(oss.str()); + + // Combine and return. + fs::path full_path = output_dir / output_file; + return full_path.string(); } public: @@ -101,10 +113,9 @@ public: void serializeBegin(Simulator& simulator) { const std::string magicCookie = magicRestartCookie_(simulator.gridView()); - fileName_ = restartFileName_(simulator.gridView(), - simulator.problem().outputDir(), - simulator.problem().name(), + fileName_ = restartFileName_(simulator, simulator.time()); + // open output file and write magic cookie outStream_.open(fileName_.c_str()); outStream_.precision(20); @@ -172,7 +183,7 @@ public: template void deserializeBegin(Simulator& simulator, Scalar t) { - fileName_ = restartFileName_(simulator.gridView(), simulator.problem().outputDir(), simulator.problem().name(), t); + fileName_ = restartFileName_(simulator, t); // open input file and read magic cookie inStream_.open(fileName_.c_str()); From 73d32fbf502760c1e969725d26c769faa6d0f4b4 Mon Sep 17 00:00:00 2001 From: hnil Date: Mon, 19 Feb 2018 13:58:53 +0100 Subject: [PATCH 11/38] added extra terminal output for deserialization --- ewoms/common/simulator.hh | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/ewoms/common/simulator.hh b/ewoms/common/simulator.hh index 1c19468ad6..bac9c53037 100644 --- a/ewoms/common/simulator.hh +++ b/ewoms/common/simulator.hh @@ -817,6 +817,10 @@ public: typedef Ewoms::Restart Restarter; Restarter res; res.deserializeBegin(*this, t); + if (gridView().comm().rank() == 0) + std::cout << "Deserialize file '" << res.fileName() << "'" + << ", next time step size: " << timeStepSize() + << "\n" << std::flush; this->deserialize(res); problem_->deserialize(res, false); model_->deserialize(res); From 7cc9c6b2359182aa63247e94d6569f907223a800 Mon Sep 17 00:00:00 2001 From: hnil Date: Fri, 25 May 2018 09:59:10 +0200 Subject: [PATCH 12/38] changed order of primary variables --- ewoms/models/blackoil/blackoilindices.hh | 4 ++-- ewoms/models/blackoil/blackoiltwophaseindices.hh | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ewoms/models/blackoil/blackoilindices.hh b/ewoms/models/blackoil/blackoilindices.hh index 9e45abc903..5bbb5de7fe 100644 --- a/ewoms/models/blackoil/blackoilindices.hh +++ b/ewoms/models/blackoil/blackoilindices.hh @@ -81,10 +81,10 @@ public: //////// //! The index of the water saturation - static const int waterSaturationIdx = PVOffset + 0; + static const int waterSaturationIdx = PVOffset + 1; //! Index of the oil pressure in a vector of primary variables - static const int pressureSwitchIdx = PVOffset + 1; + static const int pressureSwitchIdx = PVOffset + 0; /*! * \brief Index of the switching variable which determines the composition of the diff --git a/ewoms/models/blackoil/blackoiltwophaseindices.hh b/ewoms/models/blackoil/blackoiltwophaseindices.hh index a923c6c1bb..45444af34b 100644 --- a/ewoms/models/blackoil/blackoiltwophaseindices.hh +++ b/ewoms/models/blackoil/blackoiltwophaseindices.hh @@ -76,10 +76,10 @@ public: ////////////////////////////// //! The index of the water saturation. For two-phase oil gas models this is disabled. - static const int waterSaturationIdx = waterEnabled ? PVOffset + 0 : -10000; + static const int waterSaturationIdx = waterEnabled ? PVOffset + 1 : -10000; //! Index of the oil pressure in a vector of primary variables - static const int pressureSwitchIdx = waterEnabled ? PVOffset + 1 : PVOffset + 0; + static const int pressureSwitchIdx = PVOffset + 0; /*! * \brief Index of the switching variable which determines the composition of the From 10f7fca1ad2805405e0ec1a953f4ae4ada432214 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 21 Aug 2018 15:46:09 +0200 Subject: [PATCH 13/38] Update file after failed rebase. --- ewoms/disc/common/fvbaselocalresidual.hh | 33 ------------------- .../blackoil/blackoilintensivequantities.hh | 2 +- 2 files changed, 1 insertion(+), 34 deletions(-) diff --git a/ewoms/disc/common/fvbaselocalresidual.hh b/ewoms/disc/common/fvbaselocalresidual.hh index 7471ac1cb0..d9b4f526f6 100644 --- a/ewoms/disc/common/fvbaselocalresidual.hh +++ b/ewoms/disc/common/fvbaselocalresidual.hh @@ -558,38 +558,6 @@ protected: } } else { - // if the mass storage at the beginning of the time step is not cached, - // we re-calculate it from scratch. - tmp2 = 0.0; - asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/1); - Opm::Valgrind::CheckDefined(tmp2); - } - - if (elemCtx.enableStorageCache()) { - const auto& model = elemCtx.model(); - unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - if (model.newtonMethod().numIterations() == 0 && - !elemCtx.haveStashedIntensiveQuantities()) - { - // if the storage term is cached and we're in the first iteration of - // the time step, update the cache of the storage term (this assumes - // that the initial guess for the solution at the end of the time - // step is the same as the solution at the beginning of the time - // step. This is usually true, but some fancy preprocessing scheme - // might invalidate that assumption.) - for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) - tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]); - Opm::Valgrind::CheckDefined(tmp2); - - model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2); - } - else { - // if the storage term is cached and we're not looking at the first - // iteration of the time step, we take the cached data. - tmp2 = model.cachedStorage(globalDofIdx, /*timeIdx=*/1); - Opm::Valgrind::CheckDefined(tmp2); - } - } else { // if the mass storage at the beginning of the time step is not cached, // we re-calculate it from scratch. // if storage cache not enables we always use derivatives @@ -620,7 +588,6 @@ protected: } } - Opm::Valgrind::CheckDefined(residual[dofIdx]); // deal with the source term diff --git a/ewoms/models/blackoil/blackoilintensivequantities.hh b/ewoms/models/blackoil/blackoilintensivequantities.hh index 3e1ee67e0d..212879a604 100644 --- a/ewoms/models/blackoil/blackoilintensivequantities.hh +++ b/ewoms/models/blackoil/blackoilintensivequantities.hh @@ -329,7 +329,7 @@ public: asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx, focustimeidx);// strictly do not need focus time since no makeEvaluation is done asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx, focustimeidx); - asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache, focustimeidx); + asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache); // update the quantities which are required by the chosen // velocity model From 89b6ec74f417b55062da6c01f639b68c4b6ec4eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 Aug 2018 10:20:07 +0200 Subject: [PATCH 14/38] Adding focusTimeIdx parameter to energy module where necessary. --- ewoms/models/blackoil/blackoilenergymodules.hh | 8 +++++--- ewoms/models/blackoil/blackoilintensivequantities.hh | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/ewoms/models/blackoil/blackoilenergymodules.hh b/ewoms/models/blackoil/blackoilenergymodules.hh index 63af772210..c58168ec36 100644 --- a/ewoms/models/blackoil/blackoilenergymodules.hh +++ b/ewoms/models/blackoil/blackoilenergymodules.hh @@ -351,13 +351,14 @@ public: */ void updateTemperature_(const ElementContext& elemCtx, unsigned dofIdx, - unsigned timeIdx) + unsigned timeIdx, + unsigned focusTimeIdx) { auto& fs = asImp_().fluidState_; const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx); // set temperature - fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx)); + fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx, focusTimeIdx)); } /*! @@ -419,7 +420,8 @@ class BlackOilEnergyIntensiveQuantities public: void updateTemperature_(const ElementContext& elemCtx, unsigned dofIdx, - unsigned timeIdx) + unsigned timeIdx, + unsigned focusTimeIdx) { if (enableTemperature) { // even if energy is conserved, the temperature can vary over the spatial diff --git a/ewoms/models/blackoil/blackoilintensivequantities.hh b/ewoms/models/blackoil/blackoilintensivequantities.hh index 212879a604..f3d3425782 100644 --- a/ewoms/models/blackoil/blackoilintensivequantities.hh +++ b/ewoms/models/blackoil/blackoilintensivequantities.hh @@ -114,7 +114,7 @@ public: const auto& problem = elemCtx.problem(); const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx); - asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx); + asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx, focustimeidx); unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx); unsigned pvtRegionIdx = priVars.pvtRegionIndex(); From e070472c475ff1095ea074116cbffa54286443e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2018 14:32:41 +0200 Subject: [PATCH 15/38] Fix fallout from rebasing. --- ewoms/io/restart.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ewoms/io/restart.hh b/ewoms/io/restart.hh index 82cee57281..67b9fe7d11 100644 --- a/ewoms/io/restart.hh +++ b/ewoms/io/restart.hh @@ -70,7 +70,7 @@ class Restart /*! * \brief Return the restart file name. */ - template + template static const std::string restartFileName_(const Simulator& simulator, Scalar t) { // Directory. @@ -88,7 +88,7 @@ class Restart } // Filename. - int rank = gridView.comm().rank(); + int rank = simulator.gridView().comm().rank(); std::string simName = simulator.problem().name(); std::ostringstream oss; oss << simName << "_time=" << t << "_rank=" << rank << ".ers"; From 0f36a7657269467ca3f19fb527be196ac417a90c Mon Sep 17 00:00:00 2001 From: hnil Date: Mon, 5 Nov 2018 23:27:00 +0100 Subject: [PATCH 16/38] tried to correct with respect to coment in pull request 371 --- ebos/eclproblem.hh | 6 ++-- ewoms/common/simulator.hh | 2 +- ewoms/disc/common/fvbaseelementcontext.hh | 15 ++++++++++ ewoms/disc/common/fvbaselinearizer.hh | 9 +++--- ewoms/disc/common/fvbaselocalresidual.hh | 12 ++++---- ewoms/disc/common/fvbaseprimaryvariables.hh | 4 +-- .../blackoil/blackoilintensivequantities.hh | 28 +++++++++---------- ewoms/models/blackoil/blackoilmodel.hh | 2 +- .../models/blackoil/blackoilpolymermodules.hh | 4 +-- 9 files changed, 49 insertions(+), 33 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index d2e8e669f5..92d88b3dae 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -238,8 +238,8 @@ SET_SCALAR_PROP(EclBaseProblem, NewtonMaxError, 0.1); SET_INT_PROP(EclBaseProblem, NewtonMaxIterations, 14); // set no adjoint variables ad default -NEW_PROP_TAG(numAdjoint); -SET_INT_PROP(EclBaseProblem, numAdjoint, 0); +NEW_PROP_TAG(NumAdjoint); +SET_INT_PROP(EclBaseProblem, NumAdjoint, 0); // also, reduce the target for the "optimum" number of Newton iterations to 6. Note that // this is only relevant if the time step is reduced from the report step size for some @@ -629,7 +629,7 @@ public: void deserialize(Restarter& res, bool isOnRestart) { // reload the current episode/report step from the deck - //beginEpisode(/*isOnRestart=*/true); + // beginEpisode(/*isOnRestart=*/true); // should probably be removed if(isOnRestart){ beginEpisode(/*isOnRestart=*/true); diff --git a/ewoms/common/simulator.hh b/ewoms/common/simulator.hh index bac9c53037..cf2eea343e 100644 --- a/ewoms/common/simulator.hh +++ b/ewoms/common/simulator.hh @@ -554,7 +554,7 @@ public: if (verbose_) std::cout << "Deserialize from file '" << res.fileName() << "'\n" << std::flush; this->deserialize(res); - problem_->deserialize(res, true);// to things for restart + problem_->deserialize(res, true);// true make this code call beginEpisode model_->deserialize(res); res.deserializeEnd(); if (verbose_) diff --git a/ewoms/disc/common/fvbaseelementcontext.hh b/ewoms/disc/common/fvbaseelementcontext.hh index ddc7a218e7..ef7d23a13e 100644 --- a/ewoms/disc/common/fvbaseelementcontext.hh +++ b/ewoms/disc/common/fvbaseelementcontext.hh @@ -262,8 +262,18 @@ public: */ void setFocusDofIndex(unsigned dofIdx) { focusDofIdx_ = dofIdx; } + + /*! + * \brief Sets the time index on which the simulator is currently "focused" on + * + * I.e., in the case of automatic differentiation, all derivatives are with regard to + * the primary variables of that time index. Only "primary" DOFs can be + * focused on. + */ + void setFocusTimeIndex(unsigned timeIdx) { focusTimeIdx_ = timeIdx; } + /*! * \brief Returns the degree of freedom on which the simulator is currently "focused" on * @@ -272,6 +282,11 @@ public: unsigned focusDofIndex() const { return focusDofIdx_; } + /*! + * \brief Returns the time index on which the simulator is currently "focused" on + * + * \copydetails setFocusDof() + */ unsigned focusTimeIndex() const { return focusTimeIdx_; } diff --git a/ewoms/disc/common/fvbaselinearizer.hh b/ewoms/disc/common/fvbaselinearizer.hh index 5aff556588..0dbfcb1005 100644 --- a/ewoms/disc/common/fvbaselinearizer.hh +++ b/ewoms/disc/common/fvbaselinearizer.hh @@ -159,14 +159,15 @@ public: } /*! - * \brief Linearize the full system of non-linear equations. + * \brief Linearize the full system of non-linear equations about focusTimeIdx + * defualt is 0 i.e. current time normally time to be found. * * This means the spatial domain plus all auxiliary equations. */ - void linearize(unsigned focustimeIndex) + void linearize(unsigned focusTimeIdx = 0) { - linearizeDomain(focustimeIndex); - linearizeAuxiliaryEquations(focustimeIndex); + linearizeDomain(focusTimeIdx); + linearizeAuxiliaryEquations(focusTimeIdx); } /*! diff --git a/ewoms/disc/common/fvbaselocalresidual.hh b/ewoms/disc/common/fvbaselocalresidual.hh index d9b4f526f6..bf7dab525e 100644 --- a/ewoms/disc/common/fvbaselocalresidual.hh +++ b/ewoms/disc/common/fvbaselocalresidual.hh @@ -492,7 +492,7 @@ protected: { EvalVector tmp; EqVector tmp2; - EvalVector tmp2_der;// this is added to support calculation of storage term derivative focusTime!=0 for need for adjoint: always used if cachedStororage term si false + EvalVector tmp2Der;// this is added to support calculation of storage term derivative focusTime!=0 for need for adjoint: always used if cachedStororage term si false RateVector sourceRate; tmp = 0.0; @@ -561,20 +561,20 @@ protected: // if the mass storage at the beginning of the time step is not cached, // we re-calculate it from scratch. // if storage cache not enables we always use derivatives - tmp2_der = 0.0; - asImp_().computeStorage(tmp2_der, elemCtx, dofIdx, /*timeIdx=*/1); - Opm::Valgrind::CheckDefined(tmp2_der); + tmp2Der = 0.0; + asImp_().computeStorage(tmp2Der, elemCtx, dofIdx, /*timeIdx=*/1); + Opm::Valgrind::CheckDefined(tmp2Der); } if( !elemCtx.enableStorageCache() ){ // Use the implicit Euler time discretization for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { if(elemCtx.focusTimeIndex()==0){ - tmp2_der[eqIdx] .clearDerivatives();// remove derivatives for the time index + tmp2Der[eqIdx] .clearDerivatives();// remove derivatives for the time index }else{ tmp[eqIdx] .clearDerivatives(); } - tmp[eqIdx] -= tmp2_der[eqIdx]; + tmp[eqIdx] -= tmp2Der[eqIdx]; tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize(); residual[dofIdx][eqIdx] += tmp[eqIdx]; } diff --git a/ewoms/disc/common/fvbaseprimaryvariables.hh b/ewoms/disc/common/fvbaseprimaryvariables.hh index 01962ea4a2..0f9525e7f3 100644 --- a/ewoms/disc/common/fvbaseprimaryvariables.hh +++ b/ewoms/disc/common/fvbaseprimaryvariables.hh @@ -87,9 +87,9 @@ public: * it represents the a constant f = x_i. (the difference is that in the first case, * the derivative w.r.t. x_i is 1, while it is 0 in the second case. */ - Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx,unsigned focustimeidx) const + Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx,unsigned focusTimeIdx) const { - if (timeIdx == focustimeidx) + if (timeIdx == focusTimeIdx) return Toolbox::createVariable((*this)[varIdx], varIdx); else return Toolbox::createConstant((*this)[varIdx]); diff --git a/ewoms/models/blackoil/blackoilintensivequantities.hh b/ewoms/models/blackoil/blackoilintensivequantities.hh index f3d3425782..731e505273 100644 --- a/ewoms/models/blackoil/blackoilintensivequantities.hh +++ b/ewoms/models/blackoil/blackoilintensivequantities.hh @@ -109,12 +109,12 @@ public: */ void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focustimeidx) { - ParentType::update(elemCtx, dofIdx, timeIdx);// this only gives extrusion factor + ParentType::update(elemCtx, dofIdx, timeIdx); const auto& problem = elemCtx.problem(); const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx); - asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx, focustimeidx); + asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx, focusTimeIdx); unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx); unsigned pvtRegionIdx = priVars.pvtRegionIndex(); @@ -123,21 +123,21 @@ public: // extract the water and the gas saturations for convenience Evaluation Sw = 0.0; if (waterEnabled) - Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx, focustimeidx); + Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx, focusTimeIdx); Evaluation Sg = 0.0; if (compositionSwitchEnabled) { if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) // -> threephase case - Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focustimeidx); + Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focusTimeIdx); else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) { // -> gas-water case Sg = 1.0 - Sw; // deal with solvent if (enableSolvent) - Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx, focustimeidx); + Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx, focusTimeIdx); } else { @@ -154,13 +154,13 @@ public: // deal with solvent if (enableSolvent) - So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx, focustimeidx); + So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx, focusTimeIdx); fluidState_.setSaturation(waterPhaseIdx, Sw); fluidState_.setSaturation(gasPhaseIdx, Sg); fluidState_.setSaturation(oilPhaseIdx, So); - asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx, focustimeidx); + asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx, focusTimeIdx); // now we compute all phase pressures Evaluation pC[numPhases]; @@ -169,13 +169,13 @@ public: //oil is the reference phase for pressure if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) { - const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focustimeidx); + const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focusTimeIdx); for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx])); } else { - const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focustimeidx); + const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, focusTimeIdx); for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx])); } @@ -186,7 +186,7 @@ public: Opm::Valgrind::CheckDefined(mobility_); // update the Saturation functions for the blackoil solvent module. - asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx, focustimeidx); + asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx, focusTimeIdx); Scalar SoMax = elemCtx.problem().maxOilSaturation(globalSpaceIdx); @@ -224,7 +224,7 @@ public: Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(globalSpaceIdx); // oil phase, we can directly set the composition of the oil phase - const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focustimeidx); + const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focusTimeIdx); fluidState_.setRs(Opm::min(RsMax, Rs)); if (FluidSystem::enableVaporizedOil()) { @@ -245,7 +245,7 @@ public: else { assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv); - const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focustimeidx); + const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx, focusTimeIdx); fluidState_.setRv(Rv); if (FluidSystem::enableDissolvedGas()) { @@ -327,8 +327,8 @@ public: porosity_ *= 1.0 + x + 0.5*x*x; } - asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx, focustimeidx);// strictly do not need focus time since no makeEvaluation is done - asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx, focustimeidx); + asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx, focusTimeIdx); + asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx, focusTimeIdx); asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache); // update the quantities which are required by the chosen diff --git a/ewoms/models/blackoil/blackoilmodel.hh b/ewoms/models/blackoil/blackoilmodel.hh index 7bc36e1821..3798c83eea 100644 --- a/ewoms/models/blackoil/blackoilmodel.hh +++ b/ewoms/models/blackoil/blackoilmodel.hh @@ -503,7 +503,7 @@ public: } } - // this is restart spesific behavoir and to avoid this one have to hack outside for adjiont + // this is restart spesific this->solution(/*timeIdx=*/1) = this->solution(/*timeIdx=*/0); } diff --git a/ewoms/models/blackoil/blackoilpolymermodules.hh b/ewoms/models/blackoil/blackoilpolymermodules.hh index 84114b9217..655b0dcede 100644 --- a/ewoms/models/blackoil/blackoilpolymermodules.hh +++ b/ewoms/models/blackoil/blackoilpolymermodules.hh @@ -865,7 +865,7 @@ public: unsigned focustime) { const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx); - polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx, focustime); + polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx, focusTimeIdx); const Scalar cmax = PolymerModule::plymaxMaxConcentration(elemCtx, dofIdx, timeIdx); // permeability reduction due to polymer @@ -954,7 +954,7 @@ public: void polymerPropertiesUpdate_(const ElementContext& elemCtx OPM_UNUSED, unsigned scvIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED, - unsigned focustimeIdx OPM_UNUSED) + unsigned focusTimeIdxIdx OPM_UNUSED) { } From 86a2ec1827c988eb9da96cdf1c0c01cabc21e74a Mon Sep 17 00:00:00 2001 From: hnil Date: Fri, 30 Nov 2018 23:51:29 +0100 Subject: [PATCH 17/38] corrected camelcase --- ewoms/models/blackoil/blackoilintensivequantities.hh | 2 +- ewoms/models/blackoil/blackoilpolymermodules.hh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ewoms/models/blackoil/blackoilintensivequantities.hh b/ewoms/models/blackoil/blackoilintensivequantities.hh index 731e505273..48527214f9 100644 --- a/ewoms/models/blackoil/blackoilintensivequantities.hh +++ b/ewoms/models/blackoil/blackoilintensivequantities.hh @@ -107,7 +107,7 @@ public: /*! * \copydoc IntensiveQuantities::update */ - void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focustimeidx) + void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focusTimeIdx) { ParentType::update(elemCtx, dofIdx, timeIdx); diff --git a/ewoms/models/blackoil/blackoilpolymermodules.hh b/ewoms/models/blackoil/blackoilpolymermodules.hh index 655b0dcede..6f3cca2e02 100644 --- a/ewoms/models/blackoil/blackoilpolymermodules.hh +++ b/ewoms/models/blackoil/blackoilpolymermodules.hh @@ -862,7 +862,7 @@ public: void polymerPropertiesUpdate_(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, - unsigned focustime) + unsigned focusTimeIdx) { const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx); polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx, focusTimeIdx); From 2c81e37287bd23baa56090c01b929161df704904 Mon Sep 17 00:00:00 2001 From: hnil Date: Sun, 2 Dec 2018 21:28:09 +0100 Subject: [PATCH 18/38] Enabled reseting of Storage cache: it is critical that all is deleted: maybe a setEnableStoargeCache should be reintroduced in the simulator --- ewoms/disc/common/fvbasediscretization.hh | 9 +++++++++ ewoms/disc/common/fvbaseelementcontext.hh | 2 +- ewoms/disc/common/fvbaselinearizer.hh | 7 +++++++ ewoms/disc/common/fvbaselocalresidual.hh | 2 +- 4 files changed, 18 insertions(+), 2 deletions(-) diff --git a/ewoms/disc/common/fvbasediscretization.hh b/ewoms/disc/common/fvbasediscretization.hh index 9178aea2fe..41568d754a 100644 --- a/ewoms/disc/common/fvbasediscretization.hh +++ b/ewoms/disc/common/fvbasediscretization.hh @@ -762,6 +762,15 @@ public: bool enableStorageCache() const { return enableStorageCache_; } + /*! + * \brief Set the value of enable storage cache + * + * Be aware that calling the *CachedStorage() methods if the storage cache is + * disabled will crash the program. + */ + void setEnableStorageCache(bool enableStorageCache) + { enableStorageCache_= enableStorageCache; } + /*! * \brief Retrieve an entry of the cache for the storage term. * diff --git a/ewoms/disc/common/fvbaseelementcontext.hh b/ewoms/disc/common/fvbaseelementcontext.hh index ef7d23a13e..c5bb224867 100644 --- a/ewoms/disc/common/fvbaseelementcontext.hh +++ b/ewoms/disc/common/fvbaseelementcontext.hh @@ -97,7 +97,7 @@ public: { // remember the simulator object simulatorPtr_ = &simulator; - enableStorageCache_ = EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache); + enableStorageCache_ = simulator.model().enableStorageCache();//EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache); stashedDofIdx_ = -1; focusDofIdx_ = -1; focusTimeIdx_ = 0; diff --git a/ewoms/disc/common/fvbaselinearizer.hh b/ewoms/disc/common/fvbaselinearizer.hh index 0dbfcb1005..7b4b438422 100644 --- a/ewoms/disc/common/fvbaselinearizer.hh +++ b/ewoms/disc/common/fvbaselinearizer.hh @@ -143,6 +143,12 @@ public: simulatorPtr_ = &simulator; delete matrix_; // <- note that this even works for nullpointers! matrix_ = 0; + auto it = elementCtx_.begin(); + const auto& endIt = elementCtx_.end(); + for (; it != endIt; ++it){ + delete *it; + } + elementCtx_.resize(0); } /*! @@ -327,6 +333,7 @@ private: elementCtx_[threadId] = new ElementContext(simulator_()); } + // Construct the BCRS matrix for the Jacobian of the residual function void createMatrix_() { diff --git a/ewoms/disc/common/fvbaselocalresidual.hh b/ewoms/disc/common/fvbaselocalresidual.hh index bf7dab525e..56b2e598d4 100644 --- a/ewoms/disc/common/fvbaselocalresidual.hh +++ b/ewoms/disc/common/fvbaselocalresidual.hh @@ -566,7 +566,7 @@ protected: Opm::Valgrind::CheckDefined(tmp2Der); } - if( !elemCtx.enableStorageCache() ){ + if( !elemCtx.enableStorageCache() or elemCtx.focusTimeIndex()>0){ // Use the implicit Euler time discretization for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { if(elemCtx.focusTimeIndex()==0){ From 89295d9a4fd2763ad205aad43c94db0f5a9eba5f Mon Sep 17 00:00:00 2001 From: hnil Date: Fri, 7 Dec 2018 17:05:20 +0100 Subject: [PATCH 19/38] added possibility to avoid removing auxModule. Fis may be wrong if topology of wells is changing: --- ebos/eclproblem.hh | 2 +- ebos/eclwellmanager.hh | 6 +++--- ewoms/common/simulator.hh | 9 ++++++--- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 92d88b3dae..a8578de7d0 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -635,7 +635,7 @@ public: beginEpisode(/*isOnRestart=*/true); } // deserialize the wells - wellModel_.deserialize(res); + wellModel_.deserialize(res, isOnRestart); } /*! diff --git a/ebos/eclwellmanager.hh b/ebos/eclwellmanager.hh index a78f5d1d34..40daa1be7d 100644 --- a/ebos/eclwellmanager.hh +++ b/ebos/eclwellmanager.hh @@ -577,10 +577,10 @@ public: * It is the inverse of the serialize() method. */ template - void deserialize(Restarter& res OPM_UNUSED) + void deserialize(Restarter& res, bool wasRestarted=true) { // initialize the wells for the current episode - beginEpisode(simulator_.vanguard().eclState(), simulator_.vanguard().schedule(), /*wasRestarted=*/true); + beginEpisode(simulator_.vanguard().eclState(), simulator_.vanguard().schedule(), wasRestarted); } /*! @@ -689,7 +689,7 @@ protected: } } - void computeWellConnectionsMap_(unsigned reportStepIdx OPM_UNUSED, WellConnectionsMap& cartesianIdxToConnectionMap) + void computeWellConnectionsMap_(unsigned reportStepIdx, WellConnectionsMap& cartesianIdxToConnectionMap) { const auto& deckSchedule = simulator_.vanguard().schedule(); diff --git a/ewoms/common/simulator.hh b/ewoms/common/simulator.hh index cf2eea343e..d4eb7770b8 100644 --- a/ewoms/common/simulator.hh +++ b/ewoms/common/simulator.hh @@ -812,7 +812,7 @@ public: res.serializeEnd(); } - void deserializeAll(Scalar t) + void deserializeAll(Scalar t,bool only_reservoir) { typedef Ewoms::Restart Restarter; Restarter res; @@ -822,8 +822,11 @@ public: << ", next time step size: " << timeStepSize() << "\n" << std::flush; this->deserialize(res); - problem_->deserialize(res, false); - model_->deserialize(res); + if( not(only_reservoir) ){ + problem_->deserialize(res, false); + } + model_->deserialize(res); + res.deserializeEnd(); } /*! From af756b7f76077f9cebecbbec5bd26d0dbb171d5e Mon Sep 17 00:00:00 2001 From: hnil Date: Fri, 7 Dec 2018 21:37:39 +0100 Subject: [PATCH 20/38] code comples after merge --- ebos/eclproblem.hh | 16 ++++++++-------- ewoms/disc/common/fvbaselinearizer.hh | 2 +- ewoms/nonlinear/newtonmethod.hh | 15 ++++----------- 3 files changed, 13 insertions(+), 20 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index c04dc70428..b14fb41fd0 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -618,9 +618,9 @@ public: beginEpisode(/*isOnRestart=*/true); } // deserialize the wells - wellModel_.deserialize(res, isOnRestart); + //wellModel_.deserialize(res, isOnRestart); // deserialize the aquifer - aquiferModel_.deserialize(res); + //aquiferModel_.deserialize(res); } /*! @@ -629,12 +629,12 @@ public: * * The file format used here is ad-hoc. */ - template - void serialize(Restarter& res) - { - wellModel_.serialize(res); - aquiferModel_.serialize(res); - } + // template + // void serialize(Restarter& res) + // { + // wellModel_.serialize(res); + // //aquiferModel_.serialize(res); + // } /*! * \brief Called by the simulator before an episode begins. diff --git a/ewoms/disc/common/fvbaselinearizer.hh b/ewoms/disc/common/fvbaselinearizer.hh index 4923d094be..9a091e961b 100644 --- a/ewoms/disc/common/fvbaselinearizer.hh +++ b/ewoms/disc/common/fvbaselinearizer.hh @@ -246,7 +246,7 @@ public: for (unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx) { bool succeeded = true; try { - model.auxiliaryModule(auxModIdx)->linearize(*jacobian_, residual_, focusTimeIdx); + model.auxiliaryModule(auxModIdx)->linearize(*jacobian_, residual_);//, focusTimeIdx); } catch (const std::exception& e) { succeeded = false; diff --git a/ewoms/nonlinear/newtonmethod.hh b/ewoms/nonlinear/newtonmethod.hh index fb3dd00811..f7e17a151d 100644 --- a/ewoms/nonlinear/newtonmethod.hh +++ b/ewoms/nonlinear/newtonmethod.hh @@ -641,23 +641,16 @@ protected: * \brief Linearize the global non-linear system of equations associated with the * spatial domain. */ - void linearizeDomain_() -<<<<<<< HEAD + void linearizeDomain_(int focusTimeIdx) { - model().linearizer().linearizeDomain(); + model().linearizer().linearizeDomain(focusTimeIdx); } - void linearizeAuxiliaryEquations_() + void linearizeAuxiliaryEquations_(int focusTimeIdx) { - model().linearizer().linearizeAuxiliaryEquations(); + model().linearizer().linearizeAuxiliaryEquations(focusTimeIdx); model().linearizer().finalize(); } -======= - { model().linearizer().linearizeDomain(/*focus time*/ 0); } - - void linearizeAuxiliaryEquations_() - { model().linearizer().linearizeAuxiliaryEquations(/*focus time*/ 0); } ->>>>>>> adjoint_rebased void preSolve_(const SolutionVector& currentSolution OPM_UNUSED, const GlobalEqVector& currentResidual) From a73c2b9afdb25f16bead7154b24f45dfb2a3adfa Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 11 Dec 2018 15:40:17 +0100 Subject: [PATCH 21/38] working version --- ebos/eclproblem.hh | 12 ++++++------ ewoms/disc/common/fvbaselocalresidual.hh | 6 ++++-- ewoms/linear/matrixblock.hh | 2 +- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index b14fb41fd0..a666656176 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -629,12 +629,12 @@ public: * * The file format used here is ad-hoc. */ - // template - // void serialize(Restarter& res) - // { - // wellModel_.serialize(res); - // //aquiferModel_.serialize(res); - // } + template + void serialize(Restarter& res) + { + //wellModel_.serialize(res); + //aquiferModel_.serialize(res); + } /*! * \brief Called by the simulator before an episode begins. diff --git a/ewoms/disc/common/fvbaselocalresidual.hh b/ewoms/disc/common/fvbaselocalresidual.hh index 56b2e598d4..6b19f331b5 100644 --- a/ewoms/disc/common/fvbaselocalresidual.hh +++ b/ewoms/disc/common/fvbaselocalresidual.hh @@ -591,8 +591,10 @@ protected: Opm::Valgrind::CheckDefined(residual[dofIdx]); // deal with the source term - asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0); - + //assumes sourceterm is independent of prevois state + if(elemCtx.focusTimeIndex()==0){ + asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0); + } // if the model uses extensive quantities in its storage term, and we use // automatic differention and current DOF is also not the one we currently // focus on, the storage term does not need any derivatives! diff --git a/ewoms/linear/matrixblock.hh b/ewoms/linear/matrixblock.hh index 647a0aab8d..8eb44a0da5 100644 --- a/ewoms/linear/matrixblock.hh +++ b/ewoms/linear/matrixblock.hh @@ -29,7 +29,7 @@ #include #include #include - +#include #include namespace Ewoms { From 8db4520187c0338bfcbe2e1cadd150bddb2e6779 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 12 Dec 2018 09:11:49 +0100 Subject: [PATCH 22/38] added spesialization need for getting matrix market work with Ewoms::Block --- ewoms/linear/matrixmarket_ewoms.hh | 198 +++++++++++++++++++++++++++++ 1 file changed, 198 insertions(+) create mode 100644 ewoms/linear/matrixmarket_ewoms.hh diff --git a/ewoms/linear/matrixmarket_ewoms.hh b/ewoms/linear/matrixmarket_ewoms.hh new file mode 100644 index 0000000000..915eb08523 --- /dev/null +++ b/ewoms/linear/matrixmarket_ewoms.hh @@ -0,0 +1,198 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef MATRIXMARKET_EWOMS_HH +#define MATRIXMARKET_EWOMS_HH +#include +#include +namespace Dune +{ + + + namespace MatrixMarketImpl + { + + + template + struct mm_header_printer,A> > + { + static void print(std::ostream& os) + { + os<<"%%MatrixMarket matrix coordinate "; + os<::str()<<" general"< + struct mm_header_printer > + { + static void print(std::ostream& os) + { + os<<"%%MatrixMarket matrix array "; + os<::str()<<" general"< + struct mm_block_structure_header,A> > + { + typedef BCRSMatrix,A> M; + + static void print(std::ostream& os, const M&) + { + os<<"% ISTL_STRUCT blocked "; + os< + struct mm_block_structure_header > + { + typedef Ewoms::MatrixBlock M; + + static void print(std::ostream& os, const M& m) + {} + }; + + + + template + void readSparseEntries(Dune::BCRSMatrix,A>& matrix, + std::istream& file, std::size_t entries, + const MMHeader& mmHeader, const D&) + { + typedef Dune::BCRSMatrix,A> Matrix; + // First path + // store entries together with column index in a speparate + // data structure + std::vector > > rows(matrix.N()*brows); + + for(; entries>0; --entries) { + std::size_t row; + IndexData data; + skipComments(file); + file>>row; + --row; // Index was 1 based. + assert(row/bcols>data; + assert(data.index/bcols >::const_iterator Siter; + for(Siter siter=rows[brow].begin(), send=rows[brow].end(); + siter != send; ++siter, ++nnz) + iter.insert(siter->index/bcols); + } + } + + //Set the matrix values + matrix=0; + + MatrixValuesSetter Setter; + + Setter(rows, matrix); + } + } // end namespace MatrixMarketImpl + + + template + void readMatrixMarket(Dune::BCRSMatrix,A>& matrix, + std::istream& istr) + { + using namespace MatrixMarketImpl; + + MMHeader header; + if(!readMatrixMarketBanner(istr, header)) { + std::cerr << "First line was not a correct Matrix Market banner. Using default:\n" + << "%%MatrixMarket matrix coordinate real general"<> rows; + + if(lineFeed(istr)) + throw MatrixMarketFormatError(); + istr >> cols; + + if(lineFeed(istr)) + throw MatrixMarketFormatError(); + + istr >>entries; + + std::size_t nnz, blockrows, blockcols; + + std::tie(blockrows, blockcols, nnz) = calculateNNZ(rows, cols, entries, header); + + istr.ignore(std::numeric_limits::max(),'\n'); + + + matrix.setSize(blockrows, blockcols); + matrix.setBuildMode(Dune::BCRSMatrix,A>::row_wise); + + if(header.type==array_type) + DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!"); + + readSparseEntries(matrix, istr, entries, header, NumericWrapper()); + } + + + + template + struct mm_multipliers,A> > + { + enum { + rows = i, + cols = j + }; + }; + + template + void mm_print_entry(const Ewoms::MatrixBlock& entry, + typename Ewoms::MatrixBlock::size_type rowidx, + typename Ewoms::MatrixBlock::size_type colidx, + std::ostream& ostr) + { + typedef typename Ewoms::MatrixBlock::const_iterator riterator; + typedef typename Ewoms::MatrixBlock::ConstColIterator citerator; + for(riterator row=entry.begin(); row != entry.end(); ++row, ++rowidx) { + int coli=colidx; + for(citerator col = row->begin(); col != row->end(); ++col, ++coli) + ostr<< rowidx<<" "< Date: Wed, 12 Dec 2018 10:46:11 +0100 Subject: [PATCH 23/38] introduced hack to avoid calling any source terms in the evaluating with focusTimeIdex=1 --- ewoms/disc/common/fvbaselocalresidual.hh | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/ewoms/disc/common/fvbaselocalresidual.hh b/ewoms/disc/common/fvbaselocalresidual.hh index 6b19f331b5..851b8a67cb 100644 --- a/ewoms/disc/common/fvbaselocalresidual.hh +++ b/ewoms/disc/common/fvbaselocalresidual.hh @@ -163,16 +163,17 @@ public: assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0)); residual = 0.0; - + if(elemCtx.focusTimeIndex()==0){ // evaluate the flux terms - asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0); - + asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0); + } // evaluate the storage and the source terms - asImp_().evalVolumeTerms_(residual, elemCtx); + asImp_().evalVolumeTerms_(residual, elemCtx); + if(elemCtx.focusTimeIndex()==0){ // evaluate the boundary conditions - asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0); - + asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0); + } if (useVolumetricResidual) { // make the residual volume specific (i.e., make it incorrect mass per cubic // meter instead of total mass) From 3d625d339ce61a3f8d620d286b41115c8459921b Mon Sep 17 00:00:00 2001 From: hnil Date: Mon, 17 Dec 2018 20:10:49 +0100 Subject: [PATCH 24/38] fixed compilatino of pure ewoms after introdusing focusTimeIdx --- ewoms/nonlinear/newtonmethod.hh | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/ewoms/nonlinear/newtonmethod.hh b/ewoms/nonlinear/newtonmethod.hh index f7e17a151d..6ec844f566 100644 --- a/ewoms/nonlinear/newtonmethod.hh +++ b/ewoms/nonlinear/newtonmethod.hh @@ -365,9 +365,11 @@ public: << std::flush; } + // + int focusTimeIdx = 0; // do the actual linearization linearizeTimer_.start(); - asImp_().linearizeDomain_(); + asImp_().linearizeDomain_(focusTimeIdx); linearizeTimer_.stop(); // notify the implementation of the successful linearization on order to @@ -380,7 +382,7 @@ public: asImp_().preSolve_(currentSolution, residual); updateTimer_.stop(); - asImp_().linearizeAuxiliaryEquations_(); + asImp_().linearizeAuxiliaryEquations_(focusTimeIdx); if (!asImp_().proceed_()) { if (asImp_().verbose_() && isatty(fileno(stdout))) From b0356c41cbf9ba7a1909305be826acd430388489 Mon Sep 17 00:00:00 2001 From: hnil Date: Fri, 4 Jan 2019 12:54:57 +0100 Subject: [PATCH 25/38] changed use fo time from simulator --- ebos/eclwriter.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index b84be5fbce..e6bf7da643 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -147,7 +147,7 @@ public: */ void writeOutput(bool isSubStep) { - Scalar curTime = simulator_.time() + simulator_.timeStepSize(); + Scalar curTime = simulator_.time();// + simulator_.timeStepSize(); Scalar totalSolverTime = simulator_.executionTimer().realTimeElapsed(); Scalar nextStepSize = simulator_.problem().nextTimeStepSize(); From c5c4d679b6f658b9bea829d48d0a6bc35339e150 Mon Sep 17 00:00:00 2001 From: hnil Date: Fri, 4 Jan 2019 15:03:42 +0100 Subject: [PATCH 26/38] changed black oil indices back for less small changes in tests --- ewoms/models/blackoil/blackoilindices.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ewoms/models/blackoil/blackoilindices.hh b/ewoms/models/blackoil/blackoilindices.hh index 5dbadf74e4..413b2d15d4 100644 --- a/ewoms/models/blackoil/blackoilindices.hh +++ b/ewoms/models/blackoil/blackoilindices.hh @@ -79,10 +79,10 @@ struct BlackOilIndices //////// //! The index of the water saturation - static const int waterSaturationIdx = PVOffset + 1; + static const int waterSaturationIdx = PVOffset + 0; //! Index of the oil pressure in a vector of primary variables - static const int pressureSwitchIdx = PVOffset + 0; + static const int pressureSwitchIdx = PVOffset + 1; /*! * \brief Index of the switching variable which determines the composition of the From 9bc92711b4e34a7069f417c4258a74419fb9ab73 Mon Sep 17 00:00:00 2001 From: hnil Date: Mon, 7 Jan 2019 09:42:05 +0100 Subject: [PATCH 27/38] moved back index change to get test ok --- ewoms/models/blackoil/blackoiltwophaseindices.hh | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ewoms/models/blackoil/blackoiltwophaseindices.hh b/ewoms/models/blackoil/blackoiltwophaseindices.hh index b9adc54868..81ac6ba62e 100644 --- a/ewoms/models/blackoil/blackoiltwophaseindices.hh +++ b/ewoms/models/blackoil/blackoiltwophaseindices.hh @@ -74,10 +74,11 @@ struct BlackOilTwoPhaseIndices ////////////////////////////// //! The index of the water saturation. For two-phase oil gas models this is disabled. - static const int waterSaturationIdx = waterEnabled ? PVOffset + 1 : -10000; + static const int waterSaturationIdx = waterEnabled ? PVOffset + 0 : -10000; //! Index of the oil pressure in a vector of primary variables - static const int pressureSwitchIdx = PVOffset + 0; + static const int pressureSwitchIdx = waterEnabled ? PVOffset + 1 : PVOffset + 0; + //static const int pressureSwitchIdx = PVOffset + 0; /*! * \brief Index of the switching variable which determines the composition of the From ac9ad80cea64e123889d6656fd3884dfc1fd3ce8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 14 Jan 2019 17:16:20 +0100 Subject: [PATCH 28/38] Use default argument to avoid changing API. --- ebos/eclproblem.hh | 2 +- ewoms/common/simulator.hh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 49d4925787..3f103fd0c8 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -643,7 +643,7 @@ public: * \param res The deserializer object */ template - void deserialize(Restarter& res, bool isOnRestart) + void deserialize(Restarter& res, bool isOnRestart = true) { // reload the current episode/report step from the deck // beginEpisode(/*isOnRestart=*/true); diff --git a/ewoms/common/simulator.hh b/ewoms/common/simulator.hh index d4eb7770b8..e3cb773b5b 100644 --- a/ewoms/common/simulator.hh +++ b/ewoms/common/simulator.hh @@ -554,7 +554,7 @@ public: if (verbose_) std::cout << "Deserialize from file '" << res.fileName() << "'\n" << std::flush; this->deserialize(res); - problem_->deserialize(res, true);// true make this code call beginEpisode + problem_->deserialize(res); model_->deserialize(res); res.deserializeEnd(); if (verbose_) From 91993880e3f505c8808684be87a189f3dfb49570 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 14 Jan 2019 17:17:06 +0100 Subject: [PATCH 29/38] Enable compiling also for grids without the name() method. --- ewoms/io/restart.hh | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/ewoms/io/restart.hh b/ewoms/io/restart.hh index 67b9fe7d11..cd22b72332 100644 --- a/ewoms/io/restart.hh +++ b/ewoms/io/restart.hh @@ -35,6 +35,30 @@ namespace Ewoms { + +// This small trio of functions is necessary to provide a generic +// interface to getting grid names, when only some grids provide the +// name() method. +template +auto getGridNameImp(const GridView& gridView, int) + -> decltype(gridView.grid().name()) +{ + return gridView.grid().name(); +}; +template +auto getGridNameImp(const GridView&, long long) + -> std::string +{ + return "GridWithNoNameMethod"; +}; +template +auto getGridName(const GridView& gridView) + -> decltype(getGridNameImp(gridView, 0)) +{ + return getGridNameImp(gridView, 0); +}; + + /*! * \brief Load or save a state of a problem to/from the harddisk. */ @@ -47,7 +71,7 @@ class Restart template static const std::string magicRestartCookie_(const GridView& gridView) { - static const std::string gridName = gridView.grid().name(); + static const std::string gridName = getGridName(gridView); static const int dim = GridView::dimension; int numVertices = gridView.size(dim); From f6a16871a92bcbd29a3cfb6ce451d71a800548d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Jan 2019 14:51:18 +0100 Subject: [PATCH 30/38] Silence shadowing warning. --- ewoms/common/parametersystem.hh | 2 -- 1 file changed, 2 deletions(-) diff --git a/ewoms/common/parametersystem.hh b/ewoms/common/parametersystem.hh index d532d30258..617b5a59b2 100644 --- a/ewoms/common/parametersystem.hh +++ b/ewoms/common/parametersystem.hh @@ -930,7 +930,6 @@ public: check_(Dune::className(), propTagName, paramName); #endif - typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta; if (errorIfNotRegistered) { if (ParamsMeta::registrationOpen()) throw std::runtime_error("Parameters can only checked after _all_ of them have " @@ -1007,7 +1006,6 @@ private: check_(Dune::className(), propTagName, paramName); #endif - typedef typename GET_PROP(TypeTag, ParameterMetaData) ParamsMeta; if (errorIfNotRegistered) { if (ParamsMeta::registrationOpen()) throw std::runtime_error("Parameters can only retieved after _all_ of them have " From 5ad25b8ec4055be9a6558eb825dd36a8b1d4b0fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Jan 2019 13:39:48 +0100 Subject: [PATCH 31/38] Add default argument to preserve existing API. --- ewoms/disc/common/fvbaseprimaryvariables.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ewoms/disc/common/fvbaseprimaryvariables.hh b/ewoms/disc/common/fvbaseprimaryvariables.hh index 0f9525e7f3..521b0aea52 100644 --- a/ewoms/disc/common/fvbaseprimaryvariables.hh +++ b/ewoms/disc/common/fvbaseprimaryvariables.hh @@ -87,7 +87,7 @@ public: * it represents the a constant f = x_i. (the difference is that in the first case, * the derivative w.r.t. x_i is 1, while it is 0 in the second case. */ - Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx,unsigned focusTimeIdx) const + Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx, unsigned focusTimeIdx = 0) const { if (timeIdx == focusTimeIdx) return Toolbox::createVariable((*this)[varIdx], varIdx); From 1ad1bd9775611471929297c129eb3988cb6bd540 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 Jan 2019 14:41:57 +0100 Subject: [PATCH 32/38] Use Toolbox::value() instead of clearDerivatives(). This makes the code compile also when the local Evaluation type is double. --- ewoms/disc/common/fvbaselocalresidual.hh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ewoms/disc/common/fvbaselocalresidual.hh b/ewoms/disc/common/fvbaselocalresidual.hh index 55a3da2604..823059cd1a 100644 --- a/ewoms/disc/common/fvbaselocalresidual.hh +++ b/ewoms/disc/common/fvbaselocalresidual.hh @@ -583,10 +583,10 @@ protected: if( !elemCtx.enableStorageCache() or elemCtx.focusTimeIndex()>0){ // Use the implicit Euler time discretization for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { - if(elemCtx.focusTimeIndex()==0){ - tmp2Der[eqIdx] .clearDerivatives();// remove derivatives for the time index - }else{ - tmp[eqIdx] .clearDerivatives(); + if (elemCtx.focusTimeIndex() == 0) { + tmp2Der[eqIdx] = Toolbox::value(tmp2Der[eqIdx]); + } else { + tmp[eqIdx] = Toolbox::value(tmp[eqIdx]); } tmp[eqIdx] -= tmp2Der[eqIdx]; tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize(); From 31f80b8ec6cd507a4c76448c0b4968fef8eee309 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 Jan 2019 14:43:37 +0100 Subject: [PATCH 33/38] Silence unused variable warning. --- ewoms/disc/common/fvbaselinearizer.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ewoms/disc/common/fvbaselinearizer.hh b/ewoms/disc/common/fvbaselinearizer.hh index 9a091e961b..15c7dbff64 100644 --- a/ewoms/disc/common/fvbaselinearizer.hh +++ b/ewoms/disc/common/fvbaselinearizer.hh @@ -236,7 +236,7 @@ public: * \brief Linearize the part of the non-linear system of equations that is associated * with the spatial domain. */ - void linearizeAuxiliaryEquations(unsigned focusTimeIdx = 0) + void linearizeAuxiliaryEquations(unsigned focusTimeIdx OPM_UNUSED = 0) { // flush possible local caches into matrix structure jacobian_->commit(); From 1348d05690dda59ee379cfee11bd87e1c88a8e80 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 16 Jan 2019 14:43:56 +0100 Subject: [PATCH 34/38] Add focusTimeIdx to all ...IntensiveQuantities classes. --- ewoms/disc/common/fvbaseintensivequantities.hh | 3 ++- ewoms/models/blackoil/blackoilintensivequantities.hh | 2 +- ewoms/models/flash/flashintensivequantities.hh | 4 ++-- ewoms/models/immiscible/immiscibleintensivequantities.hh | 4 ++-- ewoms/models/ncp/ncpintensivequantities.hh | 5 +++-- ewoms/models/pvs/pvsintensivequantities.hh | 4 ++-- ewoms/models/richards/richardsintensivequantities.hh | 4 ++-- 7 files changed, 14 insertions(+), 12 deletions(-) diff --git a/ewoms/disc/common/fvbaseintensivequantities.hh b/ewoms/disc/common/fvbaseintensivequantities.hh index 7ed89f2b6c..1928dfca83 100644 --- a/ewoms/disc/common/fvbaseintensivequantities.hh +++ b/ewoms/disc/common/fvbaseintensivequantities.hh @@ -67,7 +67,8 @@ public: */ void update(const ElementContext& elemCtx, unsigned dofIdx, - unsigned timeIdx) + unsigned timeIdx, + unsigned focusTimeIdx OPM_UNUSED) { extrusionFactor_ = elemCtx.problem().extrusionFactor(elemCtx, dofIdx, timeIdx); } /*! diff --git a/ewoms/models/blackoil/blackoilintensivequantities.hh b/ewoms/models/blackoil/blackoilintensivequantities.hh index 7cd8d923fd..906c653f97 100644 --- a/ewoms/models/blackoil/blackoilintensivequantities.hh +++ b/ewoms/models/blackoil/blackoilintensivequantities.hh @@ -111,7 +111,7 @@ public: */ void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focusTimeIdx) { - ParentType::update(elemCtx, dofIdx, timeIdx); + ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx); const auto& problem = elemCtx.problem(); const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx); diff --git a/ewoms/models/flash/flashintensivequantities.hh b/ewoms/models/flash/flashintensivequantities.hh index cd9e632010..a2cc8a4d7d 100644 --- a/ewoms/models/flash/flashintensivequantities.hh +++ b/ewoms/models/flash/flashintensivequantities.hh @@ -99,9 +99,9 @@ public: /*! * \copydoc IntensiveQuantities::update */ - void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx) + void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focusTimeIdx) { - ParentType::update(elemCtx, dofIdx, timeIdx); + ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx); EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx); const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx); diff --git a/ewoms/models/immiscible/immiscibleintensivequantities.hh b/ewoms/models/immiscible/immiscibleintensivequantities.hh index cad51acf22..468f8fef4e 100644 --- a/ewoms/models/immiscible/immiscibleintensivequantities.hh +++ b/ewoms/models/immiscible/immiscibleintensivequantities.hh @@ -90,9 +90,9 @@ public: /*! * \copydoc IntensiveQuantities::update */ - void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx) + void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focusTimeIdx) { - ParentType::update(elemCtx, dofIdx, timeIdx); + ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx); EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx); // material law parameters diff --git a/ewoms/models/ncp/ncpintensivequantities.hh b/ewoms/models/ncp/ncpintensivequantities.hh index c125e9a43e..dc6ede7fc7 100644 --- a/ewoms/models/ncp/ncpintensivequantities.hh +++ b/ewoms/models/ncp/ncpintensivequantities.hh @@ -101,9 +101,10 @@ public: */ void update(const ElementContext& elemCtx, unsigned dofIdx, - unsigned timeIdx) + unsigned timeIdx, + unsigned focusTimeIdx) { - ParentType::update(elemCtx, dofIdx, timeIdx); + ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx); ParentType::checkDefined(); typename FluidSystem::template ParameterCache paramCache; diff --git a/ewoms/models/pvs/pvsintensivequantities.hh b/ewoms/models/pvs/pvsintensivequantities.hh index 3c64276a58..b62fd5959f 100644 --- a/ewoms/models/pvs/pvsintensivequantities.hh +++ b/ewoms/models/pvs/pvsintensivequantities.hh @@ -105,9 +105,9 @@ public: /*! * \copydoc ImmiscibleIntensiveQuantities::update */ - void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx) + void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focusTimeIdx) { - ParentType::update(elemCtx, dofIdx, timeIdx); + ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx); EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx); const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx); diff --git a/ewoms/models/richards/richardsintensivequantities.hh b/ewoms/models/richards/richardsintensivequantities.hh index 4c696c264a..1cb146e734 100644 --- a/ewoms/models/richards/richardsintensivequantities.hh +++ b/ewoms/models/richards/richardsintensivequantities.hh @@ -84,9 +84,9 @@ public: /*! * \copydoc IntensiveQuantities::update */ - void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx) + void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned focusTimeIdx) { - ParentType::update(elemCtx, dofIdx, timeIdx); + ParentType::update(elemCtx, dofIdx, timeIdx, focusTimeIdx); const auto& T = elemCtx.problem().temperature(elemCtx, dofIdx, timeIdx); fluidState_.setTemperature(T); From c9bad00105dd46a3fb0a32a0fb1d1a9e7ef3599c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 17 Jan 2019 10:06:27 +0100 Subject: [PATCH 35/38] Make timestep-past-episode logic optional in Simulator class. For adjoint backwards run, this logic is wrong. Flow provides its own timesteps anyway, so disabling it for Flow allows all existing ewoms tests to succeed, while allowing adjoint runs. --- ebos/eclproblem.hh | 3 +++ ewoms/common/simulator.hh | 20 +++++++++++--------- ewoms/disc/common/fvbasediscretization.hh | 3 +++ 3 files changed, 17 insertions(+), 9 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 3f103fd0c8..c6f7b31279 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -336,6 +336,9 @@ SET_BOOL_PROP(EclBaseProblem, EnableThermalFluxBoundaries, false); SET_BOOL_PROP(EclBaseProblem, EnableTracerModel, false); +// timesteps are managed outside the Simulator class +SET_BOOL_PROP(EclBaseProblem, SimulatorManageTimeStep, false); + END_PROPERTIES namespace Ewoms { diff --git a/ewoms/common/simulator.hh b/ewoms/common/simulator.hh index e3cb773b5b..79a0221ccb 100644 --- a/ewoms/common/simulator.hh +++ b/ewoms/common/simulator.hh @@ -56,6 +56,7 @@ NEW_PROP_TAG(EndTime); NEW_PROP_TAG(RestartTime); NEW_PROP_TAG(InitialTimeStepSize); NEW_PROP_TAG(PredeterminedTimeStepsFile); +NEW_PROP_TAG(SimulatorManageTimeStep); END_PROPERTIES @@ -82,6 +83,8 @@ class Simulator typedef typename GET_PROP_TYPE(TypeTag, Model) Model; typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; + enum { simulatorManageTimeStep = GET_PROP_VALUE(TypeTag, SimulatorManageTimeStep) }; + public: // do not allow to copy simulators around Simulator(const Simulator& ) = delete; @@ -350,15 +353,14 @@ public: */ Scalar timeStepSize() const { - return timeStepSize_; - - /* This should be handeled by ghe time stepper - Scalar maximumTimeStepSize = - std::min(episodeMaxTimeStepSize(), - std::max( Scalar(0), endTime() - this->time())); - - return std::min(timeStepSize_, maximumTimeStepSize); - */ + if (simulatorManageTimeStep) { + Scalar maximumTimeStepSize = + std::min(episodeMaxTimeStepSize(), + std::max( Scalar(0), endTime() - this->time())); + return std::min(timeStepSize_, maximumTimeStepSize); + } else { + return timeStepSize_; + } } /*! diff --git a/ewoms/disc/common/fvbasediscretization.hh b/ewoms/disc/common/fvbasediscretization.hh index 082def77b0..d729a82926 100644 --- a/ewoms/disc/common/fvbasediscretization.hh +++ b/ewoms/disc/common/fvbasediscretization.hh @@ -279,6 +279,9 @@ SET_BOOL_PROP(FvBaseDiscretization, ExtensiveStorageTerm, false); // use volumetric residuals is default SET_BOOL_PROP(FvBaseDiscretization, UseVolumetricResidual, true); +// by default the Simulator class ensures timesteps do not exceed simulation time etc. +SET_BOOL_PROP(FvBaseDiscretization, SimulatorManageTimeStep, true); + END_PROPERTIES From b0b56703567ba5302b1104c3102ef2f075417a3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 18 Jan 2019 15:34:56 +0100 Subject: [PATCH 36/38] Output is written for the end time of steps. --- ebos/eclwriter.hh | 2 +- ewoms/io/restart.hh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index 0c0cb9a784..449bec7a29 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -149,7 +149,7 @@ public: */ void writeOutput(bool isSubStep) { - Scalar curTime = simulator_.time();// + simulator_.timeStepSize(); + Scalar curTime = simulator_.time() + simulator_.timeStepSize(); Scalar totalSolverTime = simulator_.executionTimer().realTimeElapsed(); Scalar nextStepSize = simulator_.problem().nextTimeStepSize(); diff --git a/ewoms/io/restart.hh b/ewoms/io/restart.hh index cd22b72332..478e85f97a 100644 --- a/ewoms/io/restart.hh +++ b/ewoms/io/restart.hh @@ -138,7 +138,7 @@ public: { const std::string magicCookie = magicRestartCookie_(simulator.gridView()); fileName_ = restartFileName_(simulator, - simulator.time()); + simulator.time() + simulator.timeStepSize()); // open output file and write magic cookie outStream_.open(fileName_.c_str()); From 00b288ce9bdfc253dd7dbbab62e93db3fa805e5e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 Jan 2019 12:36:22 +0100 Subject: [PATCH 37/38] Add defaulted atEndOfStep arguments to serialization. --- ewoms/common/simulator.hh | 4 ++-- ewoms/io/restart.hh | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/ewoms/common/simulator.hh b/ewoms/common/simulator.hh index 79a0221ccb..85f74092cf 100644 --- a/ewoms/common/simulator.hh +++ b/ewoms/common/simulator.hh @@ -798,11 +798,11 @@ public: * name and uses the extension .ers. (Ewoms ReStart * file.) See Ewoms::Restart for details. */ - void serialize() + void serialize(const bool atEndOfStep = true) { typedef Ewoms::Restart Restarter; Restarter res; - res.serializeBegin(*this); + res.serializeBegin(*this, atEndOfStep); if (gridView().comm().rank() == 0) std::cout << "Serialize to file '" << res.fileName() << "'" << ", next time step size: " << timeStepSize() diff --git a/ewoms/io/restart.hh b/ewoms/io/restart.hh index 478e85f97a..c98f17b718 100644 --- a/ewoms/io/restart.hh +++ b/ewoms/io/restart.hh @@ -134,11 +134,11 @@ public: * \brief Write the current state of the model to disk. */ template - void serializeBegin(Simulator& simulator) + void serializeBegin(Simulator& simulator, const bool atEndOfStep = true) { const std::string magicCookie = magicRestartCookie_(simulator.gridView()); - fileName_ = restartFileName_(simulator, - simulator.time() + simulator.timeStepSize()); + const double serializationTime = atEndOfStep ? simulator.time() + simulator.timeStepSize() : simulator.time(); + fileName_ = restartFileName_(simulator, serializationTime); // open output file and write magic cookie outStream_.open(fileName_.c_str()); From 6c46be5c9b3019d1c234acb7b1630c331926893d Mon Sep 17 00:00:00 2001 From: hnil Date: Fri, 31 May 2019 21:38:18 +0200 Subject: [PATCH 38/38] added comment on hack to get adjoint to work --- ewoms/disc/common/fvbaselocalresidual.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ewoms/disc/common/fvbaselocalresidual.hh b/ewoms/disc/common/fvbaselocalresidual.hh index 581630cf1d..607ebb1613 100644 --- a/ewoms/disc/common/fvbaselocalresidual.hh +++ b/ewoms/disc/common/fvbaselocalresidual.hh @@ -163,14 +163,14 @@ public: assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0)); residual = 0.0; - if(elemCtx.focusTimeIndex()==0){ + if(elemCtx.focusTimeIndex()==0){//NB: This is only valid if fluxeterm sis purely implicite, This was nesseary to not get false derivatives form wells in adjiont. // evaluate the flux terms asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0); } // evaluate the storage and the source terms asImp_().evalVolumeTerms_(residual, elemCtx); - if(elemCtx.focusTimeIndex()==0){ + if(elemCtx.focusTimeIndex()==0){//NB: This is only valid if fluxeterm sis purely implicite, This was nesseary to not get false derivatives form wells in adjiont. // evaluate the boundary conditions asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0); }