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

Merge pull request #292 from blattms/parallel-newton-convergence

Browse Source 
Adds parallel solver support for sim_fibo_ad_cp
This commit is contained in:
Atgeirr Flø Rasmussen
2015-02-12 11:24:15 +01:00
parent 02b682ea52 d1239d0a35
commit d2b3e52472
11 changed files with 532 additions and 74 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@@ -28,6 +28,7 @@
list (APPEND MAIN_SOURCE_FILES
opm/autodiff/BlackoilPropsAd.cpp
opm/autodiff/BlackoilPropsAdInterface.cpp
opm/autodiff/ExtractParallelGridInformationToISTL.cpp
opm/autodiff/NewtonIterationBlackoilCPR.cpp
opm/autodiff/NewtonIterationBlackoilSimple.cpp
opm/autodiff/GridHelpers.cpp
@@ -101,6 +102,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/CPRPreconditioner.hpp
opm/autodiff/fastSparseProduct.hpp
opm/autodiff/DuneMatrix.hpp
opm/autodiff/ExtractParallelGridInformationToISTL.hpp
opm/autodiff/GeoProps.hpp
opm/autodiff/GridHelpers.hpp
opm/autodiff/ImpesTPFAAD.hpp

8
examples/sim_fibo_ad_cp.cpp
View File

@@ -61,6 +61,7 @@
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/autodiff/ExtractParallelGridInformationToISTL.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
@@ -210,10 +211,13 @@ try
// Solver for Newton iterations.
std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
boost::any parallel_information;
Opm::extractParallelGridInformationToISTL(*grid, parallel_information);
if (param.getDefault("use_cpr", true)) {
fis_solver.reset(new NewtonIterationBlackoilCPR(param));
fis_solver.reset(new NewtonIterationBlackoilCPR(param, parallel_information));
} else {
fis_solver.reset(new NewtonIterationBlackoilSimple(param));
fis_solver.reset(new NewtonIterationBlackoilSimple(param, parallel_information));
}
// Write parameters used for later reference.

296
opm/autodiff/CPRPreconditioner.hpp
View File

@@ -1,6 +1,8 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 IRIS AS.
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
This file is part of the Open Porous Media project (OPM).
@@ -22,6 +24,7 @@
#define OPM_CPRPRECONDITIONER_HEADER_INCLUDED
#include <memory>
#include <type_traits>
#include <opm/core/utility/platform_dependent/disable_warnings.h>
@@ -43,8 +46,219 @@
namespace Opm
{
namespace
{
//! \brief A custom deleter for the parallel preconditioners.
//!
//! In dune-istl they hold a reference to the sequential preconditioner.
//! In CPRPreconditioner we use unique_ptr for the memory management.
//! Ergo we need to construct the sequential preconditioner with new and
//! make sure that it gets deleted together with the enclosing parallel
//! preconditioner. Therefore this deleter stores a pointer to it and deletes
//! it during destruction.
template<class PREC>
class ParallelPreconditionerDeleter
{
public:
ParallelPreconditionerDeleter()
: ilu_()
{}
ParallelPreconditionerDeleter(PREC& ilu)
: ilu_(&ilu){}
template<class T>
void operator()(T* pt)
{
delete pt;
delete ilu_;
}
private:
PREC* ilu_;
};
///
/// \brief A traits class for selecting the types of the preconditioner.
///
/// \tparam M The type of the matrix.
/// \tparam X The type of the domain of the linear problem.
/// \tparam Y The type of the range of the linear problem.
/// \tparam P The type of the parallel information.
////
template<class M, class X, class Y, class P>
struct CPRSelector
{
/// \brief The information about the parallelization and communication
typedef Dune::Amg::SequentialInformation ParallelInformation;
/// \brief The operator type;
typedef Dune::MatrixAdapter<M, X, Y> Operator;
/// \brief The type of the preconditioner used for the elliptic part.
typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
/// \brief The type of the unique pointer to the preconditioner of the elliptic part.
typedef std::unique_ptr<EllipticPreconditioner> EllipticPreconditionerPointer;
/// \brief creates an Operator from the matrix
/// \param M The matrix to use.
/// \param p The parallel information to use.
static Operator* makeOperator(const M& m, const P&)
{
return new Operator(m);
}
};
#if HAVE_MPI
/// \copydoc CPRSelector<M,X,X,Y,P>
template<class M, class X, class Y, class I1, class I2>
struct CPRSelector<M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
{
/// \brief The information about the parallelization and communication
typedef Dune::OwnerOverlapCopyCommunication<I1,I2> ParallelInformation;
/// \brief The operator type;
typedef Dune::OverlappingSchwarzOperator<M,X,X,ParallelInformation> Operator;
/// \brief The type of the preconditioner used for the elliptic part.
typedef Dune::BlockPreconditioner<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
EllipticPreconditioner;
/// \brief The type of the unique pointer to the preconditioner of the elliptic part.
typedef std::unique_ptr<EllipticPreconditioner,
ParallelPreconditionerDeleter<Dune::SeqILU0<M,X,X> > >
EllipticPreconditionerPointer;
/// \brief creates an Operator from the matrix
/// \param M The matrix to use.
/// \param p The parallel information to use.
static Operator* makeOperator(const M& m, const ParallelInformation& p)
{
return new Operator(m, p);
}
};
//! \brief Creates the deleter needed for the parallel ILU preconditioners.
//! \tparam ILU The type of the underlying sequential ILU preconditioner.
//! \tparam I1 The global index type.
//! \tparam I2 The local index type.
//! \param ilu A reference to the wrapped preconditioner
//! \param p The parallel information for template parameter deduction.
template<class ILU, class I1, class I2>
ParallelPreconditionerDeleter<ILU>
createParallelDeleter(ILU& ilu, const Dune::OwnerOverlapCopyCommunication<I1,I2>& p)
{
(void) p;
return ParallelPreconditionerDeleter<ILU>(ilu);
}
#endif
//! \brief Creates and initializes a unique pointer to an sequential ILU0 preconditioner.
//! \param A The matrix of the linear system to solve.
//! \param relax The relaxation factor to use.
template<class M, class X>
std::shared_ptr<Dune::SeqILU0<M,X,X> >
createILU0Ptr(const M& A, double relax, const Dune::Amg::SequentialInformation&)
{
return std::shared_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>( A, relax) );
}
//! \brief Creates and initializes a shared pointer to an ILUn preconditioner.
//! \param A The matrix of the linear system to solve.
//! \param ilu_n The n parameter for the extension of the nonzero pattern.
//! \param relax The relaxation factor to use.
template<class M, class X>
std::shared_ptr<Dune::SeqILUn<M,X,X> >
createILUnPtr(const M& A, int ilu_n, double relax, const Dune::Amg::SequentialInformation&)
{
return std::shared_ptr<Dune::SeqILUn<M,X,X> >(new Dune::SeqILUn<M,X,X>( A, ilu_n, relax) );
}
#if HAVE_MPI
template<class ILU, class I1, class I2>
struct SelectParallelILUSharedPtr
{
typedef std::shared_ptr<
Dune::BlockPreconditioner<
typename ILU::range_type,
typename ILU::domain_type,
Dune::OwnerOverlapCopyCommunication<I1,I2>,
ILU
>
> type;
};
//! \brief Creates and initializes a shared pointer to an ILUn preconditioner.
//! \param A The matrix of the linear system to solve.
//! \param relax The relaxation factor to use.
/// \param comm The object describing the parallelization information and communication.
template<class M, class X, class I1, class I2>
typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>::type
createILU0Ptr(const M& A, double relax,
const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
{
typedef Dune::BlockPreconditioner<
X,
X,
Dune::OwnerOverlapCopyCommunication<I1,I2>,
Dune::SeqILU0<M,X,X>
> PointerType;
Dune::SeqILU0<M,X,X>* ilu = new Dune::SeqILU0<M,X,X>(A, relax);
return typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>
::type ( new PointerType(*ilu, comm), createParallelDeleter(*ilu, comm));
}
//! \brief Creates and initializes a shared pointer to an ILUn preconditioner.
//! \param A The matrix of the linear system to solve.
//! \param ilu_n The n parameter for the extension of the nonzero pattern.
//! \param relax The relaxation factor to use.
/// \param comm The object describing the parallelization information and communication.
template<class M, class X, class I1, class I2>
typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
createILUnPtr(const M& A, int ilu_n, double relax,
const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
{
typedef Dune::BlockPreconditioner<
X,
X,
Dune::OwnerOverlapCopyCommunication<I1,I2>,
Dune::SeqILUn<M,X,X>
> PointerType;
Dune::SeqILUn<M,X,X>* ilu = new Dune::SeqILUn<M,X,X>( A, ilu_n, relax);
return typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
(new PointerType(*ilu, comm),createParallelDeleter(*ilu, comm));
}
#endif
/// \brief Creates the elliptic preconditioner (ILU0)
/// \param Ae The matrix of the elliptic system.
/// \param relax The relaxation parameter for ILU0
template<class M, class X=typename M::range_type>
std::unique_ptr<Dune::SeqILU0<M,X,X> >
createEllipticPreconditionerPointer(const M& Ae, double relax,
const Dune::Amg::SequentialInformation&)
{
return std::unique_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>(Ae, relax));
}
#if HAVE_MPI
/// \brief Creates the elliptic preconditioner (ILU0)
/// \param Ae The matrix of the elliptic system.
/// \param relax The relaxation parameter for ILU0.
/// \param comm The object describing the parallelization information and communication.
template<class M, class X=typename M::range_type, class I1, class I2>
typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
::EllipticPreconditionerPointer
createEllipticPreconditionerPointer(const M& Ae, double relax,
const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
{
typedef Dune::BlockPreconditioner<X, X,
Dune::OwnerOverlapCopyCommunication<I1,I2>,
Dune::SeqILU0<M,X,X> >
ParallelPreconditioner;
Dune::SeqILU0<M,X,X>* ilu=new Dune::SeqILU0<M,X,X>(Ae, relax);
typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
::EllipticPreconditionerPointer EllipticPreconditionerPointer;
return EllipticPreconditionerPointer(new ParallelPreconditioner(*ilu, comm),
createParallelDeleter(*ilu, comm));
}
#endif
} // end namespace
/*!
\brief Sequential CPR preconditioner.
\brief CPR preconditioner.
This is a two-stage preconditioner, combining an elliptic-type
partial solution with ILU0 for the whole system.
@@ -52,14 +266,21 @@ namespace Opm
\tparam M The matrix type to operate on
\tparam X Type of the update
\tparam Y Type of the defect
\tparam P Type of the parallel information. If not provided
this will be Dune::Amg::SequentialInformation.
The preconditioner is parallel if this is
Dune::OwnerOverlapCopyCommunication<int,int>
*/
template<class M, class X, class Y>
template<class M, class X, class Y,
class P=Dune::Amg::SequentialInformation>
class CPRPreconditioner : public Dune::Preconditioner<X,Y>
{
// prohibit copying for now
CPRPreconditioner( const CPRPreconditioner& );
public:
//! \brief The type describing the parallel information
typedef P ParallelInformation;
//! \brief The matrix type the preconditioner is for.
typedef typename Dune::remove_const<M>::type matrix_type;
//! \brief The domain type of the preconditioner.
@@ -72,21 +293,28 @@ namespace Opm
// define the category
enum {
//! \brief The category the preconditioner is part of.
category = Dune::SolverCategory::sequential
category = std::is_same<P,Dune::Amg::SequentialInformation>::value?
Dune::SolverCategory::sequential:Dune::SolverCategory::overlapping
};
//! \brief Elliptic Operator
typedef Dune::MatrixAdapter<M,X,X> Operator;
typedef typename CPRSelector<M,X,X,P>::Operator Operator;
//! \brief preconditioner for the whole system (here either ILU(0) or ILU(n)
typedef Dune::Preconditioner<X,X> WholeSystemPreconditioner;
//! \brief ilu-0 preconditioner for the elliptic system
typedef Dune::SeqILU0<M,X,X> EllipticPreconditioner;
//! \brief the ilu-0 preconditioner used the for the elliptic system
typedef typename CPRSelector<M,X,X,P>::EllipticPreconditioner
EllipticPreconditioner;
//! \brief type of the unique pointer to the ilu-0 preconditioner
//! used the for the elliptic system
typedef typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
EllipticPreconditionerPointer;
//! \brief amg preconditioner for the elliptic system
typedef EllipticPreconditioner Smoother;
typedef Dune::Amg::AMG<Operator, X, Smoother> AMG;
typedef Dune::Amg::AMG<Operator, X, Smoother, P> AMG;
/*! \brief Constructor.
@@ -96,32 +324,36 @@ namespace Opm
\param relax The ILU0 relaxation factor.
\param useAMG if true, AMG is used as a preconditioner for the elliptic sub-system, otherwise ilu-0 (default)
\param useBiCG if true, BiCG solver is used (default), otherwise CG solver
\param paralleInformation The information about the parallelization, if this is a
parallel run
*/
CPRPreconditioner (const M& A, const M& Ae, const field_type relax,
const unsigned int ilu_n,
const bool useAMG,
const bool useBiCG )
const bool useBiCG,
const ParallelInformation& comm=ParallelInformation())
: A_(A),
Ae_(Ae),
de_( Ae_.N() ),
ve_( Ae_.M() ),
dmodified_( A_.N() ),
opAe_( Ae_ ),
opAe_(CPRSelector<M,X,Y,P>::makeOperator(Ae_, comm)),
precond_(), // ilu0 preconditioner for elliptic system
amg_(), // amg preconditioner for elliptic system
pre_(), // copy A will be made be the preconditioner
vilu_( A_.N() ),
relax_(relax),
use_bicg_solver_( useBiCG )
use_bicg_solver_( useBiCG ),
comm_(comm)
{
// create appropriate preconditioner for elliptic system
createPreconditioner( useAMG );
createPreconditioner( useAMG, comm );
if( ilu_n == 0 ) {
pre_.reset( new Dune::SeqILU0<M,X,X>( A_, relax_) );
pre_ = createILU0Ptr<M,X>( A_, relax_, comm );
}
else {
pre_.reset( new Dune::SeqILUn<M,X,X>( A_, ilu_n, relax_) );
pre_ = createILUnPtr<M,X>( A_, ilu_n, relax_, comm );
}
}
@@ -193,17 +425,22 @@ namespace Opm
// operator result containing iterations etc.
Dune::InverseOperatorResult result;
// sequential scalar product
Dune::SeqScalarProduct<X> sp;
// the scalar product chooser
typedef Dune::ScalarProductChooser<X,ParallelInformation,category>
ScalarProductChooser;
// the scalar product.
std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
sp(ScalarProductChooser::construct(comm_));
if( amg_ )
{
// Solve system with AMG
if( use_bicg_solver_ ) {
Dune::BiCGSTABSolver<X> linsolve(opAe_, sp, (*amg_), tolerance, maxit, verbosity);
Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
linsolve.apply(x, de, result);
}
else {
Dune::CGSolver<X> linsolve(opAe_, sp, (*amg_), tolerance, maxit, verbosity);
Dune::CGSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
linsolve.apply(x, de, result);
}
}
@@ -212,11 +449,11 @@ namespace Opm
assert( precond_ );
// Solve system with ILU-0
if( use_bicg_solver_ ) {
Dune::BiCGSTABSolver<X> linsolve(opAe_, sp, (*precond_), tolerance, maxit, verbosity);
Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
linsolve.apply(x, de, result);
}
else {
Dune::CGSolver<X> linsolve(opAe_, sp, (*precond_), tolerance, maxit, verbosity);
Dune::CGSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
linsolve.apply(x, de, result);
}
@@ -236,15 +473,20 @@ namespace Opm
Y de_, ve_, dmodified_;
//! \brief elliptic operator
Operator opAe_;
std::unique_ptr<Operator> opAe_;
//! \brief ILU0 preconditioner for the elliptic system
std::unique_ptr< EllipticPreconditioner > precond_;
EllipticPreconditionerPointer precond_;
//! \brief AMG preconditioner with ILU0 smoother
std::unique_ptr< AMG > amg_;
//! \brief The preconditioner for the whole system
std::unique_ptr< WholeSystemPreconditioner > pre_;
//!
//! We have to use a shared_ptr instead of a unique_ptr
//! as we need to use a custom allocator based on dynamic
//! information. But for unique_ptr the type of this deleter
//! has to be available at coompile time.
std::shared_ptr< WholeSystemPreconditioner > pre_;
//! \brief temporary variables for ILU solve
Y vilu_;
@@ -255,8 +497,10 @@ namespace Opm
//! \brief true if ISTL BiCGSTABSolver is used, otherwise ISTL CGSolver is used
const bool use_bicg_solver_;
//! \brief The information about the parallelization
const P& comm_;
protected:
void createPreconditioner( const bool amg )
void createPreconditioner( const bool amg, const P& comm )
{
if( amg )
{
@@ -275,10 +519,12 @@ namespace Opm
criterion.setAlpha(.67);
criterion.setBeta(1.0e-6);
criterion.setMaxLevel(10);
amg_ = std::unique_ptr< AMG > (new AMG(opAe_, criterion, smootherArgs));
amg_ = std::unique_ptr< AMG > (new AMG(*opAe_, criterion, smootherArgs));
}
else
precond_ = std::unique_ptr< EllipticPreconditioner > (new EllipticPreconditioner( Ae_, relax_ ));
{
precond_ = createEllipticPreconditionerPointer<M,X>( Ae_, relax_, comm);
}
}
};

49
opm/autodiff/ExtractParallelGridInformationToISTL.cpp Normal file
View File

@@ -0,0 +1,49 @@
/*
Copyright 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
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/>.
*/
#include <config.h>
#include "ExtractParallelGridInformationToISTL.hpp"
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <dune/common/version.hh>
#include <dune/common/shared_ptr.hh>
namespace Opm
{
#if defined(HAVE_DUNE_CORNERPOINT)
#if defined(HAVE_MPI) && defined(HAVE_DUNE_ISTL) && DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
// Extracts the information about the data decomposition from the grid for dune-istl
void extractParallelGridInformationToISTL(const Dune::CpGrid& grid, boost::any& anyComm)
{
if(grid.comm().size()>1)
{
// this is a parallel run with distributed data.
Dune::CpGrid& mgrid=const_cast<Dune::CpGrid&>(grid);
Dune::CpGrid::ParallelIndexSet& idx=mgrid.getCellIndexSet();
Dune::CpGrid::RemoteIndices& ridx=mgrid.getCellRemoteIndices();
anyComm=boost::any(Opm::ParallelISTLInformation(Dune::stackobject_to_shared_ptr(idx),
Dune::stackobject_to_shared_ptr(ridx),
grid.comm()));
}
}
#else
// Missing support for MPI or dune-istl -> do nothing.
void extractParallelGridInformationToISTL(const Dune::CpGrid&, boost::any&)
{}
#endif //defined(HAVE_MPI) && defined(HAVE_DUNE_ISTL) && DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
#endif //defined(HAVE_DUNE_CORNERPOINT)
} // end namespace Opm

44
opm/autodiff/ExtractParallelGridInformationToISTL.hpp Normal file
View File

@@ -0,0 +1,44 @@
/*
Copyright 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
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_EXTRACTPARALLELGRIDINFORMATIONTOISTL_HEADER_INCLUDED
#define OPM_EXTRACTPARALLELGRIDINFORMATIONTOISTL_HEADER_INCLUDED
#ifdef HAVE_DUNE_CORNERPOINT
#include<dune/grid/CpGrid.hpp>
#include<boost/any.hpp>
namespace Opm
{
/// \brief Extracts the information about the data decomposition from the grid for dune-istl
///
/// In the case that grid is a parallel grid this method will query it to get the information
/// about the data decompoisition and convert it to the format expected by the linear algebra
/// of dune-istl.
/// \warn if there is no support for dune-istl and MPI then this functio does not do anything.
/// \param[in] grid The grid to inspect.
/// \param[out] anyComm The handle to to store the information in. If grid is a parallel grid
/// then this will ecapsulate an instance of ParallelISTLInformation.
void extractParallelGridInformationToISTL(const Dune::CpGrid& grid, boost::any& anyComm);
} // end namespace Opm
#endif //defined(HAVE_DUNE_CORNERPOINT)
#endif // OPM_EXTRACTPARALLELGRIDINFORMATIONTOISTL_HEADER_INCLUDED

63
opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
View File

@@ -1,5 +1,7 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
This file is part of the Open Porous Media project (OPM).
@@ -29,6 +31,7 @@
#include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
@@ -1898,20 +1901,60 @@ namespace {
int nc) const
{
// Do the global reductions
for ( int idx=0; idx<MaxNumPhases; ++idx )
#if HAVE_MPI
if(linsolver_.parallelInformation().type()==typeid(ParallelISTLInformation))
{
if (active_[idx]) {
B_avg[idx] = B.col(idx).sum()/nc;
maxCoeff[idx]=tempV.col(idx).maxCoeff();
R_sum[idx] = R.col(idx).sum();
}
else
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
// Compute the global number of cells and porevolume
std::vector<int> v(nc, 1);
auto nc_and_pv = std::tuple<int, double>(0, 0.0);
auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
Opm::Reduction::makeGlobalSumFunctor<double>());
auto nc_and_pv_containers = std::make_tuple(v, geo_.poreVolume());
info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
for ( int idx=0; idx<MaxNumPhases; ++idx )
{
R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.;
if (active_[idx]) {
auto values = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
auto containers = std::make_tuple(B.col(idx),
tempV.col(idx),
R.col(idx));
auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
Opm::Reduction::makeGlobalMaxFunctor<double>(),
Opm::Reduction::makeGlobalSumFunctor<double>());
info.computeReduction(containers, operators, values);
B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv);
maxCoeff[idx] = std::get<1>(values);
R_sum[idx] = std::get<2>(values);
}
else
{
R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.0;
}
}
// Compute pore volume
return std::get<1>(nc_and_pv);
}
else
#endif
{
for ( int idx=0; idx<MaxNumPhases; ++idx )
{
if (active_[idx]) {
B_avg[idx] = B.col(idx).sum()/nc;
maxCoeff[idx]=tempV.col(idx).maxCoeff();
R_sum[idx] = R.col(idx).sum();
}
else
{
R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.0;
}
}
// Compute total pore volume
return geo_.poreVolume().sum();
}
// Compute total pore volume
return geo_.poreVolume().sum();
}
template<class T>

67
opm/autodiff/NewtonIterationBlackoilCPR.cpp
View File

@@ -1,5 +1,7 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
This file is part of the Open Porous Media project (OPM).
@@ -22,17 +24,14 @@
#include <opm/autodiff/DuneMatrix.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/autodiff/CPRPreconditioner.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/core/utility/platform_dependent/disable_warnings.h>
#include <dune/istl/bvector.hh>
// #include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/io.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/preconditioners.hh>
@@ -57,11 +56,8 @@ namespace Opm
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Dune::FieldVector<double, 1 > VectorBlockType;
typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> Vector;
namespace {
@@ -114,8 +110,9 @@ namespace Opm
/// Construct a system solver.
NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param)
: iterations_( 0 )
NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param,
const boost::any& parallelInformation)
: iterations_( 0 ), parallelInformation_(parallelInformation)
{
cpr_relax_ = param.getDefault("cpr_relax", 1.0);
cpr_ilu_n_ = param.getDefault("cpr_ilu_n", 0);
@@ -186,34 +183,38 @@ namespace Opm
// Create ISTL matrix for elliptic part.
DuneMatrix istlAe( A.topLeftCorner(nc, nc) );
// Construct operator, scalar product and vectors needed.
typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator;
Operator opA(istlA);
Dune::SeqScalarProduct<Vector> sp;
// Right hand side.
Vector istlb(opA.getmat().N());
Vector istlb(istlA.N());
std::copy_n(b.data(), istlb.size(), istlb.begin());
// System solution
Vector x(opA.getmat().M());
Vector x(istlA.M());
x = 0.0;
// Construct preconditioner.
// typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
typedef Opm::CPRPreconditioner<Mat,Vector,Vector> Preconditioner;
Preconditioner precond(istlA, istlAe, cpr_relax_, cpr_ilu_n_, cpr_use_amg_, cpr_use_bicgstab_);
// Construct linear solver.
const double tolerance = 1e-3;
const int maxit = 150;
const int verbosity = 0;
const int restart = 40;
Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond, tolerance, restart, maxit, verbosity);
// Solve system.
Dune::InverseOperatorResult result;
linsolve.apply(x, istlb, result);
#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());
info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices());
// Construct operator, scalar product and vectors needed.
typedef Dune::OverlappingSchwarzOperator<Mat,Vector,Vector,Comm> Operator;
Operator opA(istlA, istlComm);
constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, istlAe, 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, istlAe, x, istlb, info, result);
}
// store number of iterations
iterations_ = result.iterations;
@@ -235,6 +236,10 @@ namespace Opm
return dx;
}
const boost::any& NewtonIterationBlackoilCPR::parallelInformation() const
{
return parallelInformation_;
}

50
opm/autodiff/NewtonIterationBlackoilCPR.hpp
View File

@@ -20,10 +20,14 @@
#ifndef OPM_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED
#define OPM_NEWTONITERATIONBLACKOILCPR_HEADER_INCLUDED
#include <opm/autodiff/DuneMatrix.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/autodiff/CPRPreconditioner.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/bvector.hh>
#include <memory>
namespace Opm
@@ -37,7 +41,12 @@ namespace Opm
/// in Fully Implicit Reservoir Simulations" by Gries et al (SPE 163608).
class NewtonIterationBlackoilCPR : public NewtonIterationBlackoilInterface
{
typedef Dune::FieldVector<double, 1 > VectorBlockType;
typedef Dune::FieldMatrix<double, 1, 1> MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> Vector;
public:
/// Construct a system solver.
/// \param[in] param parameters controlling the behaviour of
/// the preconditioning and choice of
@@ -47,7 +56,10 @@ namespace Opm
/// cpr_ilu_n (default 0) use ILU(n) for preconditioning of the linear system
/// cpr_use_amg (default false) if true, use AMG preconditioner for elliptic part
/// cpr_use_bicgstab (default true) if true, use BiCGStab (else use CG) for elliptic part
NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param);
/// \param[in] parallelInformation In the case of a parallel run
/// with dune-istl the information about the parallelization.
NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param,
const boost::any& parallelInformation=boost::any());
/// Solve the system of linear equations Ax = b, with A being the
/// combined derivative matrix of the residual and b
@@ -58,12 +70,46 @@ namespace Opm
/// \copydoc NewtonIterationBlackoilInterface::iterations
virtual int iterations () const { return iterations_; }
/// \copydoc NewtonIterationBlackoilInterface::parallelInformation
virtual const boost::any& parallelInformation() const;
private:
/// \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 P>
void constructPreconditionerAndSolve(O& opA, DuneMatrix& istlAe,
Vector& x, Vector& istlb,
const P& parallelInformation,
Dune::InverseOperatorResult& result) const
{
typedef Dune::ScalarProductChooser<Vector,P,category> ScalarProductChooser;
std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
sp(ScalarProductChooser::construct(parallelInformation));
// Construct preconditioner.
// typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
typedef Opm::CPRPreconditioner<Mat,Vector,Vector,P> Preconditioner;
Preconditioner precond(opA.getmat(), istlAe, cpr_relax_, cpr_ilu_n_, cpr_use_amg_, cpr_use_bicgstab_, parallelInformation);
// Construct linear solver.
const double tolerance = 1e-3;
const int maxit = 150;
const int verbosity = 0;
const int restart = 40;
Dune::RestartedGMResSolver<Vector> linsolve(opA, *sp, precond, tolerance, restart, maxit, verbosity);
// Solve system.
linsolve.apply(x, istlb, result);
}
mutable int iterations_;
double cpr_relax_;
unsigned int cpr_ilu_n_;
bool cpr_use_amg_;
bool cpr_use_bicgstab_;
const boost::any& parallelInformation_;
};
} // namespace Opm

5
opm/autodiff/NewtonIterationBlackoilInterface.hpp
View File

@@ -23,6 +23,7 @@
#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <boost/any.hpp>
namespace Opm
{
@@ -42,6 +43,10 @@ namespace Opm
/// \return number of linear iterations used during last call of computeNewtonIncrement
virtual int iterations () const = 0;
/// \brief Get the information about the parallelization of the grid.
virtual const boost::any& parallelInformation() const = 0;
};
} // namespace Opm

13
opm/autodiff/NewtonIterationBlackoilSimple.cpp
View File

@@ -29,8 +29,11 @@ namespace Opm
/// Construct a system solver.
/// \param[in] linsolver linear solver to use
NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param)
: iterations_( 0 )
/// \param[in] parallelInformation In the case of a parallel run
/// with dune-istl the information about the parallelization.
NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param,
const boost::any& parallelInformation)
: iterations_( 0 ), parallelInformation_(parallelInformation)
{
linsolver_.reset(new LinearSolverFactory(param));
}
@@ -58,7 +61,7 @@ namespace Opm
Opm::LinearSolverInterface::LinearSolverReport rep
= linsolver_->solve(matr.rows(), matr.nonZeros(),
matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
total_residual.value().data(), dx.data());
total_residual.value().data(), dx.data(), parallelInformation_);
// store iterations
iterations_ = rep.iterations;
@@ -71,5 +74,9 @@ namespace Opm
return dx;
}
const boost::any& NewtonIterationBlackoilSimple::parallelInformation() const
{
return parallelInformation_;
}
} // namespace Opm

9
opm/autodiff/NewtonIterationBlackoilSimple.hpp
View File

@@ -40,7 +40,10 @@ namespace Opm
/// Construct a system solver.
/// \param[in] param parameters controlling the behaviour and
/// choice of linear solver.
NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param);
/// \param[in] parallelInformation In the case of a parallel run
/// with dune-istl the information about the parallelization.
NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param,
const boost::any& parallelInformation=boost::any());
/// Solve the system of linear equations Ax = b, with A being the
/// combined derivative matrix of the residual and b
@@ -52,9 +55,13 @@ namespace Opm
/// \copydoc NewtonIterationBlackoilInterface::iterations
virtual int iterations () const { return iterations_; }
/// \copydoc NewtonIterationBlackoilInterface::parallelInformation
virtual const boost::any& parallelInformation() const;
private:
std::unique_ptr<LinearSolverInterface> linsolver_;
mutable int iterations_;
const boost::any& parallelInformation_;
};
} // namespace Opm