diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 961294063..b3b34f6b3 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/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 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 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 - // are invalid - ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0); - - // solver has succeed i.e. no need for restart. - isRestart_ = false; + // since the solution was changed, the cache for the intensive quantities are invalid + // ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0); } + + 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 << " " < - void applyWellModel(const X& x, Y& y ) + void applyWellModelAdd(const X& x, Y& y ) { - wellModel().apply(x, y); + wellModel().apply(x, y); + } + + template + 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 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 SequentialOperator; - SequentialOperator& sOpA = static_cast< SequentialOperator& > (opA); - istlSolver().solve( sOpA, x, ebosResid ); + typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, false > Operator; + Operator opA(ebosJac, const_cast< ThisType& > (*this) ); + istlSolver().solve( opA, x, ebosResid ); } // recover wells. @@ -464,7 +477,7 @@ namespace Opm { Adapts a matrix to the assembled linear operator interface */ - template + template class WellModelMatrixAdapter : public Dune::AssembledLinearOperator { typedef Dune::AssembledLinearOperator BaseType; @@ -478,16 +491,18 @@ namespace Opm { #if HAVE_MPI typedef Dune::OwnerOverlapCopyCommunication communication_type; #else - typedef Dune::CollectiveCommunication 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 OverlappingWellModelMatrixAdapter : public WellModelMatrixAdapter - { - 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 diff --git a/opm/autodiff/ISTLSolver.hpp b/opm/autodiff/ISTLSolver.hpp index b7ca977b1..a5132c4a1 100644 --- a/opm/autodiff/ISTLSolver.hpp +++ b/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 - void constructPreconditionerAndSolve(O& opA, + template + 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 + template 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 + 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 ); } diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index bed4516bc..ea5161d90 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/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 {