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

Added parallel CPR that reuses the matrix hierarchy.

Browse Source 
Currently this just compiles but still segfaults in parallel.
This commit is contained in:
Markus Blatt 2019-04-09 14:56:00 +02:00
parent 07e495c9da
commit 6538b59d9e
3 changed files with 186 additions and 180 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
flow/flow_blackoil_dunecpr.cpp
View File

@ -27,8 +27,6 @@
BEGIN_PROPERTIES
NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem));
NEW_PROP_TAG(FluidState);
NEW_PROP_TAG(CprSmootherFine);
NEW_PROP_TAG(CprSmootherCoarse);
//SET_TYPE_PROP(EclBaseProblem, Problem, Ewoms::EclProblem<TypeTag>);
SET_PROP(EclFlowProblemSimple, FluidState)
{
@ -47,33 +45,6 @@ SET_PROP(EclFlowProblemSimple, FluidState)
//typedef Opm::BlackOilFluidSystemSimple<Scalar> type;
typedef Opm::BlackOilFluidState<Evaluation, FluidSystem, enableTemperature, enableEnergy, compositionSwitchEnabled, Indices::numPhases > type;
};
SET_PROP(EclFlowProblemSimple, CprSmootherFine)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
typedef typename SparseMatrixAdapter::IstlMatrix Matrix;
typedef Dune::Amg::SequentialInformation POrComm;
public:
typedef Opm::ParallelOverlappingILU0<Matrix,Vector,Vector, POrComm> type;
//typedef Dune::SeqILU0<Matrix,Vector,Vector, POrComm> type;
};
SET_PROP(EclFlowProblemSimple, CprSmootherCoarse)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
typedef typename SparseMatrixAdapter::IstlMatrix Matrix;
typedef Dune::Amg::SequentialInformation POrComm;
public:
typedef Opm::ParallelOverlappingILU0<Matrix,Vector,Vector, POrComm> type;
//typedef Dune::SeqILU0<Matrix,Vector,Vector, POrComm> type;
};
SET_BOOL_PROP(EclFlowProblemSimple,MatrixAddWellContributions,true);
SET_INT_PROP(EclFlowProblemSimple,LinearSolverVerbosity,1);

1
opm/autodiff/ISTLSolverEbos.hpp
View File

@ -76,6 +76,7 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M)
// Implementation for ISTL-matrix based operator
//=====================================================================
/*!
\brief Adapter to turn a matrix into a linear operator.

158
opm/autodiff/ISTLSolverEbosCpr.hpp
View File

@ -24,10 +24,6 @@
#include <opm/autodiff/BlackoilAmgCpr.hpp>
#include <utility>
#include <memory>
BEGIN_PROPERTIES
NEW_PROP_TAG(CprSmootherFine);
NEW_PROP_TAG(CprSmootherCoarse);
END_PROPERTIES
namespace Opm
{
@ -60,23 +56,25 @@ namespace Opm
enum { pressureVarIndex = SuperClass::pressureVarIndex };
// New properties in this subclass.
using CprSmootherFine = typename GET_PROP_TYPE(TypeTag, CprSmootherFine);
using CprSmootherCoarse = typename GET_PROP_TYPE(TypeTag, CprSmootherCoarse);
using OperatorSerial = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false>;
using POrComm = Dune::Amg::SequentialInformation;
using Preconditioner = Dune::Preconditioner<Vector, Vector>;
using MatrixAdapter = Dune::MatrixAdapter<Matrix,Vector, Vector>;
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
#else
static constexpr int category = Dune::SolverCategory::sequential;
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
#endif
using CouplingMetric = Opm::Amg::Element<pressureEqnIndex,pressureVarIndex>;
using CritBase = Dune::Amg::SymmetricCriterion<Matrix, CouplingMetric>;
using Criterion = Dune::Amg::CoarsenCriterion<CritBase>;
using BlackoilAmgType = BlackoilAmgCpr<MatrixAdapter,CprSmootherFine, CprSmootherCoarse, Criterion, POrComm, pressureEqnIndex, pressureVarIndex>;
using ParallelMatrixAdapter = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Dune::OwnerOverlapCopyCommunication<int,int> >;
using CprSmootherFine = Opm::ParallelOverlappingILU0<Matrix, Vector, Vector, Dune::Amg::SequentialInformation>;
using CprSmootherCoarse = CprSmootherFine;
using BlackoilAmgType = BlackoilAmgCpr<MatrixAdapter,CprSmootherFine, CprSmootherCoarse, Criterion, Dune::Amg::SequentialInformation,
pressureEqnIndex, pressureVarIndex>;
using ParallelCprSmootherFine = Opm::ParallelOverlappingILU0<Matrix, Vector, Vector, Dune::OwnerOverlapCopyCommunication<int,int> >;
using ParallelCprSmootherCoarse = ParallelCprSmootherFine;
using ParallelBlackoilAmgType = BlackoilAmgCpr<ParallelMatrixAdapter, ParallelCprSmootherFine, ParallelCprSmootherCoarse, Criterion,
Dune::OwnerOverlapCopyCommunication<int,int>, pressureEqnIndex, pressureVarIndex>;
using OperatorSerial = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false>;
using OperatorParallel = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true>;
public:
static void registerParameters()
@ -88,12 +86,17 @@ namespace Opm
/// \param[in] parallelInformation In the case of a parallel run
/// with dune-istl the information about the parallelization.
explicit ISTLSolverEbosCpr(const Simulator& simulator)
: SuperClass(simulator)
: SuperClass(simulator), oldMat()
{
extractParallelGridInformationToISTL(this->simulator_.vanguard().grid(), this->parallelInformation_);
detail::findOverlapRowsAndColumns(this->simulator_.vanguard().grid(), this->overlapRowAndColumns_);
}
void prepare(const SparseMatrixAdapter& M, Vector& b)
{
if (oldMat != nullptr)
std::cout << "old was "<<oldMat<<" new is "<<&M.istlMatrix()<<std::endl;
oldMat = &M.istlMatrix();
int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) {
SuperClass::matrix_.reset(new Matrix(M.istlMatrix()));
@ -106,41 +109,65 @@ namespace Opm
#if HAVE_MPI
if( this->isParallel() ) {
#if 0
// Copied from superclass.
typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true> Operator;
auto ebosJacIgnoreOverlap = Matrix(*(this->matrix_));
//remove ghost rows in local matrix
this->makeOverlapRowsInvalid(ebosJacIgnoreOverlap);
//remove ghost rows in local matrix without doing a copy.
this->makeOverlapRowsInvalid(*(this->matrix_));
if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) {
//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,
this->parallelInformation_ );
assert( opA.comm() );
//SuperClass::solve( opA, x, *(this->rhs_), *(opA.comm()) );
Dune::OwnerOverlapCopyCommunication<int,int>* comm = opA.comm();
const size_t size = opA.getmat().N();
opAParallel_.reset(new OperatorParallel(*(this->matrix_), *(this->matrix_), wellModel,
this->parallelInformation_ ));
}
typedef OperatorParallel LinearOperator;
using POrComm = Dune::OwnerOverlapCopyCommunication<int,int>;
POrComm* comm = opAParallel_->comm();
if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) {
assert(comm->indexSet().size()==0);
const size_t size = opAParallel_->getmat().N();
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>( this->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(),
info.copyValuesTo(comm->indexSet(), comm->remoteIndices(),
size, 1);
// Construct operator, scalar product and vectors needed.
Dune::InverseOperatorResult result;
SuperClass::constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, *(this->rhs_), *comm, result);
SuperClass::checkConvergence(result);
#endif // if 0
}
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::overlapping;
auto sp = Dune::createScalarProduct<Vector,POrComm>(*comm, category);
sp_ = std::move(sp);
#else
constexpr int category = Dune::SolverCategory::overlapping;
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
SPPointer sp(ScalarProductChooser::construct(info));
sp_ = std::move(sp);
#endif
using AMGOperator = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, POrComm>;
// If clause is always execute as as Linearoperator is WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false|true>;
if( ! std::is_same< LinearOperator, AMGOperator > :: value &&
( newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_) ) ) {
// create new operator in case linear operator and matrix operator differ
opA_.reset( new AMGOperator( opAParallel_->getmat(), *comm ));
}
prepareSolver(*opAParallel_, *comm);
} else
#endif
{
const WellModel& wellModel = this->simulator_.problem().wellModel();
if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) {
opASerial_.reset(new OperatorSerial(*(this->matrix_), *(this->matrix_), wellModel));
typedef Dune::Amg::SequentialInformation POrComm;
}
using POrComm = Dune::Amg::SequentialInformation;
POrComm parallelInformation_arg;
typedef OperatorSerial LinearOperator;
@ -155,18 +182,28 @@ namespace Opm
SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
sp_ = std::move(sp);
#endif
Vector& istlb = *(this->rhs_);
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
if( ! std::is_same< LinearOperator, MatrixAdapter > :: value ) {
// If clause is always execute as as Linearoperator is WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false|true>;
if( ! std::is_same< LinearOperator, MatrixAdapter > :: value &&
( newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_) ) ) {
// create new operator in case linear operator and matrix operator differ
opA_.reset( new MatrixAdapter( opASerial_->getmat()));//, parallelInformation_arg ) );
}
prepareSolver(*opASerial_, parallelInformation_arg);
}
}
template<typename Operator, typename Comm>
void prepareSolver(Operator& wellOpA, Comm& comm)
{
Vector& istlb = *(this->rhs_);
comm.copyOwnerToAll(istlb, istlb);
const double relax = this->parameters_.ilu_relaxation_;
const MILU_VARIANT ilu_milu = this->parameters_.ilu_milu_;
using Matrix = typename MatrixAdapter::matrix_type;
POrComm& comm = parallelInformation_arg;
const int verbosity = ( this->parameters_.cpr_solver_verbose_ &&
comm.communicator().rank()==0 ) ? 1 : 0;
@ -184,8 +221,11 @@ namespace Opm
//criterion.setGamma(1); // //1 V cycle 2 WW
// Since DUNE 2.2 we also need to pass the smoother args instead of steps directly
typedef typename BlackoilAmgType::Smoother Smoother;
typedef typename BlackoilAmgType::Smoother Smoother;
using AmgType = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
BlackoilAmgType, ParallelBlackoilAmgType>::type;
using OperatorType = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
MatrixAdapter, ParallelMatrixAdapter>::type;
typedef typename AmgType::Smoother Smoother;
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
SmootherArgs smootherArgs;
smootherArgs.iterations = 1;
@ -193,7 +233,7 @@ namespace Opm
const Opm::CPRParameter& params(this->parameters_); // strange conversion
ISTLUtility::setILUParameters(smootherArgs, ilu_milu);
MatrixAdapter& opARef = *opA_;
auto& opARef = reinterpret_cast<OperatorType&>(*opA_);
int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
double dt = this->simulator_.timeStepSize();
bool update_preconditioner = false;
@ -213,12 +253,12 @@ namespace Opm
}
if ( update_preconditioner or (amg_== 0) ) {
amg_.reset( new BlackoilAmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) );
amg_.reset( new AmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) );
} else {
if (this->parameters_.cpr_solver_verbose_) {
std::cout << " Only update amg solver " << std::endl;
}
amg_->updatePreconditioner(opARef, smootherArgs, comm);
reinterpret_cast<AmgType*>(amg_.get())->updatePreconditioner(opARef, smootherArgs, comm);
}
// Solve.
//SuperClass::solve(linearOperator, x, istlb, *sp, *amg, result);
@ -229,21 +269,14 @@ namespace Opm
verbosity_linsolve = this->parameters_.linear_solver_verbosity_;
}
LinearOperator& opASerialRef = *opASerial_;
linsolve_.reset(new Dune::BiCGSTABSolver<Vector>(opASerialRef, *sp_, *amg_,
linsolve_.reset(new Dune::BiCGSTABSolver<Vector>(wellOpA, *sp_, *amg_,
this->parameters_.linear_solver_reduction_,
this->parameters_.linear_solver_maxiter_,
verbosity_linsolve));
}
}
bool solve(Vector& x)
{
if (this->isParallel()) {
// for now only call the superclass
bool converged = SuperClass::solve(x);
return converged;
} else {
// Solve system.
Dune::InverseOperatorResult result;
Vector& istlb = *(this->rhs_);
@ -252,7 +285,6 @@ namespace Opm
if (this->parameters_.scale_linear_system_) {
this->scaleSolution(x);
}
}
return this->converged_;
}
@ -264,16 +296,18 @@ namespace Opm
#else
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
#endif
std::unique_ptr< MatrixAdapter > opA_;
///! \brief The dune-istl operator (either serial or parallel
std::unique_ptr< Dune::LinearOperator<Vector, Vector> > opA_;
///! \brief Serial well matrix adapter
std::unique_ptr< OperatorSerial > opASerial_;
std::unique_ptr< BlackoilAmgType > amg_;
///! \brief Parallel well matrix adapter
std::unique_ptr< OperatorParallel > opAParallel_;
///! \brief The preconditoner to use (either serial or parallel CPR with AMG)
std::unique_ptr< Preconditioner > amg_;
SPPointer sp_;
std::shared_ptr< Dune::BiCGSTABSolver<Vector> > linsolve_;
//std::shared_ptr< Dune::LinearOperator<Vector,Vector> > op_;
//std::shared_ptr< Dune::Preconditioner<Vector,Vector> > prec_;
//std::shared_ptr< Dune::ScalarProduct<Vector> > sp_;
//Vector solution_;
const void* oldMat;
}; // end ISTLSolver
} // namespace Opm