Merge pull request #1739 from andlaus/misc_fixes

Miscellaneous fixes and improvements

ebos/eclequilinitializer.hh

@@ -31,6 +31,7 @@
 #include "equil/initstateequil.hh"
 
 #include <ewoms/common/propertysystem.hh>
+#include <ewoms/models/blackoil/blackoilproperties.hh>
 
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
 #include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>

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)

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))
         {

ebos/eclwellmanager.hh

@@ -406,24 +406,20 @@ public:
             wells_[wellIdx]->beginIterationPreProcess();
 
         // call the accumulation routines
-        ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(simulator_.vanguard().gridView());
+        ElementContext elemCtx(simulator_);
-#ifdef _OPENMP
+        const auto gridView = simulator_.vanguard().gridView();
-#pragma omp parallel
+        auto elemIt = gridView.template begin</*codim=*/0>();
-#endif
+        const auto& elemEndIt = gridView.template end</*codim=*/0>();
-        {
+        for (; elemIt != elemEndIt; ++elemIt) {
-            ElementContext elemCtx(simulator_);
+            const Element& elem = *elemIt;
-            auto elemIt = threadedElemIt.beginParallel();
+            if (elem.partitionType() != Dune::InteriorEntity)
-            for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
+                continue;
-                const Element& elem = *elemIt;
-                if (elem.partitionType() != Dune::InteriorEntity)
-                    continue;
 
-                elemCtx.updatePrimaryStencil(elem);
+            elemCtx.updatePrimaryStencil(elem);
-                elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+            elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
 
-                for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx)
+            for (size_t wellIdx = 0; wellIdx < wellSize; ++wellIdx)
-                    wells_[wellIdx]->beginIterationAccumulate(elemCtx, /*timeIdx=*/0);
+                wells_[wellIdx]->beginIterationAccumulate(elemCtx, /*timeIdx=*/0);
-            }
         }
 
         // call the postprocessing routines

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);

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:

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 " <<

opm/autodiff/WellInterface.hpp

@@ -47,6 +47,9 @@
 #include <opm/simulators/timestepping/ConvergenceReport.hpp>
 #include <opm/simulators/DeferredLogger.hpp>
 
+#include <ewoms/models/blackoil/blackoilpolymermodules.hh>
+#include <ewoms/models/blackoil/blackoilsolventmodules.hh>
+
 #include<dune/common/fmatrix.hh>
 #include<dune/istl/bcrsmatrix.hh>
 #include<dune/istl/matrixmatrix.hh>