From 0a9d6a0760e028597c2c2e640dacee5f49df0319 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 1 Mar 2019 10:36:29 +0100 Subject: [PATCH 1/8] include missing header files this makes the well model and the equil initializer header more autonomous. --- ebos/eclequilinitializer.hh | 1 + opm/autodiff/WellInterface.hpp | 3 +++ 2 files changed, 4 insertions(+) diff --git a/ebos/eclequilinitializer.hh b/ebos/eclequilinitializer.hh index a30bb85c4..f1a70de13 100644 --- a/ebos/eclequilinitializer.hh +++ b/ebos/eclequilinitializer.hh @@ -31,6 +31,7 @@ #include "equil/initstateequil.hh" #include +#include #include #include diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index a64de1a5c..2ebb86ae4 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -47,6 +47,9 @@ #include #include +#include +#include + #include #include #include From 43dd9928b4d044fbe4477d652f0e0c98af20f47e Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 1 Mar 2019 10:36:29 +0100 Subject: [PATCH 2/8] EclNewtonMethod: Fix convergence criterion if the residual is volumetric --- ebos/eclnewtonmethod.hh | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/ebos/eclnewtonmethod.hh b/ebos/eclnewtonmethod.hh index fff1577c7..d4b721da4 100644 --- a/ebos/eclnewtonmethod.hh +++ b/ebos/eclnewtonmethod.hh @@ -159,16 +159,25 @@ public: } const auto& r = currentResidual[dofIdx]; - const double pvValue = + Scalar pvValue = this->simulator_.problem().porosity(dofIdx) * this->model().dofTotalVolume(dofIdx); sumPv += pvValue; bool cnvViolated = false; + Scalar dofVolume = this->model().dofTotalVolume(dofIdx); + for (unsigned eqIdx = 0; eqIdx < r.size(); ++eqIdx) { Scalar tmpError = r[eqIdx] * dt * this->model().eqWeight(dofIdx, eqIdx) / pvValue; Scalar tmpError2 = r[eqIdx] * this->model().eqWeight(dofIdx, eqIdx); + // in the case of a volumetric formulation, the residual in the above is + // per cubic meter + if (GET_PROP_VALUE(TypeTag, UseVolumetricResidual)) { + tmpError *= dofVolume; + tmpError2 *= dofVolume; + } + this->error_ = Opm::max(std::abs(tmpError), this->error_); if (std::abs(tmpError) > this->tolerance_) From 8e5f1279f3f58959aeb476a8c6777c3cf4ad25f5 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 1 Mar 2019 10:36:29 +0100 Subject: [PATCH 3/8] EclWellManager: do not use parallelism the speedup gained by parallelism here are simply not worth the headaches. note that `flow` is unaffected by this because it uses `Opm::BlackoilWellModel`. --- ebos/eclwellmanager.hh | 28 ++++++++++++---------------- 1 file changed, 12 insertions(+), 16 deletions(-) diff --git a/ebos/eclwellmanager.hh b/ebos/eclwellmanager.hh index 4adc841e5..7d850f2ac 100644 --- a/ebos/eclwellmanager.hh +++ b/ebos/eclwellmanager.hh @@ -406,24 +406,20 @@ public: wells_[wellIdx]->beginIterationPreProcess(); // call the accumulation routines - ThreadedEntityIterator threadedElemIt(simulator_.vanguard().gridView()); -#ifdef _OPENMP -#pragma omp parallel -#endif - { - ElementContext elemCtx(simulator_); - auto elemIt = threadedElemIt.beginParallel(); - for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) { - const Element& elem = *elemIt; - if (elem.partitionType() != Dune::InteriorEntity) - continue; + ElementContext elemCtx(simulator_); + const auto gridView = simulator_.vanguard().gridView(); + auto elemIt = gridView.template begin(); + const auto& elemEndIt = gridView.template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const Element& elem = *elemIt; + if (elem.partitionType() != Dune::InteriorEntity) + continue; - elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx) - wells_[wellIdx]->beginIterationAccumulate(elemCtx, /*timeIdx=*/0); - } + for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx) + wells_[wellIdx]->beginIterationAccumulate(elemCtx, /*timeIdx=*/0); } // call the postprocessing routines From 3ca8b8f2854e97b3e9f89dcd5c9dc61557cbc424 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 1 Mar 2019 10:36:29 +0100 Subject: [PATCH 4/8] ebos: improve the messages printed during the run the convergence behaviour can now be understood and the report step information is printed, too. This does not affect `flow`, becase it implements its own newton and time stepping routines. --- ebos/eclnewtonmethod.hh | 2 ++ ebos/eclproblem.hh | 12 ++++++++++++ 2 files changed, 14 insertions(+) diff --git a/ebos/eclnewtonmethod.hh b/ebos/eclnewtonmethod.hh index d4b721da4..5471387ed 100644 --- a/ebos/eclnewtonmethod.hh +++ b/ebos/eclnewtonmethod.hh @@ -210,6 +210,8 @@ public: Scalar y = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonSumToleranceExponent); sumTolerance_ = x*std::pow(sumPv, y); + this->endIterMsg() << " (max: " << this->tolerance_ << ", violated for " << errorPvFraction_*100 << "% of the pore volume), aggegate error: " << errorSum_ << " (max: " << sumTolerance_ << ")"; + // make sure that the error never grows beyond the maximum // allowed one if (this->error_ > newtonMaxError) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 72365a2f9..ce8ff3af3 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -707,6 +707,18 @@ public: // The first thing to do in the morning of an episode is update update the // eclState and the deck if they need to be changed. int nextEpisodeIdx = simulator.episodeIndex(); + + if (this->gridView().comm().rank() == 0) { + boost::posix_time::ptime curDateTime = + boost::posix_time::from_time_t(timeMap.getStartTime(nextEpisodeIdx+1)); + std::cout << "Report step " << nextEpisodeIdx + 2 + << "/" << timeMap.numTimesteps() + << " at day " << timeMap.getTimePassedUntil(nextEpisodeIdx+1)/(24*3600) + << "/" << timeMap.getTotalTime()/(24*3600) + << ", date = " << curDateTime.date() + << "\n "; + } + if (nextEpisodeIdx > 0 && events.hasEvent(Opm::ScheduleEvents::GEO_MODIFIER, nextEpisodeIdx)) { From b60305082db45525864bd273deda4a589dc17c8e Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 1 Mar 2019 10:36:29 +0100 Subject: [PATCH 5/8] EclBaseProblem: fix the parent type tags for the EclBaseProblem tag in the non-default cases this bitrot a bit because it was never seen by the compiler. (I still did not check if `ebos` compiles and works if `CpGrid` is replaced by dune-alugrid or `PolyhedralGrid`.) --- ebos/eclproblem.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index ce8ff3af3..be896e4ad 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -111,10 +111,10 @@ class EclProblem; BEGIN_PROPERTIES #if EBOS_USE_ALUGRID -NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclAluGridVanguard, EclOutputBlackOil)); +NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclAluGridVanguard, EclOutputBlackOil, VtkEclTracer)); #else NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclCpGridVanguard, EclOutputBlackOil, VtkEclTracer)); -//NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclPolyhedralGridVanguard, EclOutputBlackOil)); +//NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclPolyhedralGridVanguard, EclOutputBlackOil, VtkEclTracer)); #endif // The class which deals with ECL wells From 21f5c1fbd3a63765b24975d34eecf015c7f88521 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 1 Mar 2019 10:36:29 +0100 Subject: [PATCH 6/8] ISTLSolverEbos: throw `NumericalIssue` instead of `LinearSolverProblem` the former is caught by `ebos`, while the latter isn't. Alternatively, this can be fixed by deriving `LinearSolverProblem` from `NumericalIssue`, if preferred. --- opm/autodiff/ISTLSolverEbos.hpp | 2 +- opm/autodiff/ParallelOverlappingILU0.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index 52ad5674d..e00fd6a16 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -534,7 +534,7 @@ protected: // Check for failure of linear solver. if (!parameters_.ignoreConvergenceFailure_ && !result.converged) { const std::string msg("Convergence failure for linear solver."); - OPM_THROW_NOLOG(LinearSolverProblem, msg); + OPM_THROW_NOLOG(NumericalIssue, msg); } } protected: diff --git a/opm/autodiff/ParallelOverlappingILU0.hpp b/opm/autodiff/ParallelOverlappingILU0.hpp index caf0b41be..2698841ef 100644 --- a/opm/autodiff/ParallelOverlappingILU0.hpp +++ b/opm/autodiff/ParallelOverlappingILU0.hpp @@ -935,7 +935,7 @@ protected: milun_decomposition( A, iluIteration, milu, *ILU, *reorderer, *inverseReorderer ); } } - catch ( const Dune::MatrixBlockError& error ) + catch (const Dune::MatrixBlockError& error) { message = error.what(); std::cerr<<"Exception occured on process " << rank << " during " << From c37bc1cf38e9fffc828c3c73134af11536bff95f Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 1 Mar 2019 10:36:29 +0100 Subject: [PATCH 7/8] flow: hide a few unused eWoms parameters these parameters where introduced with support for the TUNIING keyword in `ebos`. since `flow` implements its own time stepping these parameters are unused and should thus be hidden from view in it. --- opm/autodiff/FlowMainEbos.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index eb4d20bad..c7b7eb7bc 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -151,6 +151,11 @@ namespace Opm EWOMS_HIDE_PARAM(TypeTag, NewtonTargetIterations); EWOMS_HIDE_PARAM(TypeTag, NewtonVerbose); EWOMS_HIDE_PARAM(TypeTag, NewtonWriteConvergence); + EWOMS_HIDE_PARAM(TypeTag, EclNewtonSumTolerance); + EWOMS_HIDE_PARAM(TypeTag, EclNewtonSumToleranceExponent); + EWOMS_HIDE_PARAM(TypeTag, EclNewtonStrictIterations); + EWOMS_HIDE_PARAM(TypeTag, EclNewtonRelaxedVolumeFraction); + EWOMS_HIDE_PARAM(TypeTag, EclNewtonRelaxedTolerance); // the default eWoms checkpoint/restart mechanism does not work with flow EWOMS_HIDE_PARAM(TypeTag, RestartTime); From d7e74e7c4e450ffe979bfcd9951078dbc201210b Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Mon, 4 Mar 2019 13:59:50 +0100 Subject: [PATCH 8/8] fix review comment by [at]totto82 --- ebos/eclproblem.hh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index be896e4ad..dd237c19b 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -370,6 +370,7 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; enum { numPhases = FluidSystem::numPhases }; enum { numComponents = FluidSystem::numComponents }; + enum { enableExperiments = GET_PROP_VALUE(TypeTag, EnableExperiments) }; enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) }; enum { enablePolymerMolarWeight = GET_PROP_VALUE(TypeTag, EnablePolymerMW) }; @@ -708,7 +709,7 @@ public: // eclState and the deck if they need to be changed. int nextEpisodeIdx = simulator.episodeIndex(); - if (this->gridView().comm().rank() == 0) { + if (enableExperiments && this->gridView().comm().rank() == 0) { boost::posix_time::ptime curDateTime = boost::posix_time::from_time_t(timeMap.getStartTime(nextEpisodeIdx+1)); std::cout << "Report step " << nextEpisodeIdx + 2