OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

[feature] make flow and flow_ebos use the same linear solver setup.

Browse Source
This commit is contained in:
Robert Kloefkorn 2016-11-02 13:10:44 +01:00
parent 01bb7ee4d7
commit 4ff23191eb
6 changed files with 526 additions and 286 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -176,6 +176,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/GridHelpers.hpp
opm/autodiff/GridInit.hpp
opm/autodiff/ImpesTPFAAD.hpp
opm/autodiff/ISTLSolver.hpp
opm/autodiff/moduleVersion.hpp
opm/autodiff/NewtonIterationBlackoilCPR.hpp
opm/autodiff/NewtonIterationBlackoilInterface.hpp

155
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -59,7 +59,7 @@
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#include <dune/istl/solvers.hh>
#include <opm/autodiff/ISTLSolver.hpp>
#include <opm/common/data/SimulationDataContainer.hpp>
#include <cassert>
@ -124,6 +124,8 @@ namespace Opm {
typedef Dune::FieldMatrix<Scalar, 3, 3 > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> BVector;
typedef ISTLSolver< MatrixBlockType, VectorBlockType > ISTLSolverType;
//typedef typename SolutionVector :: value_type PrimaryVariables ;
// --------- Public methods ---------
@ -146,17 +148,17 @@ namespace Opm {
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const StandardWellsDense<FluidSystem, BlackoilIndices>& well_model,
const StandardWellsDense<FluidSystem, BlackoilIndices>& well_model,
const NewtonIterationBlackoilInterface& linsolver,
const bool terminal_output)
: ebosSimulator_(ebosSimulator)
, grid_(ebosSimulator_.gridManager().grid())
, istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) )
, fluid_ (fluid)
, geo_ (geo)
, vfp_properties_(
eclState().getTableManager().getVFPInjTables(),
eclState().getTableManager().getVFPProdTables())
, linsolver_ (linsolver)
, active_(detail::activePhases(fluid.phaseUsage()))
, has_disgas_(FluidSystem::enableDissolvedGas())
, has_vapoil_(FluidSystem::enableVaporizedOil())
@ -174,6 +176,11 @@ namespace Opm {
well_model_.init(&fluid_, &active_, &vfp_properties_, gravity, depth, pv);
wellModel().setWellsActive( localWellsActive() );
global_nc_ = Opm::AutoDiffGrid::numCells(grid_);
if( ! istlSolver_ )
{
OPM_THROW(std::logic_error,"solver down cast to ISTLSolver failed");
}
}
const EclipseState& eclState() const
@ -222,7 +229,6 @@ namespace Opm {
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());
Dune::InverseOperatorResult result;
if (must_solve) {
// enable single precision for solvers when dt is smaller then 20 days
//residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ;
@ -232,7 +238,7 @@ namespace Opm {
const int nw = wellModel().wells().number_of_wells;
BVector x(nc);
BVector xw(nw);
solveJacobianSystem(result, x, xw);
solveJacobianSystem(x, xw);
// Stabilize the nonlinear update.
bool isOscillate = false;
@ -262,7 +268,7 @@ namespace Opm {
isRestart_ = false;
}
const bool failed = false; // Not needed in this model.
const int linear_iters = must_solve ? result.iterations : 0;
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) {
@ -335,16 +341,16 @@ namespace Opm {
}
// compute || u^n - u^n+1 ||
const double stateOld = detail::euclidianNormSquared( p0.begin(), p0.end(), 1, linsolver_.parallelInformation() ) +
const double stateOld = detail::euclidianNormSquared( p0.begin(), p0.end(), 1, istlSolver().parallelInformation() ) +
detail::euclidianNormSquared( sat0.begin(), sat0.end(),
current.numPhases(),
linsolver_.parallelInformation() );
istlSolver().parallelInformation() );
// compute || u^n+1 ||
const double stateNew = detail::euclidianNormSquared( current.pressure().begin(), current.pressure().end(), 1, linsolver_.parallelInformation() ) +
const double stateNew = detail::euclidianNormSquared( current.pressure().begin(), current.pressure().end(), 1, istlSolver().parallelInformation() ) +
detail::euclidianNormSquared( current.saturation().begin(), current.saturation().end(),
current.numPhases(),
linsolver_.parallelInformation() );
istlSolver().parallelInformation() );
if( stateNew > 0.0 ) {
return stateOld / stateNew ;
@ -366,7 +372,7 @@ namespace Opm {
/// Number of linear iterations used in last call to solveJacobianSystem().
int linearIterationsLastSolve() const
{
return linsolver_.iterations();
return istlSolver().iterations();
}
template <class X, class Y>
@ -378,35 +384,33 @@ namespace Opm {
/// Solve the Jacobian system Jx = r where J is the Jacobian and
/// r is the residual.
void solveJacobianSystem(Dune::InverseOperatorResult& result, BVector& x, BVector& xw) const
void solveJacobianSystem(BVector& x, BVector& xw) const
{
typedef double Scalar;
typedef Dune::FieldVector<Scalar, 3 > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, 3, 3 > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> BVector;
const auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
auto& ebosResid = ebosSimulator_.model().linearizer().residual();
typedef WellModelMatrixAdapter<Mat,BVector,BVector, ThisType> Operator;
Operator opA(ebosJac, const_cast< ThisType& > (*this));
const double relax = 0.9;
typedef Dune::SeqILU0<Mat, BVector, BVector> SeqPreconditioner;
SeqPreconditioner precond(opA.getmat(), relax);
Dune::SeqScalarProduct<BVector> sp;
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);
Dune::BiCGSTABSolver<BVector> linsolve(opA, sp, precond,
0.01,
100,
false);
// Solve system.
// set initial guess
x = 0.0;
linsolve.apply(x, ebosResid, result);
typedef typename Operator :: communication_type Comm;
Comm* comm = opA.comm();
// Solve system.
if( comm )
{
istlSolver().solve( opA, x, ebosResid, *comm );
}
else
{
typedef WellModelMatrixAdapter<Mat,BVector,BVector, ThisType> SequentialOperator;
SequentialOperator& sOpA = static_cast< SequentialOperator& > (opA);
istlSolver().solve( sOpA, x, ebosResid );
}
// recover wells.
xw = 0.0;
@ -423,32 +427,84 @@ namespace Opm {
Adapts a matrix to the assembled linear operator interface
*/
template<class M, class X, class Y, class WellModel>
class WellModelMatrixAdapter : public Dune::MatrixAdapter<M,X,Y>
class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
{
typedef Dune::MatrixAdapter<M,X,Y> BaseType;
typedef Dune::AssembledLinearOperator<M,X,Y> BaseType;
public:
//! export types
typedef M matrix_type;
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
enum {
//! \brief The solver category.
category=Dune::SolverCategory::sequential
};
//! constructor: just store a reference to a matrix
explicit WellModelMatrixAdapter (const M& A, WellModel& wellMod ) : BaseType( A ), wellMod_( wellMod ) {}
//! apply operator to x: \f$ y = A(x) \f$
virtual void apply (const X& x, Y& y) const
WellModelMatrixAdapter (const M& A, WellModel& wellMod, const boost::any& parallelInformation )
: A_( A ), wellMod_( wellMod ), comm_()
{
BaseType::apply( x, y );
wellMod_.applyWellModel(x, y );
#if HAVE_MPI
if( parallelInformation.type() == typeid(ParallelISTLInformation) )
{
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>( parallelInformation);
comm_.reset( new communication_type( info.communicator() ) );
}
#endif
}
private:
virtual void apply( const X& x, Y& y ) const
{
A_.mv( x, y );
wellMod_.applyWellModel(x, y );
if( comm_ )
comm_->project( y );
}
virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
{
A_.usmv(alpha,x,y);
wellMod_.applyWellModel(x, y );
if( comm_ )
comm_->project( y );
}
virtual const matrix_type& getmat() const { return A_; }
communication_type* comm()
{
return comm_.operator->();
}
protected:
const matrix_type& A_ ;
WellModel& wellMod_;
std::unique_ptr< communication_type > comm_;
};
/** @} end documentation */
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.
@ -793,16 +849,22 @@ namespace Opm {
const Simulator& ebosSimulator() const
{ return ebosSimulator_; }
protected:
protected:
const ISTLSolverType& istlSolver() const
{
assert( istlSolver_ );
return *istlSolver_;
}
// --------- Data members ---------
Simulator& ebosSimulator_;
const Grid& grid_;
const Grid& grid_;
const ISTLSolverType* istlSolver_;
const BlackoilPropsAdInterface& fluid_;
const DerivedGeology& geo_;
VFPProperties vfp_properties_;
const NewtonIterationBlackoilInterface& linsolver_;
// For each canonical phase -> true if active
const std::vector<bool> active_;
// Size = # active phases. Maps active -> canonical phase indices.
@ -830,7 +892,6 @@ namespace Opm {
public:
/// return the StandardWells object
StandardWellsDense<FluidSystem, BlackoilIndices>& wellModel() { return well_model_; }
const StandardWellsDense<FluidSystem, BlackoilIndices>& wellModel() const { return well_model_; }

6
opm/autodiff/FlowMain.hpp
View File

@ -152,7 +152,7 @@ namespace Opm
asImpl().setupOutputWriter();
asImpl().setupLinearSolver();
asImpl().createSimulator();
// Run.
auto ret = asImpl().runSimulator();
@ -380,12 +380,12 @@ namespace Opm
// Setup OpmLog backend with output_dir.
// Setup OpmLog backend with output_dir.
void setupLogging()
{
std::string deck_filename = param_.get<std::string>("deck_filename");
// create logFile
using boost::filesystem::path;
using boost::filesystem::path;
path fpath(deck_filename);
std::string baseName;
std::ostringstream debugFileStream;

9
opm/autodiff/FlowMainEbos.hpp
View File

@ -77,6 +77,15 @@ namespace Opm
}
}
// Setup linear solver.
// Writes to:
// fis_solver_
void setupLinearSolver()
{
typedef typename BlackoilModelEbos :: ISTLSolverType ISTLSolverType;
Base::fis_solver_.reset( new ISTLSolverType( Base::param_, Base::parallel_information_ ) );
}
/// This is the main function of Flow.
// Create simulator instance.
// Writes to:

394
opm/autodiff/ISTLSolver.hpp Normal file
View File

@ -0,0 +1,394 @@
/*
Copyright 2016 IRIS AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ISTLSOLVER_HEADER_INCLUDED
#define OPM_ISTLSOLVER_HEADER_INCLUDED
#include <opm/autodiff/AdditionalObjectDeleter.hpp>
#include <opm/autodiff/CPRPreconditioner.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
#include <opm/autodiff/NewtonIterationUtilities.hpp>
#include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/amg.hh>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
namespace Dune
{
namespace ISTLUtility {
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,1,1> &matrix)
{
FieldMatrix<K,1,1> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,2,2> &matrix)
{
FieldMatrix<K,2,2> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
{
FieldMatrix<K,3,3> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling matrix.invert
template <typename K, int n>
static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)
{
matrix.invert();
}
} // end ISTLUtility
template <class Scalar, int n, int m>
class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
{
public:
typedef Dune::FieldMatrix<Scalar, n, m> BaseType;
using BaseType :: operator= ;
using BaseType :: rows;
using BaseType :: cols;
explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {}
void invert()
{
ISTLUtility::invertMatrix( *this );
}
const BaseType& asBase() const { return static_cast< const BaseType& > (*this); }
BaseType& asBase() { return static_cast< BaseType& > (*this); }
};
template<class K, int n, int m>
void
print_row (std::ostream& s, const MatrixBlock<K,n,m>& A,
typename FieldMatrix<K,n,m>::size_type I,
typename FieldMatrix<K,n,m>::size_type J,
typename FieldMatrix<K,n,m>::size_type therow, int width,
int precision)
{
print_row(s, A.asBase(), I, J, therow, width, precision);
}
template<class K, int n, int m>
K& firstmatrixelement (MatrixBlock<K,n,m>& A)
{
return firstmatrixelement( A.asBase() );
}
template<typename Scalar, int n, int m>
struct MatrixDimension< MatrixBlock< Scalar, n, m > >
: public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType >
{
};
} // end namespace Dune
namespace Opm
{
/// 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
/// number of cell variables np .
template < class MatrixBlockType, class VectorBlockType >
class ISTLSolver : public NewtonIterationBlackoilInterface
{
typedef typename MatrixBlockType :: field_type Scalar;
typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
typedef Dune::BlockVector<VectorBlockType> Vector;
public:
typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType;
typedef NewtonIterationBlackoilInterface :: SolutionVector SolutionVector;
/// Construct a system solver.
/// \param[in] param parameters controlling the behaviour of the linear solvers
/// \param[in] parallelInformation In the case of a parallel run
/// with dune-istl the information about the parallelization.
ISTLSolver(const NewtonIterationBlackoilInterleavedParameters& param,
const boost::any& parallelInformation_arg=boost::any())
: iterations_( 0 ),
parallelInformation_(parallelInformation_arg),
isIORank_(isIORank(parallelInformation_arg)),
parameters_( param )
{
}
/// Construct a system solver.
/// \param[in] param ParameterGroup controlling the behaviour of the linear solvers
/// \param[in] parallelInformation In the case of a parallel run
/// with dune-istl the information about the parallelization.
ISTLSolver(const parameter::ParameterGroup& param,
const boost::any& parallelInformation_arg=boost::any())
: iterations_( 0 ),
parallelInformation_(parallelInformation_arg),
isIORank_(isIORank(parallelInformation_arg)),
parameters_( param )
{
}
// dummy method that is not implemented for this class
SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
{
OPM_THROW(std::logic_error,"This method is not implemented");
return SolutionVector();
}
/// Solve the system of linear equations Ax = b, with A being the
/// combined derivative matrix of the residual and b
/// being the residual itself.
/// \param[in] residual residual object containing A and b.
/// \return the solution x
/// \copydoc NewtonIterationBlackoilInterface::iterations
int iterations () const { return iterations_; }
/// \copydoc NewtonIterationBlackoilInterface::parallelInformation
const boost::any& parallelInformation() const { return parallelInformation_; }
public:
/// \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>
void constructPreconditionerAndSolve(O& opA,
Vector& x, Vector& istlb,
const POrComm& parallelInformation_arg,
Dune::InverseOperatorResult& result) const
{
// Construct scalar product.
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
// Communicate if parallel.
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
#if ! HAVE_UMFPACK
const bool useAmg = false ;
if( useAmg )
{
typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType;
typedef typename CPRSelectorType::AMG AMG;
std::unique_ptr< AMG > amg;
// Construct preconditioner.
constructAMGPrecond(opA, parallelInformation_arg, amg);
// Solve.
solve(opA, x, istlb, *sp, *amg, result);
}
else
#endif
{
// Construct preconditioner.
auto precond = constructPrecond(opA, parallelInformation_arg);
// Solve.
solve(opA, x, istlb, *sp, *precond, result);
}
}
typedef Dune::SeqILU0<Matrix, Vector, Vector> SeqPreconditioner;
template <class Operator>
std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
{
const double relax = 0.9;
std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), relax));
return precond;
}
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
typedef ParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
template <class Operator>
std::unique_ptr<ParPreconditioner>
constructPrecond(Operator& opA, const Comm& comm) const
{
typedef std::unique_ptr<ParPreconditioner> Pointer;
const double relax = 0.9;
return Pointer(new ParPreconditioner(opA.getmat(), comm, relax));
}
#endif
template <class Operator, class POrComm, class AMG >
void
constructAMGPrecond(Operator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg ) const
{
const double relax = 1.0;
ISTLUtility::createAMGPreconditionerPointer( opA, relax, comm, amg );
}
/// \brief Solve the system using the given preconditioner and scalar product.
template <class Operator, class ScalarProd, class Precond>
void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const
{
// TODO: Revise when linear solvers interface opm-core is done
// Construct linear solver.
// GMRes solver
if ( parameters_.newton_use_gmres_ ) {
Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond,
parameters_.linear_solver_reduction_,
parameters_.linear_solver_restart_,
parameters_.linear_solver_maxiter_,
parameters_.linear_solver_verbosity_);
// Solve system.
linsolve.apply(x, istlb, result);
}
else { // BiCGstab solver
Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond,
parameters_.linear_solver_reduction_,
parameters_.linear_solver_maxiter_,
parameters_.linear_solver_verbosity_);
// Solve system.
linsolve.apply(x, istlb, result);
}
}
/// Solve the linear system Ax = b, with A being the
/// combined derivative matrix of the residual and b
/// being the residual itself.
/// \param[in] A matrix A
/// \param[inout] x solution to be computed x
/// \param[in] b right hand side b
void solve(Matrix& A, Vector& x, Vector& b ) const
{
// Parallel version is deactivated until we figure out how to do it properly.
#if HAVE_MPI
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
{
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
Comm istlComm(info.communicator());
// Construct operator, scalar product and vectors needed.
typedef Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector,Comm> Operator;
Operator opA(A, istlComm);
solve( opA, x, b, istlComm );
}
else
#endif
{
// Construct operator, scalar product and vectors needed.
Dune::MatrixAdapter< Matrix, Vector, Vector> opA( A );
solve( opA, x, b );
}
}
/// Solve the linear system Ax = b, with A being the
/// combined derivative matrix of the residual and b
/// being the residual itself.
/// \param[in] A matrix A
/// \param[inout] x solution to be computed x
/// \param[in] b right hand side b
template <class Operator, class Comm >
void solve(Operator& opA, Vector& x, Vector& b, Comm& comm) const
{
Dune::InverseOperatorResult result;
// Parallel version is deactivated until we figure out how to do it properly.
#if HAVE_MPI
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
{
const size_t size = opA.getmat().N();
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
// As we use a dune-istl with block size np the number of components
// per parallel is only one.
info.copyValuesTo(comm.indexSet(), comm.remoteIndices(),
size, 1);
// Construct operator, scalar product and vectors needed.
constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, b, comm, result);
}
else
#endif
{
OPM_THROW(std::logic_error,"this method if for parallel solve only");
}
checkConvergence( result );
}
/// Solve the linear system Ax = b, with A being the
/// combined derivative matrix of the residual and b
/// being the residual itself.
/// \param[in] A matrix A
/// \param[inout] x solution to be computed x
/// \param[in] b right hand side b
template <class Operator>
void solve(Operator& opA, Vector& x, Vector& b ) const
{
Dune::InverseOperatorResult result;
// Construct operator, scalar product and vectors needed.
Dune::Amg::SequentialInformation info;
constructPreconditionerAndSolve(opA, x, b, info, result);
checkConvergence( result );
}
void checkConvergence( const Dune::InverseOperatorResult& result ) const
{
// store number of iterations
iterations_ = result.iterations;
// Check for failure of linear solver.
if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
const std::string msg("Convergence failure for linear solver.");
if (isIORank_) {
OpmLog::problem(msg);
}
OPM_THROW_NOLOG(LinearSolverProblem, msg);
}
}
protected:
mutable int iterations_;
boost::any parallelInformation_;
bool isIORank_;
NewtonIterationBlackoilInterleavedParameters parameters_;
}; // end ISTLSolver
} // namespace Opm
#endif

247
opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
View File

@ -34,13 +34,9 @@
#include <opm/common/Exceptions.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/autodiff/ISTLSolver.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/amg.hh>
#if HAVE_UMFPACK
#include <Eigen/UmfPackSupport>
@ -49,90 +45,6 @@
#endif
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
namespace Dune
{
namespace ISTLUtility {
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,1,1> &matrix)
{
FieldMatrix<K,1,1> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,2,2> &matrix)
{
FieldMatrix<K,2,2> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling FMatrixHelp::invert
template <typename K>
static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
{
FieldMatrix<K,3,3> A ( matrix );
FMatrixHelp::invertMatrix(A, matrix );
}
//! invert matrix by calling matrix.invert
template <typename K, int n>
static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)
{
matrix.invert();
}
} // end ISTLUtility
template <class Scalar, int n, int m>
class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
{
public:
typedef Dune::FieldMatrix<Scalar, n, m> BaseType;
using BaseType :: operator= ;
using BaseType :: rows;
using BaseType :: cols;
explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {}
void invert()
{
ISTLUtility::invertMatrix( *this );
}
const BaseType& asBase() const { return static_cast< const BaseType& > (*this); }
BaseType& asBase() { return static_cast< BaseType& > (*this); }
};
template<class K, int n, int m>
void
print_row (std::ostream& s, const MatrixBlock<K,n,m>& A,
typename FieldMatrix<K,n,m>::size_type I,
typename FieldMatrix<K,n,m>::size_type J,
typename FieldMatrix<K,n,m>::size_type therow, int width,
int precision)
{
print_row(s, A.asBase(), I, J, therow, width, precision);
}
template<class K, int n, int m>
K& firstmatrixelement (MatrixBlock<K,n,m>& A)
{
return firstmatrixelement( A.asBase() );
}
template<typename Scalar, int n, int m>
struct MatrixDimension< MatrixBlock< Scalar, n, m > >
: public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType >
{
};
} // end namespace Dune
namespace Opm
{
@ -161,9 +73,12 @@ namespace Opm
typedef Dune::FieldVector<Scalar, np > VectorBlockType;
typedef Dune::MatrixBlock<Scalar, np, np > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> Vector;
typedef Opm::ISTLSolver< MatrixBlockType, VectorBlockType > ISTLSolverType;
public:
typedef NewtonIterationBlackoilInterface :: SolutionVector SolutionVector;
/// Construct a system solver.
@ -172,9 +87,7 @@ namespace Opm
/// with dune-istl the information about the parallelization.
NewtonIterationBlackoilInterleavedImpl(const NewtonIterationBlackoilInterleavedParameters& param,
const boost::any& parallelInformation_arg=boost::any())
: iterations_( 0 ),
parallelInformation_(parallelInformation_arg),
isIORank_(isIORank(parallelInformation_arg)),
: istlSolver_( param, parallelInformation_arg ),
parameters_( param )
{
}
@ -186,110 +99,12 @@ namespace Opm
/// \return the solution x
/// \copydoc NewtonIterationBlackoilInterface::iterations
int iterations () const { return iterations_; }
int iterations () const { return istlSolver_.iterations(); }
/// \copydoc NewtonIterationBlackoilInterface::parallelInformation
const boost::any& parallelInformation() const { return parallelInformation_; }
const boost::any& parallelInformation() const { return istlSolver_.parallelInformation(); }
public:
/// \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>
void constructPreconditionerAndSolve(O& opA,
Vector& x, Vector& istlb,
const POrComm& parallelInformation_arg,
Dune::InverseOperatorResult& result) const
{
// Construct scalar product.
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
// Communicate if parallel.
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
#if ! HAVE_UMFPACK
const bool useAmg = false ;
if( useAmg )
{
typedef ISTLUtility::CPRSelector< Mat, Vector, Vector, POrComm> CPRSelectorType;
typedef typename CPRSelectorType::AMG AMG;
std::unique_ptr< AMG > amg;
// Construct preconditioner.
constructAMGPrecond(opA, parallelInformation_arg, amg);
// Solve.
solve(opA, x, istlb, *sp, *amg, result);
}
else
#endif
{
// Construct preconditioner.
auto precond = constructPrecond(opA, parallelInformation_arg);
// Solve.
solve(opA, x, istlb, *sp, *precond, result);
}
}
typedef Dune::SeqILU0<Mat, Vector, Vector> SeqPreconditioner;
template <class Operator>
std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
{
const double relax = 0.9;
std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), relax));
return precond;
}
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
typedef ParallelOverlappingILU0<Mat,Vector,Vector,Comm> ParPreconditioner;
template <class Operator>
std::unique_ptr<ParPreconditioner>
constructPrecond(Operator& opA, const Comm& comm) const
{
typedef std::unique_ptr<ParPreconditioner> Pointer;
const double relax = 0.9;
return Pointer(new ParPreconditioner(opA.getmat(), comm, relax));
}
#endif
template <class Operator, class POrComm, class AMG >
void
constructAMGPrecond(Operator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg ) const
{
const double relax = 1.0;
ISTLUtility::createAMGPreconditionerPointer( opA, relax, comm, amg );
}
/// \brief Solve the system using the given preconditioner and scalar product.
template <class Operator, class ScalarProd, class Precond>
void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const
{
// TODO: Revise when linear solvers interface opm-core is done
// Construct linear solver.
// GMRes solver
if ( parameters_.newton_use_gmres_ ) {
Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond,
parameters_.linear_solver_reduction_,
parameters_.linear_solver_restart_,
parameters_.linear_solver_maxiter_,
parameters_.linear_solver_verbosity_);
// Solve system.
linsolve.apply(x, istlb, result);
}
else { // BiCGstab solver
Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond,
parameters_.linear_solver_reduction_,
parameters_.linear_solver_maxiter_,
parameters_.linear_solver_verbosity_);
// Solve system.
linsolve.apply(x, istlb, result);
}
}
void formInterleavedSystem(const std::vector<LinearisedBlackoilResidual::ADB>& eqs,
Mat& istlA) const
{
@ -455,45 +270,8 @@ namespace Opm
Vector x(istlA.M());
x = 0.0;
Dune::InverseOperatorResult result;
// Parallel version is deactivated until we figure out how to do it properly.
#if HAVE_MPI
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
{
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
Comm istlComm(info.communicator());
// As we use a dune-istl with block size np the number of components
// per parallel is only one.
info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices(),
size, 1);
// Construct operator, scalar product and vectors needed.
typedef Dune::OverlappingSchwarzOperator<Mat,Vector,Vector,Comm> Operator;
Operator opA(istlA, istlComm);
constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, istlb, istlComm, result);
}
else
#endif
{
// Construct operator, scalar product and vectors needed.
typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator;
Operator opA(istlA);
Dune::Amg::SequentialInformation info;
constructPreconditionerAndSolve(opA, x, istlb, info, result);
}
// store number of iterations
iterations_ = result.iterations;
// Check for failure of linear solver.
if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
const std::string msg("Convergence failure for linear solver.");
if (isIORank_) {
OpmLog::problem(msg);
}
OPM_THROW_NOLOG(LinearSolverProblem, msg);
}
// solve linear system using ISTL methods
istlSolver_.solve( istlA, x, istlb );
// Copy solver output to dx.
for (int i = 0; i < size; ++i) {
@ -512,10 +290,7 @@ namespace Opm
}
protected:
mutable int iterations_;
boost::any parallelInformation_;
bool isIORank_;
ISTLSolverType istlSolver_;
NewtonIterationBlackoilInterleavedParameters parameters_;
}; // end NewtonIterationBlackoilInterleavedImpl