From 439a084508946fcb1dfd87d33668f07d3965b107 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Thu, 17 Nov 2016 11:55:30 +0100 Subject: [PATCH 1/4] [bugfix][BlockoilModelEbos] fix invalidation of intensive quantities after linear solver failure. --- opm/autodiff/BlackoilModelEbos.hpp | 39 ++++++++++++++++++------------ 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 961294063..9b2fb160b 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 ) + { + // 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 << " " < 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 +1433,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 +1462,7 @@ namespace Opm { public: bool isBeginReportStep_; - bool isRestart_; - + bool invalidateIntensiveQuantitiesCache_; }; } // namespace Opm From a0da20378c8c8292fd91548f0d8f51db1711b95b Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Thu, 17 Nov 2016 13:43:39 +0100 Subject: [PATCH 2/4] [bugfix][WellModelMatrixAdapter] fix applyscaleadd method. --- opm/autodiff/BlackoilModelEbos.hpp | 20 ++++++++++++++------ opm/autodiff/StandardWellsDense.hpp | 17 +++++++++++++++++ 2 files changed, 31 insertions(+), 6 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 9b2fb160b..7d88503ee 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -308,7 +308,7 @@ namespace Opm { // ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0); } - if( converged ) + if( converged && (iteration >= nonlinear_solver.minIter()) ) { // in case of convergence we do not need to reset intensive quantities invalidateIntensiveQuantitiesCache_ = false ; @@ -423,9 +423,15 @@ namespace Opm { } template - 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); } @@ -487,7 +493,7 @@ namespace Opm { #if HAVE_MPI typedef Dune::OwnerOverlapCopyCommunication communication_type; #else - typedef Dune::CollectiveCommunication communication_type; + typedef Dune::CollectiveCommunication< Grid > communication_type; #endif enum { @@ -512,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_ ) @@ -524,7 +531,8 @@ namespace Opm { 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_ ) 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 { From 59f40ba14e853115fc9a0de6badc42b49af92356 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Thu, 17 Nov 2016 14:24:16 +0100 Subject: [PATCH 3/4] [cleanup][WellModelMatrixAdapter] use only one implementation of the matrix adapter to avoid confusion. --- opm/autodiff/BlackoilModelEbos.hpp | 44 +++++++++--------------------- 1 file changed, 13 insertions(+), 31 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 7d88503ee..b3b34f6b3 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -442,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. @@ -479,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; @@ -498,11 +496,13 @@ namespace Opm { 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 @@ -553,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 From 5fce54fd447a6e9412c5c9aaebcb5873f6511597 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Thu, 17 Nov 2016 16:17:41 +0100 Subject: [PATCH 4/4] [bugfix][ISTLSolver] make code compile with AMG when different matrix operator is used. --- opm/autodiff/ISTLSolver.hpp | 37 ++++++++++++++++++++++++++++--------- 1 file changed, 28 insertions(+), 9 deletions(-) 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 ); }