From 2d2d0c44336d151df614557fa2e973099d3e6ead Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 12 Nov 2018 11:03:54 +0100 Subject: [PATCH] Adapt to linear solver api from ewoms --- opm/autodiff/BlackoilDetails.hpp | 1 + opm/autodiff/BlackoilModelEbos.hpp | 200 +++-------------- opm/autodiff/BlackoilModelParametersEbos.hpp | 7 - opm/autodiff/BlackoilWellModel.hpp | 2 +- opm/autodiff/FlowLinearSolverParameters.hpp | 13 ++ opm/autodiff/FlowMainEbos.hpp | 43 +--- opm/autodiff/ISTLSolverEbos.hpp | 207 +++++++++++++++++- .../SimulatorFullyImplicitBlackoilEbos.hpp | 7 +- 8 files changed, 241 insertions(+), 239 deletions(-) diff --git a/opm/autodiff/BlackoilDetails.hpp b/opm/autodiff/BlackoilDetails.hpp index e809acf4b..d51793446 100644 --- a/opm/autodiff/BlackoilDetails.hpp +++ b/opm/autodiff/BlackoilDetails.hpp @@ -25,6 +25,7 @@ #define OPM_BLACKOILDETAILS_HEADER_INCLUDED #include +#include namespace Opm { namespace detail { diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 24bfec244..bef385486 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -86,6 +86,8 @@ SET_BOOL_PROP(EclFlowProblem, EnableTemperature, true); SET_BOOL_PROP(EclFlowProblem, EnableEnergy, false); SET_TYPE_PROP(EclFlowProblem, EclWellModel, Opm::BlackoilWellModel); +SET_TAG_PROP(EclFlowProblem, LinearSolverSplice, FlowIstlSolver); + END_PROPERTIES @@ -149,11 +151,9 @@ namespace Opm { BlackoilModelEbos(Simulator& ebosSimulator, const ModelParameters& param, BlackoilWellModel& well_model, - const ISTLSolverType& linsolver, const bool terminal_output) : ebosSimulator_(ebosSimulator) , grid_(ebosSimulator_.vanguard().grid()) - , istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) ) , phaseUsage_(phaseUsageFromDeck(eclState())) , has_disgas_(FluidSystem::enableDissolvedGas()) , has_vapoil_(FluidSystem::enableVaporizedOil()) @@ -169,12 +169,6 @@ namespace Opm { { // compute global sum of number of cells global_nc_ = detail::countGlobalCells(grid_); - //find rows of matrix corresponding to overlap - detail::findOverlapRowsAndColumns(grid_,overlapRowAndColumns_); - if (!istlSolver_) - { - OPM_THROW(std::logic_error,"solver down cast to ISTLSolver failed"); - } convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves. } @@ -375,16 +369,6 @@ namespace Opm { ebosSimulator_.model().linearizer().linearize(); ebosSimulator_.problem().endIteration(); - auto& ebosJac = ebosSimulator_.model().linearizer().jacobian(); - if (param_.matrix_add_well_contributions_) { - wellModel().addWellContributions(ebosJac.istlMatrix()); - } - if ( param_.preconditioner_add_well_contributions_ && - ! param_.matrix_add_well_contributions_ ) { - matrix_for_preconditioner_ .reset(new Mat(ebosJac.istlMatrix())); - wellModel().addWellContributions(*matrix_for_preconditioner_); - } - return wellModel().lastReport(); } @@ -469,179 +453,39 @@ namespace Opm { /// Number of linear iterations used in last call to solveJacobianSystem(). int linearIterationsLastSolve() const { - return istlSolver().iterations(); - } - - /// Zero out off-diagonal blocks on rows corresponding to overlap cells - /// Diagonal blocks on ovelap rows are set to diag(1e100). - void makeOverlapRowsInvalid(Mat& ebosJacIgnoreOverlap) const - { - //value to set on diagonal - MatrixBlockType diag_block(0.0); - for (int eq = 0; eq < numEq; ++eq) - diag_block[eq][eq] = 1.0e100; - - //loop over precalculated overlap rows and columns - for (auto row = overlapRowAndColumns_.begin(); row != overlapRowAndColumns_.end(); row++ ) - { - int lcell = row->first; - //diagonal block set to large value diagonal - ebosJacIgnoreOverlap[lcell][lcell] = diag_block; - - //loop over off diagonal blocks in overlap row - for (auto col = row->second.begin(); col != row->second.end(); ++col) - { - int ncell = *col; - //zero out block - ebosJacIgnoreOverlap[lcell][ncell] = 0.0; - } - } + return ebosSimulator_.model().newtonMethod().linearSolver().iterations (); } /// Solve the Jacobian system Jx = r where J is the Jacobian and /// r is the residual. - void solveJacobianSystem(BVector& x) const + void solveJacobianSystem(BVector& x) { + + auto& ebosJac = ebosSimulator_.model().linearizer().jacobian(); + auto& ebosResid = ebosSimulator_.model().linearizer().residual(); // J = [A, B; C, D], where A is the reservoir equations, B and C the interaction of well // with the reservoir and D is the wells itself. // The full system is reduced to a number of cells X number of cells system via Schur complement // A -= B^T D^-1 C - // Instead of modifying A, the Ax operator is modified. i.e Ax -= B^T D^-1 C x in the WellModelMatrixAdapter. + // If matrix_add_well_contribution is false, the Ax operator is modified. i.e Ax -= B^T D^-1 C x in the WellModelMatrixAdapter + // instead of A. // The residual is modified similarly. // r = [r, r_well], where r is the residual and r_well the well residual. // r -= B^T * D^-1 r_well - - auto& ebosJac = ebosSimulator_.model().linearizer().jacobian(); - auto& ebosResid = ebosSimulator_.model().linearizer().residual(); - wellModel().apply(ebosResid); + if (param_.matrix_add_well_contributions_) { + wellModel().addWellContributions(ebosJac.istlMatrix()); + } // set initial guess x = 0.0; - const Mat& actual_mat_for_prec = matrix_for_preconditioner_ ? *matrix_for_preconditioner_.get() : ebosJac.istlMatrix(); - // Solve system. - if( isParallel() ) - { - typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel, true > Operator; + auto& ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver(); + ebosSolver.prepare(ebosJac, ebosResid); + ebosSolver.solve(x); + } - auto ebosJacIgnoreOverlap = Mat(ebosJac.istlMatrix()); - //remove ghost rows in local matrix - makeOverlapRowsInvalid(ebosJacIgnoreOverlap); - //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables - //to be certain that correct matrix is used for preconditioning. - Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel(), - istlSolver().parallelInformation() ); - assert( opA.comm() ); - istlSolver().solve( opA, x, ebosResid, *(opA.comm()) ); - } - else - { - typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel, false > Operator; - Operator opA(ebosJac.istlMatrix(), actual_mat_for_prec, wellModel()); - istlSolver().solve( opA, x, ebosResid ); - } - } - - //===================================================================== - // Implementation for ISTL-matrix based operator - //===================================================================== - - /*! - \brief Adapter to turn a matrix into a linear operator. - - Adapts a matrix to the assembled linear operator interface - */ - template - class WellModelMatrixAdapter : public Dune::AssembledLinearOperator - { - typedef Dune::AssembledLinearOperator BaseType; - - public: - typedef M matrix_type; - typedef X domain_type; - typedef Y range_type; - typedef typename X::field_type field_type; - -#if HAVE_MPI - typedef Dune::OwnerOverlapCopyCommunication communication_type; -#else - typedef Dune::CollectiveCommunication< Grid > communication_type; -#endif - -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - Dune::SolverCategory::Category category() const override - { - return overlapping ? - Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential; - } -#else - enum { - //! \brief The solver category. - category = overlapping ? - Dune::SolverCategory::overlapping : - Dune::SolverCategory::sequential - }; -#endif - - //! constructor: just store a reference to a matrix - WellModelMatrixAdapter (const M& A, - const M& A_for_precond, - const WellModel& wellMod, - const boost::any& parallelInformation = boost::any() ) - : A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), comm_() - { -#if HAVE_MPI - if( parallelInformation.type() == typeid(ParallelISTLInformation) ) - { - const ParallelISTLInformation& info = - boost::any_cast( parallelInformation); - comm_.reset( new communication_type( info.communicator() ) ); - } -#endif - } - - virtual void apply( const X& x, Y& y ) const - { - A_.mv( x, y ); - - // add well model modification to y - wellMod_.apply(x, y ); - -#if HAVE_MPI - if( comm_ ) - comm_->project( y ); -#endif - } - - // y += \alpha * A * x - virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const - { - A_.usmv(alpha,x,y); - - // add scaled well model modification to y - wellMod_.applyScaleAdd( alpha, x, y ); - -#if HAVE_MPI - if( comm_ ) - comm_->project( y ); -#endif - } - - virtual const matrix_type& getmat() const { return A_for_precond_; } - - communication_type* comm() - { - return comm_.operator->(); - } - - protected: - const matrix_type& A_ ; - const matrix_type& A_for_precond_ ; - const WellModel& wellMod_; - std::unique_ptr< communication_type > comm_; - }; /// Apply an update to the primary variables. void updateSolution(const BVector& dx) @@ -945,6 +789,15 @@ namespace Opm { std::vector B_avg(numEq, 0.0); auto report = getReservoirConvergence(timer.currentStepLength(), iteration, B_avg, residual_norms); report += wellModel().getWellConvergence(B_avg); + + // Throw if any NaN or too large residual found. + ConvergenceReport::Severity severity = report.severityOfWorstFailure(); + if (severity == ConvergenceReport::Severity::NotANumber) { + OPM_THROW(Opm::NumericalIssue, "NaN residual found!"); + } else if (severity == ConvergenceReport::Severity::TooLarge) { + OPM_THROW(Opm::NumericalIssue, "Too large residual found!"); + } + return report; } @@ -1030,9 +883,6 @@ namespace Opm { double current_relaxation_; BVector dx_old_; - std::unique_ptr matrix_for_preconditioner_; - std::vector>> overlapRowAndColumns_; - std::vector convergence_reports_; public: /// return the StandardWells object diff --git a/opm/autodiff/BlackoilModelParametersEbos.hpp b/opm/autodiff/BlackoilModelParametersEbos.hpp index 2181db853..ca12f0732 100644 --- a/opm/autodiff/BlackoilModelParametersEbos.hpp +++ b/opm/autodiff/BlackoilModelParametersEbos.hpp @@ -48,7 +48,6 @@ NEW_PROP_TAG(SolveWelleqInitially); NEW_PROP_TAG(UpdateEquationsScaling); NEW_PROP_TAG(UseUpdateStabilization); NEW_PROP_TAG(MatrixAddWellContributions); -NEW_PROP_TAG(PreconditionerAddWellContributions); // parameters for multisegment wells NEW_PROP_TAG(TolerancePressureMsWells); @@ -72,7 +71,6 @@ SET_BOOL_PROP(FlowModelParameters, SolveWelleqInitially, true); SET_BOOL_PROP(FlowModelParameters, UpdateEquationsScaling, false); SET_BOOL_PROP(FlowModelParameters, UseUpdateStabilization, true); SET_BOOL_PROP(FlowModelParameters, MatrixAddWellContributions, false); -SET_BOOL_PROP(FlowModelParameters, PreconditionerAddWellContributions, false); SET_SCALAR_PROP(FlowModelParameters, TolerancePressureMsWells, 0.01 *1e5); SET_SCALAR_PROP(FlowModelParameters, MaxPressureChangeMsWells, 2.0 *1e5); SET_BOOL_PROP(FlowModelParameters, UseInnerIterationsMsWells, true); @@ -155,9 +153,6 @@ namespace Opm // Whether to add influences of wells between cells to the matrix and preconditioner matrix bool matrix_add_well_contributions_; - // Whether to add influences of wells between cells to the preconditioner matrix only - bool preconditioner_add_well_contributions_; - /// Construct from user parameters or defaults. BlackoilModelParametersEbos() { @@ -181,7 +176,6 @@ namespace Opm update_equations_scaling_ = EWOMS_GET_PARAM(TypeTag, bool, UpdateEquationsScaling); use_update_stabilization_ = EWOMS_GET_PARAM(TypeTag, bool, UseUpdateStabilization); matrix_add_well_contributions_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); - preconditioner_add_well_contributions_ = EWOMS_GET_PARAM(TypeTag, bool, PreconditionerAddWellContributions); deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName); } @@ -208,7 +202,6 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, bool, UpdateEquationsScaling, "Update scaling factors for mass balance equations during the run"); EWOMS_REGISTER_PARAM(TypeTag, bool, UseUpdateStabilization, "Try to detect and correct oscillations or stagnation during the Newton method"); EWOMS_REGISTER_PARAM(TypeTag, bool, MatrixAddWellContributions, "Explicitly specify the influences of wells between cells in the Jacobian and preconditioner matrices"); - EWOMS_REGISTER_PARAM(TypeTag, bool, PreconditionerAddWellContributions, "Explicitly specify the influences of wells between cells for the preconditioner matrix only"); } }; } // namespace Opm diff --git a/opm/autodiff/BlackoilWellModel.hpp b/opm/autodiff/BlackoilWellModel.hpp index 230651493..ecf5246e8 100644 --- a/opm/autodiff/BlackoilWellModel.hpp +++ b/opm/autodiff/BlackoilWellModel.hpp @@ -229,7 +229,7 @@ namespace Opm { const SimulatorReport& lastReport() const; - void addWellContributions(Mat& mat) + void addWellContributions(Mat& mat) const { for ( const auto& well: well_container_ ) { well->addWellContributions(mat); diff --git a/opm/autodiff/FlowLinearSolverParameters.hpp b/opm/autodiff/FlowLinearSolverParameters.hpp index 4bf8aa11e..8f2e1fe98 100644 --- a/opm/autodiff/FlowLinearSolverParameters.hpp +++ b/opm/autodiff/FlowLinearSolverParameters.hpp @@ -32,6 +32,12 @@ #include #include +namespace Opm { +template +class ISTLSolverEbos; +} + + BEGIN_PROPERTIES NEW_TYPE_TAG(FlowIstlSolverParams); @@ -51,6 +57,8 @@ NEW_PROP_TAG(LinearSolverRequireFullSparsityPattern); NEW_PROP_TAG(LinearSolverIgnoreConvergenceFailure); NEW_PROP_TAG(UseAmg); NEW_PROP_TAG(UseCpr); +NEW_PROP_TAG(LinearSolverBackend); +NEW_PROP_TAG(PreconditionerAddWellContributions); SET_SCALAR_PROP(FlowIstlSolverParams, LinearSolverReduction, 1e-2); SET_SCALAR_PROP(FlowIstlSolverParams, IluRelaxation, 0.9); @@ -66,6 +74,10 @@ SET_BOOL_PROP(FlowIstlSolverParams, LinearSolverRequireFullSparsityPattern, fals SET_BOOL_PROP(FlowIstlSolverParams, LinearSolverIgnoreConvergenceFailure, false); SET_BOOL_PROP(FlowIstlSolverParams, UseAmg, false); SET_BOOL_PROP(FlowIstlSolverParams, UseCpr, false); +SET_TYPE_PROP(FlowIstlSolverParams, LinearSolverBackend, Opm::ISTLSolverEbos); +SET_BOOL_PROP(FlowIstlSolverParams, PreconditionerAddWellContributions, false); + + END_PROPERTIES @@ -186,6 +198,7 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge"); EWOMS_REGISTER_PARAM(TypeTag, bool, UseAmg, "Use AMG as the linear solver's preconditioner"); EWOMS_REGISTER_PARAM(TypeTag, bool, UseCpr, "Use CPR as the linear solver's preconditioner"); + } FlowLinearSolverParameters() { reset(); } diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index 0f2ce66a7..f0dcd2606 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -41,6 +41,7 @@ #include #include +#include #include #include #include @@ -94,7 +95,6 @@ namespace Opm typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef Opm::SimulatorFullyImplicitBlackoilEbos Simulator; - typedef typename BlackoilModelEbos::ISTLSolverType ISTLSolverType; // Read the command line parameters. Throws an exception if something goes wrong. static int setupParameters_(int argc, char** argv) @@ -108,8 +108,6 @@ namespace Opm "Specify the number of report steps between two consecutive writes of restart data"); Simulator::registerParameters(); - ISTLSolverType::registerParameters(); - // register the parameters inherited from ebos Ewoms::registerAllParameters_(/*finalizeRegistration=*/false); @@ -136,14 +134,6 @@ namespace Opm EWOMS_HIDE_PARAM(TypeTag, MinTimeStepSize); EWOMS_HIDE_PARAM(TypeTag, PredeterminedTimeStepsFile); - // flow currently uses its own linear solver - EWOMS_HIDE_PARAM(TypeTag, LinearSolverMaxError); - EWOMS_HIDE_PARAM(TypeTag, LinearSolverMaxIterations); - EWOMS_HIDE_PARAM(TypeTag, LinearSolverOverlapSize); - EWOMS_HIDE_PARAM(TypeTag, LinearSolverTolerance); - EWOMS_HIDE_PARAM(TypeTag, LinearSolverVerbosity); - EWOMS_HIDE_PARAM(TypeTag, PreconditionerRelaxation); - // flow also does not use the eWoms Newton method EWOMS_HIDE_PARAM(TypeTag, NewtonMaxError); EWOMS_HIDE_PARAM(TypeTag, NewtonMaxIterations); @@ -224,7 +214,6 @@ namespace Opm setupLogging(); printPRTHeader(); runDiagnostics(); - setupLinearSolver(); createSimulator(); // do the actual work @@ -546,33 +535,6 @@ namespace Opm } } - // Setup linear solver. - // Writes to: - // linearSolver_ - void setupLinearSolver() - { - extractParallelGridInformationToISTL(grid(), parallel_information_); - auto *tmp = new ISTLSolverType(parallel_information_); - linearSolver_.reset(tmp); - - // Deactivate selection of CPR via eclipse keyword - // as this preconditioner is still considered experimental - // and fails miserably for some models. - if (output_cout_ - && eclState().getSimulationConfig().useCPR() - && !tmp->parameters().use_cpr_) - { - std::ostringstream message; - message << "Ignoring request for CPRPreconditioner " - << "via Eclipse keyword as it is considered " - <<" experimental. To activate use " - <<"\"--flow-use-cpr=true\" command " - <<"line parameter."; - OpmLog::info(message.str()); - } - - } - /// This is the main function of Flow. // Create simulator instance. // Writes to: @@ -580,7 +542,7 @@ namespace Opm void createSimulator() { // Create the simulator instance. - simulator_.reset(new Simulator(*ebosSimulator_, *linearSolver_)); + simulator_.reset(new Simulator(*ebosSimulator_)); } unsigned long long getTotalSystemMemory() @@ -602,7 +564,6 @@ namespace Opm FileOutputMode output_ = OUTPUT_ALL; bool output_to_files_ = false; boost::any parallel_information_; - std::unique_ptr linearSolver_; std::unique_ptr simulator_; std::string logFile_; }; diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index 26b7585a0..660dcc648 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -26,7 +26,8 @@ #include #include #include - +#include +#include #include #include #include @@ -51,11 +52,112 @@ NEW_PROP_TAG(Scalar); NEW_PROP_TAG(GlobalEqVector); NEW_PROP_TAG(SparseMatrixAdapter); NEW_PROP_TAG(Indices); +NEW_PROP_TAG(Simulator); +NEW_PROP_TAG(EclWellModel); END_PROPERTIES namespace Opm { +//===================================================================== +// Implementation for ISTL-matrix based operator +//===================================================================== + +/*! + \brief Adapter to turn a matrix into a linear operator. + + Adapts a matrix to the assembled linear operator interface + */ +template +class WellModelMatrixAdapter : public Dune::AssembledLinearOperator +{ + typedef Dune::AssembledLinearOperator BaseType; + +public: + typedef M matrix_type; + typedef X domain_type; + typedef Y range_type; + typedef typename X::field_type field_type; + +#if HAVE_MPI + typedef Dune::OwnerOverlapCopyCommunication communication_type; +#else + typedef Dune::CollectiveCommunication< int > communication_type; +#endif + +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::SolverCategory::Category category() const override + { + return overlapping ? + Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential; + } +#else + enum { + //! \brief The solver category. + category = overlapping ? + Dune::SolverCategory::overlapping : + Dune::SolverCategory::sequential + }; +#endif + + //! constructor: just store a reference to a matrix + WellModelMatrixAdapter (const M& A, + const M& A_for_precond, + const WellModel& wellMod, + const boost::any& parallelInformation = boost::any() ) + : A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), comm_() + { +#if HAVE_MPI + if( parallelInformation.type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& info = + boost::any_cast( parallelInformation); + comm_.reset( new communication_type( info.communicator() ) ); + } +#endif + } + + virtual void apply( const X& x, Y& y ) const + { + A_.mv( x, y ); + + // add well model modification to y + wellMod_.apply(x, y ); + +#if HAVE_MPI + if( comm_ ) + comm_->project( y ); +#endif + } + + // y += \alpha * A * x + virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const + { + A_.usmv(alpha,x,y); + + // add scaled well model modification to y + wellMod_.applyScaleAdd( alpha, x, y ); + +#if HAVE_MPI + if( comm_ ) + comm_->project( y ); +#endif + } + + virtual const matrix_type& getmat() const { return A_for_precond_; } + + communication_type* comm() + { + return comm_.operator->(); + } + +protected: + const matrix_type& A_ ; + const matrix_type& A_for_precond_ ; + const WellModel& wellMod_; + std::unique_ptr< communication_type > comm_; +}; + /// This class solves the fully implicit black-oil system by /// solving the reduced system (after eliminating well variables) /// as a block-structured matrix (one block for all cell variables) for a fixed @@ -72,10 +174,13 @@ namespace Opm typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - + typedef typename GET_PROP_TYPE(TypeTag, EclWellModel) WellModel; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename SparseMatrixAdapter::IstlMatrix Matrix; enum { pressureIndex = Indices::pressureSwitchIdx }; + static const int numEq = Indices::numEq; + public: typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType; @@ -88,16 +193,57 @@ namespace Opm /// Construct a system solver. /// \param[in] parallelInformation In the case of a parallel run /// with dune-istl the information about the parallelization. - ISTLSolverEbos(const boost::any& parallelInformation_arg=boost::any()) - : iterations_( 0 ) - , parallelInformation_(parallelInformation_arg) - , isIORank_(isIORank(parallelInformation_arg)) + ISTLSolverEbos(const Simulator& simulator) + : simulator_(simulator), + iterations_( 0 ), + converged_(false) { parameters_.template init(); + extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_); + detail::findOverlapRowsAndColumns(simulator_.vanguard().grid(),overlapRowAndColumns_); + } + + // nothing to clean here + void eraseMatrix() { + matrix_for_preconditioner_.reset(); + } + + void prepare(const SparseMatrixAdapter& M, Vector& b) { + matrix_ = &M.istlMatrix(); + rhs_ = &b; + } + + bool solve(Vector& x) { + // Solve system. + + const WellModel& wellModel = simulator_.problem().wellModel(); + + if( isParallel() ) + { + typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator; + + auto ebosJacIgnoreOverlap = Matrix(*matrix_); + //remove ghost rows in local matrix + makeOverlapRowsInvalid(ebosJacIgnoreOverlap); + + //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables + //to be certain that correct matrix is used for preconditioning. + Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel, + parallelInformation_ ); + assert( opA.comm() ); + solve( opA, x, *rhs_, *(opA.comm()) ); + } + else + { + typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false > Operator; + Operator opA(*matrix_, *matrix_, wellModel); + solve( opA, x, *rhs_ ); + } + + return converged_; + } - const FlowLinearSolverParameters& parameters() const - { return parameters_; } /// Solve the system of linear equations Ax = b, with A being the /// combined derivative matrix of the residual and b @@ -111,7 +257,7 @@ namespace Opm /// \copydoc NewtonIterationBlackoilInterface::parallelInformation const boost::any& parallelInformation() const { return parallelInformation_; } - public: + protected: /// \brief construct the CPR preconditioner and the solver. /// \tparam P The type of the parallel information. /// \param parallelInformation the information about the parallelization. @@ -373,6 +519,7 @@ namespace Opm { // store number of iterations iterations_ = result.iterations; + converged_ = result.converged; // Check for failure of linear solver. if (!parameters_.ignoreConvergenceFailure_ && !result.converged) { @@ -381,10 +528,52 @@ namespace Opm } } protected: + + bool isParallel() const { +#if HAVE_MPI + return parallelInformation_.type() == typeid(ParallelISTLInformation); +#else + return false; +#endif + } + + /// Zero out off-diagonal blocks on rows corresponding to overlap cells + /// Diagonal blocks on ovelap rows are set to diag(1e100). + void makeOverlapRowsInvalid(Matrix& ebosJacIgnoreOverlap) const + { + //value to set on diagonal + typedef Dune::FieldMatrix MatrixBlockType; + MatrixBlockType diag_block(0.0); + for (int eq = 0; eq < numEq; ++eq) + diag_block[eq][eq] = 1.0e100; + + //loop over precalculated overlap rows and columns + for (auto row = overlapRowAndColumns_.begin(); row != overlapRowAndColumns_.end(); row++ ) + { + int lcell = row->first; + //diagonal block set to large value diagonal + ebosJacIgnoreOverlap[lcell][lcell] = diag_block; + + //loop over off diagonal blocks in overlap row + for (auto col = row->second.begin(); col != row->second.end(); ++col) + { + int ncell = *col; + //zero out block + ebosJacIgnoreOverlap[lcell][ncell] = 0.0; + } + } + } + + const Simulator& simulator_; mutable int iterations_; + mutable bool converged_; boost::any parallelInformation_; bool isIORank_; + const Matrix *matrix_; + Vector *rhs_; + std::unique_ptr matrix_for_preconditioner_; + std::vector>> overlapRowAndColumns_; FlowLinearSolverParameters parameters_; }; // end ISTLSolver diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 07b5aecdd..2b8434f54 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -75,7 +75,6 @@ public: typedef typename Solver::SolverParameters SolverParameters; typedef BlackoilWellModel WellModel; typedef BlackoilAquiferModel AquiferModel; - typedef typename BlackoilModelEbos::ISTLSolverType ISTLSolverType; /// Initialise from parameters and objects to observe. @@ -99,10 +98,8 @@ public: /// \param[in] eclipse_state the object which represents an internalized ECL deck /// \param[in] output_writer /// \param[in] threshold_pressures_by_face if nonempty, threshold pressures that inhibit flow - SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator, - ISTLSolverType& linearSolver) + SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator) : ebosSimulator_(ebosSimulator) - , linearSolver_(linearSolver) { phaseUsage_ = phaseUsageFromDeck(eclState()); @@ -327,7 +324,6 @@ protected: auto model = std::unique_ptr(new Model(ebosSimulator_, modelParam_, wellModel, - linearSolver_, terminalOutput_)); return std::unique_ptr(new Solver(solverParam_, std::move(model))); @@ -375,7 +371,6 @@ protected: SolverParameters solverParam_; // Observed objects. - ISTLSolverType& linearSolver_; PhaseUsage phaseUsage_; // Misc. data bool terminalOutput_;