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/ebos/eclnewtonmethod.hh b/ebos/eclnewtonmethod.hh index fff1577c7..5471387ed 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_) @@ -201,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..dd237c19b 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 @@ -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) }; @@ -707,6 +708,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 (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 + << "/" << 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)) { 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 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); 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 " << 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