Merge pull request #932 from dr-robertk/PR/bug-hunt2

Bug fix for reusing cached intensive quantities.

opm/autodiff/BlackoilModelEbos.hpp

@@ -184,7 +184,7 @@ namespace Opm {
         , current_relaxation_(1.0)
         , dx_old_(AutoDiffGrid::numCells(grid_))
         , isBeginReportStep_(false)
-        , isRestart_(false)
+        , invalidateIntensiveQuantitiesCache_(true)
         {
             const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
             const std::vector<double> pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
@@ -258,16 +258,21 @@ namespace Opm {
                 current_relaxation_ = 1.0;
                 dx_old_ = 0.0;
             }
+
+            // reset intensive quantities cache useless other options are set
+            // further down
+            invalidateIntensiveQuantitiesCache_ = true;
+
             IterationReport iter_report = assemble(timer, iteration, reservoir_state, well_state);
             std::vector<double> residual_norms;
             const bool converged = getConvergence(timer, iteration,residual_norms);
             residual_norms_history_.push_back(residual_norms);
             bool must_solve = (iteration < nonlinear_solver.minIter()) || (!converged);
-            // is first set to true if a linear solve is needed, but then it is set to false if the solver succeed.
+
-            isRestart_ = must_solve && (iteration == nonlinear_solver.maxIter());
             // don't solve if we have reached the maximum number of iteration.
             must_solve = must_solve && (iteration < nonlinear_solver.maxIter());
             if (must_solve) {
+
                 // enable single precision for solvers when dt is smaller then 20 days
                 //residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ;
 
@@ -276,6 +281,7 @@ namespace Opm {
                 const int nw = wellModel().wells().number_of_wells;
                 BVector x(nc);
                 BVector xw(nw);
+
                 solveJacobianSystem(x, xw);
 
                 // Stabilize the nonlinear update.
@@ -298,17 +304,21 @@ namespace Opm {
                 updateState(x,reservoir_state);
                 wellModel().updateWellState(xw, well_state);
 
-                // since the solution was changed, the cache for the intensive quantities
+                // since the solution was changed, the cache for the intensive quantities are invalid
-                // are invalid
+                // ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
-                ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
-
-                // solver has succeed i.e. no need for restart.
-                isRestart_ = false;
             }
+
+            if( converged && (iteration >= nonlinear_solver.minIter()) )
+            {
+                // in case of convergence we do not need to reset intensive quantities
+                invalidateIntensiveQuantitiesCache_ = false ;
+            }
+
             const bool failed = false; // Not needed in this model.
             const int linear_iters = must_solve ? linearIterationsLastSolve() : 0;
             return IterationReport{ failed, converged, linear_iters, iter_report.well_iterations };
         }
+
         void printIf(int c, double x, double y, double eps, std::string type) {
             if (std::abs(x-y) > eps) {
                 std::cout << type << " " <<c << ": "<<x << " " << y << std::endl;
@@ -351,7 +361,6 @@ namespace Opm {
             }
             catch ( const Dune::FMatrixError& e  )
             {
-                isRestart_ = true;
                 OPM_THROW(Opm::NumericalProblem,"Well equation did not converge");
             }
 
@@ -414,9 +423,15 @@ namespace Opm {
         }
 
         template <class X, class Y>
-        void applyWellModel(const X& x, Y& y )
+        void applyWellModelAdd(const X& x, Y& y )
         {
-           wellModel().apply(x, y);
+            wellModel().apply(x, y);
+        }
+
+        template <class X, class Y>
+        void applyWellModelScaleAdd(const Scalar alpha, const X& x, Y& y )
+        {
+            wellModel().applyScaleAdd(alpha, x, y);
         }
 
 
@@ -427,27 +442,25 @@ namespace Opm {
             const auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
             auto& ebosResid = ebosSimulator_.model().linearizer().residual();
 
-            typedef OverlappingWellModelMatrixAdapter<Mat,BVector,BVector, ThisType> Operator;
-            Operator opA(ebosJac, const_cast< ThisType& > (*this), istlSolver().parallelInformation() );
-
             // apply well residual to the residual.
             wellModel().apply(ebosResid);
 
             // set initial guess
             x = 0.0;
 
-            typedef typename Operator :: communication_type Comm;
-            Comm* comm = opA.comm();
             // Solve system.
-            if( comm )
+            if( isParallel() )
             {
-                istlSolver().solve( opA, x, ebosResid, *comm );
+                typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, true > Operator;
+                Operator opA(ebosJac, const_cast< ThisType& > (*this), istlSolver().parallelInformation() );
+                assert( opA.comm() );
+                istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
             }
             else
             {
-                typedef WellModelMatrixAdapter<Mat,BVector,BVector, ThisType> SequentialOperator;
+                typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, false > Operator;
-                SequentialOperator& sOpA = static_cast< SequentialOperator& > (opA);
+                Operator opA(ebosJac, const_cast< ThisType& > (*this) );
-                istlSolver().solve( sOpA, x, ebosResid );
+                istlSolver().solve( opA, x, ebosResid );
             }
 
             // recover wells.
@@ -464,7 +477,7 @@ namespace Opm {
 
            Adapts a matrix to the assembled linear operator interface
          */
-        template<class M, class X, class Y, class WellModel>
+        template<class M, class X, class Y, class WellModel, bool overlapping >
         class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
         {
           typedef Dune::AssembledLinearOperator<M,X,Y> BaseType;
@@ -478,16 +491,18 @@ namespace Opm {
 #if HAVE_MPI
           typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
 #else
-          typedef Dune::CollectiveCommunication<int> communication_type;
+          typedef Dune::CollectiveCommunication< Grid > communication_type;
 #endif
 
           enum {
             //! \brief The solver category.
-            category=Dune::SolverCategory::sequential
+            category = overlapping ?
+                Dune::SolverCategory::overlapping :
+                Dune::SolverCategory::sequential
           };
 
           //! constructor: just store a reference to a matrix
-          WellModelMatrixAdapter (const M& A, WellModel& wellMod, const boost::any& parallelInformation )
+          WellModelMatrixAdapter (const M& A, WellModel& wellMod, const boost::any& parallelInformation = boost::any() )
               : A_( A ), wellMod_( wellMod ), comm_()
           {
 #if HAVE_MPI
@@ -503,7 +518,8 @@ namespace Opm {
           virtual void apply( const X& x, Y& y ) const
           {
             A_.mv( x, y );
-            wellMod_.applyWellModel(x, y );
+            // add well model modification to y
+            wellMod_.applyWellModelAdd(x, y );
 
 #if HAVE_MPI
             if( comm_ )
@@ -511,10 +527,12 @@ namespace Opm {
 #endif
           }
 
+          // y += \alpha * A * x
           virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
           {
             A_.usmv(alpha,x,y);
-            wellMod_.applyWellModel(x, y );
+            // add scaled well model modification to y
+            wellMod_.applyWellModelScaleAdd( alpha, x, y );
 
 #if HAVE_MPI
             if( comm_ )
@@ -535,24 +553,6 @@ namespace Opm {
           std::unique_ptr< communication_type > comm_;
         };
 
-        template<class M, class X, class Y, class WellModel>
-        class OverlappingWellModelMatrixAdapter : public WellModelMatrixAdapter<M,X,Y,WellModel>
-        {
-        public:
-          typedef WellModelMatrixAdapter< M,X,Y,WellModel > BaseType;
-
-          enum {
-            //! \brief The solver category.
-            category=Dune::SolverCategory::overlapping
-          };
-
-          //! constructor: just store a reference to a matrix
-          OverlappingWellModelMatrixAdapter(const M& A, WellModel& wellMod, const boost::any& parallelInformation )
-              : BaseType( A, wellMod, parallelInformation )
-          {}
-        };
-
-
         /// Apply an update to the primary variables, chopped if appropriate.
         /// \param[in]      dx                updates to apply to primary variables
         /// \param[in, out] reservoir_state   reservoir state variables
@@ -950,13 +950,11 @@ namespace Opm {
                 if (std::isnan(mass_balance_residual[phaseIdx])
                     || std::isnan(CNV[phaseIdx])
                     || (phaseIdx < np && std::isnan(well_flux_residual[phaseIdx]))) {
-                    isRestart_ = true;
                     OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName);
                 }
                 if (mass_balance_residual[phaseIdx] > maxResidualAllowed()
                     || CNV[phaseIdx] > maxResidualAllowed()
                     || (phaseIdx < np && well_flux_residual[phaseIdx] > maxResidualAllowed())) {
-                    isRestart_ = true;
                     OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName);
                 }
             }
@@ -1425,7 +1423,7 @@ namespace Opm {
                 ebosSimulator_.problem().beginTimeStep();
             }
             // if the last step failed we want to recalculate the IntesiveQuantities.
-            if (isRestart_) {
+            if ( invalidateIntensiveQuantitiesCache_ ) {
                 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
             }
 
@@ -1454,8 +1452,7 @@ namespace Opm {
 
     public:
         bool isBeginReportStep_;
-        bool isRestart_;
+        bool invalidateIntensiveQuantitiesCache_;
-
     };
 } // namespace Opm
 

opm/autodiff/ISTLSolver.hpp

@@ -191,8 +191,8 @@ namespace Opm
         /// \brief construct the CPR preconditioner and the solver.
         /// \tparam P The type of the parallel information.
         /// \param parallelInformation the information about the parallelization.
-        template<int category=Dune::SolverCategory::sequential, class O, class POrComm>
+        template<int category=Dune::SolverCategory::sequential, class LinearOperator, class POrComm>
-        void constructPreconditionerAndSolve(O& opA,
+        void constructPreconditionerAndSolve(LinearOperator& linearOperator,
                                              Vector& x, Vector& istlb,
                                              const POrComm& parallelInformation_arg,
                                              Dune::InverseOperatorResult& result) const
@@ -211,21 +211,33 @@ namespace Opm
             {
                 typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm>  CPRSelectorType;
                 typedef typename CPRSelectorType::AMG AMG;
+                typedef typename CPRSelectorType::Operator MatrixOperator;
+
                 std::unique_ptr< AMG > amg;
+                std::unique_ptr< MatrixOperator > opA;
+
+                if( ! std::is_same< LinearOperator, MatrixOperator > :: value )
+                {
+                    // create new operator in case linear operator and matrix operator differ
+                    opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) );
+                }
+
+                const double relax = 1.0;
+
                 // Construct preconditioner.
-                constructAMGPrecond(opA, parallelInformation_arg, amg);
+                constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
 
                 // Solve.
-                solve(opA, x, istlb, *sp, *amg, result);
+                solve(linearOperator, x, istlb, *sp, *amg, result);
             }
             else
 #endif
             {
                 // Construct preconditioner.
-                auto precond = constructPrecond(opA, parallelInformation_arg);
+                auto precond = constructPrecond(linearOperator, parallelInformation_arg);
 
                 // Solve.
-                solve(opA, x, istlb, *sp, *precond, result);
+                solve(linearOperator, x, istlb, *sp, *precond, result);
             }
         }
 
@@ -252,11 +264,18 @@ namespace Opm
         }
 #endif
 
-        template <class Operator, class POrComm, class AMG >
+        template <class LinearOperator, class MatrixOperator, class POrComm, class AMG >
         void
-        constructAMGPrecond(Operator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg ) const
+        constructAMGPrecond(LinearOperator& linearOperator, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const
+        {
+            ISTLUtility::createAMGPreconditionerPointer( *opA, relax, comm, amg );
+        }
+
+
+        template <class MatrixOperator, class POrComm, class AMG >
+        void
+        constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const
         {
-            const double relax = 1.0;
             ISTLUtility::createAMGPreconditionerPointer( opA, relax, comm, amg );
         }
 

opm/autodiff/StandardWellsDense.hpp

@@ -336,6 +336,22 @@ namespace Opm {
                 duneB_.mmtv(invDCx,Ax);
             }
 
+            // apply well model with scaling of alpha
+            void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax)
+            {
+                if ( ! localWellsActive() ) {
+                    return;
+                }
+
+                if( scaleAddRes_.size() != Ax.size() ) {
+                    scaleAddRes_.resize( Ax.size() );
+                }
+
+                scaleAddRes_ = 0.0;
+                apply( x, scaleAddRes_ );
+                Ax.axpy( alpha, scaleAddRes_ );
+            }
+
             // xw = inv(D)*(rw - C*x)
             void recoverVariable(const BVector& x, BVector& xw) const {
                 if ( ! localWellsActive() ) {
@@ -1418,6 +1434,7 @@ namespace Opm {
 
             mutable BVector Cx_;
             mutable BVector invDrw_;
+            mutable BVector scaleAddRes_;
 
             // protected methods
             EvalWell getBhp(const int wellIdx) const {