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

Used the correct parallel information for the whole system.

Browse Source 
The whole system consists out of three equations per cell. Using
the parallel index set of the grid cells for it is wrong. Therefore
we use PR OPM/opm-core#803 to set up an  additional parallel index set
for the whole system and use this for the communication that is needed e.g.
by the ILU preconditioner.
This commit is contained in:
Markus Blatt 2015-05-19 19:41:32 +02:00
parent cb4970c9a6
commit 764e1e26c1
3 changed files with 25 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
opm/autodiff/CPRPreconditioner.hpp
View File

@ -371,28 +371,31 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
parallel run
*/
CPRPreconditioner (const CPRParameter& param, const M& A, const M& Ae,
const ParallelInformation& comm=ParallelInformation())
const ParallelInformation& comm=ParallelInformation(),
const ParallelInformation& commAe=ParallelInformation())
: param_( param ),
A_(A),
Ae_(Ae),
de_( Ae_.N() ),
ve_( Ae_.M() ),
dmodified_( A_.N() ),
opAe_(CPRSelector<M,X,Y,P>::makeOperator(Ae_, comm)),
opAe_(CPRSelector<M,X,Y,P>::makeOperator(Ae_, commAe)),
precond_(), // ilu0 preconditioner for elliptic system
amg_(), // amg preconditioner for elliptic system
pre_(), // copy A will be made be the preconditioner
vilu_( A_.N() ),
comm_(comm)
comm_(comm),
commAe_(commAe)
{
// create appropriate preconditioner for elliptic system
createEllipticPreconditioner( param_.cpr_use_amg_, comm );
createEllipticPreconditioner( param_.cpr_use_amg_, commAe_ );
// create the preconditioner for the whole system.
if( param_.cpr_ilu_n_ == 0 ) {
pre_ = createILU0Ptr<M,X>( A_, param_.cpr_relax_, comm );
pre_ = createILU0Ptr<M,X>( A_, param_.cpr_relax_, comm_ );
}
else {
pre_ = createILUnPtr<M,X>( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm );
pre_ = createILUnPtr<M,X>( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm_ );
}
}
@ -430,6 +433,8 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
// dmodified = d - A * vfull
dmodified_ = d;
A_.mmv(v, dmodified_);
// A is not parallel, do communication manually.
comm_.copyOwnerToAll(dmodified_, dmodified_);
// Apply Preconditioner for whole system (relax will be applied already)
pre_->apply( vilu_, dmodified_);
@ -470,7 +475,7 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
ScalarProductChooser;
// the scalar product.
std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
sp(ScalarProductChooser::construct(comm_));
sp(ScalarProductChooser::construct(commAe_));
if( amg_ )
{
@ -534,8 +539,11 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
//! \brief temporary variables for ILU solve
Y vilu_;
//! \brief The information about the parallelization
//! \brief The information about the parallelization of the whole system.
const P& comm_;
//! \brief The information about the parallelization of the elliptic part
//! of the system
const P& commAe_;
protected:
void createEllipticPreconditioner( const bool amg, const P& comm )
{

9
opm/autodiff/NewtonIterationBlackoilCPR.cpp
View File

@ -188,11 +188,14 @@ namespace Opm
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
Comm istlComm(info.communicator());
info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices());
Comm istlAeComm(info.communicator());
info.copyValuesTo(istlAeComm.indexSet(), istlAeComm.remoteIndices());
info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices(),
istlAe.N(), istlA.N()/istlAe.N());
// 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);
constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, istlAe, x, istlb, istlComm, istlAeComm, result);
}
else
#endif
@ -201,7 +204,7 @@ namespace Opm
typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator;
Operator opA(istlA);
Dune::Amg::SequentialInformation info;
constructPreconditionerAndSolve(opA, istlAe, x, istlb, info, result);
constructPreconditionerAndSolve(opA, istlAe, x, istlb, info, info, result);
}
// store number of iterations

4
opm/autodiff/NewtonIterationBlackoilCPR.hpp
View File

@ -85,6 +85,7 @@ namespace Opm
void constructPreconditionerAndSolve(O& opA, DuneMatrix& istlAe,
Vector& x, Vector& istlb,
const P& parallelInformation,
const P& parallelInformationAe,
Dune::InverseOperatorResult& result) const
{
typedef Dune::ScalarProductChooser<Vector,P,category> ScalarProductChooser;
@ -94,7 +95,8 @@ namespace Opm
// typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
typedef Opm::CPRPreconditioner<Mat,Vector,Vector,P> Preconditioner;
parallelInformation.copyOwnerToAll(istlb, istlb);
Preconditioner precond(cpr_param_, opA.getmat(), istlAe, parallelInformation);
Preconditioner precond(cpr_param_, opA.getmat(), istlAe, parallelInformation,
parallelInformationAe);
// TODO: Revise when linear solvers interface opm-core is done
// Construct linear solver.