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 #547 from blattms/parallel-ilu

Browse Source 
Parallel overlapping ILU0
This commit is contained in:
Atgeirr Flø Rasmussen
2015-12-02 11:47:13 +01:00
parent c91e144a1d cb9d566ec2
commit 241da90531
7 changed files with 436 additions and 43 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
opm/autodiff/CPRPreconditioner.hpp
View File

@@ -48,6 +48,8 @@
#include <opm/common/Exceptions.hpp>
#include <opm/autodiff/AdditionalObjectDeleter.hpp>
#include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
namespace Opm
{
namespace
@@ -95,7 +97,11 @@ struct CPRSelector<M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
/// \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> >
//typedef Dune::BlockPreconditioner<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
//EllipticPreconditioner;
//typedef ParallelRestrictedOverlappingSchwarz<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
//EllipticPreconditioner;
typedef ParallelOverlappingILU0<M,X, X, ParallelInformation>
EllipticPreconditioner;
/// \brief The type of the unique pointer to the preconditioner of the elliptic part.
typedef std::unique_ptr<EllipticPreconditioner,
@@ -252,16 +258,13 @@ typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
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;
typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
::EllipticPreconditioner ParallelPreconditioner;
Dune::SeqILU0<M,X,X>* ilu=new Dune::SeqILU0<M,X,X>(Ae, relax);
//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));
return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax));
}
#endif
} // end namespace

40
opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
View File

@@ -27,6 +27,8 @@
#include <opm/autodiff/AdditionalObjectDeleter.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>
@@ -141,43 +143,15 @@ namespace Opm
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
typedef Dune::BlockPreconditioner<Vector, Vector, Comm, SeqPreconditioner> ParPreconditioner;
typedef ParallelOverlappingILU0<Mat,Vector,Vector,Comm> ParPreconditioner;
template <class Operator>
std::unique_ptr<ParPreconditioner,
AdditionalObjectDeleter<SeqPreconditioner> >
std::unique_ptr<ParPreconditioner>
constructPrecond(Operator& opA, const Comm& comm) const
{
typedef AdditionalObjectDeleter<SeqPreconditioner> Deleter;
typedef std::unique_ptr<ParPreconditioner, Deleter> Pointer;
int ilu_setup_successful = 1;
std::string message;
typedef std::unique_ptr<ParPreconditioner> Pointer;
const double relax = 1.0;
SeqPreconditioner* seq_precond = nullptr;
try {
seq_precond = new SeqPreconditioner(opA.getmat(), relax);
}
catch ( Dune::MatrixBlockError error )
{
message = error.what();
std::cerr<<"Exception occured on process " <<
comm.communicator().rank() << " during " <<
"setup of ILU0 preconditioner with message: " <<
message<<std::endl;
ilu_setup_successful = 0;
}
// Check whether there was a problem on some process
if ( comm.communicator().min(ilu_setup_successful) == 0 )
{
if ( seq_precond ) // not null if constructor succeeded
{
// prevent memory leak
delete seq_precond;
throw Dune::MatrixBlockError();
}
}
return Pointer(new ParPreconditioner(*seq_precond, comm),
Deleter(*seq_precond));
return Pointer(new ParPreconditioner(opA.getmat(), comm, relax));
}
#endif

1
opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
View File

@@ -25,7 +25,6 @@
#ifndef OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED
#define OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>

206
opm/autodiff/ParallelOverlappingILU0.hpp Normal file
View File

@@ -0,0 +1,206 @@
/*
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 Statoil 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_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
#define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
#include <dune/istl/preconditioner.hh>
#include <dune/istl/paamg/smoother.hh>
namespace Opm
{
template<class M, class X, class Y, class C>
class ParallelOverlappingILU0;
} // end namespace Opm
namespace Dune
{
namespace Amg
{
/// \brief Tells AMG how to construct the Opm::ParallelOverlappingILU0 smoother
/// \tparam Matrix The type of the Matrix.
/// \tparam Domain The type of the Vector representing the domain.
/// \tparam Range The type of the Vector representing the range.
/// \tparam ParallelInfo The type of the parallel information object
/// used, e.g. Dune::OwnerOverlapCommunication
template<class Matrix, class Domain, class Range, class ParallelInfo>
struct ConstructionTraits<Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> >
{
typedef Dune::SeqILU0<Matrix,Domain,Range> T;
typedef DefaultParallelConstructionArgs<T,ParallelInfo> Arguments;
typedef ConstructionTraits<T> SeqConstructionTraits;
static inline Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* construct(Arguments& args)
{
return new Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>(args.getMatrix(),
args.getComm(),
args.getArgs().relaxationFactor);
}
static inline void deconstruct(Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* bp)
{
delete bp;
}
};
} // end namespace Amg
} // end namespace Dune
namespace Opm
{
/// \brief A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as
///
/// This preconditioner differs from a ParallelRestrictedOverlappingSchwarz with
/// Dune:SeqILU0 in the follwing way:
/// During apply we make sure that the current residual is consistent (i.e.
/// each process knows the same value for each index. The we solve
/// Ly= d for y and make y consistent again. Last we solve Ux = y and
/// make sure that x is consistent.
/// In contrast for ParallelRestrictedOverlappingSchwarz we solve (LU)x = d for x
/// without forcing consistency between the two steps.
/// \tparam Matrix The type of the Matrix.
/// \tparam Domain The type of the Vector representing the domain.
/// \tparam Range The type of the Vector representing the range.
/// \tparam ParallelInfo The type of the parallel information object
/// used, e.g. Dune::OwnerOverlapCommunication
template<class Matrix, class Domain, class Range, class ParallelInfo>
class ParallelOverlappingILU0
: public Dune::Preconditioner<Domain,Range> {
public:
//! \brief The matrix type the preconditioner is for.
typedef typename Dune::remove_const<Matrix>::type matrix_type;
//! \brief The domain type of the preconditioner.
typedef Domain domain_type;
//! \brief The range type of the preconditioner.
typedef Range range_type;
//! \brief The field type of the preconditioner.
typedef typename Domain::field_type field_type;
// define the category
enum {
//! \brief The category the preconditioner is part of.
category=Dune::SolverCategory::overlapping
};
/*! \brief Constructor.
Constructor gets all parameters to operate the prec.
\param A The matrix to operate on.
\param w The relaxation factor.
*/
ParallelOverlappingILU0 (const Matrix& A, const ParallelInfo& comm,
field_type w)
: ilu_(A), comm_(comm), w_(w)
{
int ilu_setup_successful = 1;
std::string message;
try
{
bilu0_decomposition(ilu_);
}
catch ( Dune::MatrixBlockError error )
{
message = error.what();
std::cerr<<"Exception occured on process " <<
comm.communicator().rank() << " during " <<
"setup of ILU0 preconditioner with message: " <<
message<<std::endl;
ilu_setup_successful = 0;
}
// Check whether there was a problem on some process
if ( comm.communicator().min(ilu_setup_successful) == 0 )
{
throw Dune::MatrixBlockError();
}
}
/*!
\brief Prepare the preconditioner.
\copydoc Preconditioner::pre(X&,Y&)
*/
virtual void pre (Domain& x, Range& b)
{
DUNE_UNUSED_PARAMETER(x);
DUNE_UNUSED_PARAMETER(b);
}
/*!
\brief Apply the preconditoner.
\copydoc Preconditioner::apply(X&,const Y&)
*/
virtual void apply (Domain& v, const Range& d)
{
Range& md = const_cast<Range&>(d);
comm_.copyOwnerToAll(md,md);
auto endrow=ilu_.end();
for ( auto row = ilu_.begin(); row != endrow; ++row )
{
auto rhs(d[row.index()]);
for ( auto col = row->begin(); col.index() < row.index(); ++col )
{
col->mmv(v[col.index()],rhs);
}
v[row.index()] = rhs;
}
comm_.copyOwnerToAll(v, v);
auto rendrow = ilu_.beforeBegin();
for( auto row = ilu_.beforeEnd(); row != rendrow; --row)
{
auto rhs(v[row.index()]);
auto col = row->beforeEnd();
for( ; col.index() > row.index(); --col)
{
col->mmv(v[col.index()], rhs);
}
v[row.index()] = 0;
col->umv(rhs, v[row.index()]);
}
comm_.copyOwnerToAll(v, v);
v *= w_;
}
/*!
\brief Clean up.
\copydoc Preconditioner::post(X&)
*/
virtual void post (Range& x)
{
DUNE_UNUSED_PARAMETER(x);
}
private:
//! \brief The ILU0 decomposition of the matrix.
Matrix ilu_;
const ParallelInfo& comm_;
//! \brief The relaxation factor to use.
field_type w_;
};
} // end namespace Opm
#endif

209
opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp Normal file
View File

@@ -0,0 +1,209 @@
/*
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 Statoil 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_PARALLELRESTRICTEDADDITIVESCHWARZ_HEADER_INCLUDED
#define OPM_PARALLELRESTRICTEDADDITIVESCHWARZ_HEADER_INCLUDED
#include <dune/istl/preconditioner.hh>
#include <dune/istl/paamg/smoother.hh>
namespace Opm
{
template<class X, class Y, class C, class T>
class ParallelRestrictedOverlappingSchwarz;
} // end namespace Opm
namespace Dune
{
namespace Amg
{
/// \brief Tells AMG how to construct the Opm::ParallelOverlappingILU0 smoother
/// \tparam Domain The type of the Vector representing the domain.
/// \tparam Range The type of the Vector representing the range.
/// \tparam ParallelInfo The type of the parallel information object
/// used, e.g. Dune::OwnerOverlapCommunication
/// \tparam SeqPreconditioner The underlying sequential preconditioner to use.
template<class Range, class Domain, class ParallelInfo, class SeqPreconditioner>
struct ConstructionTraits<Opm::ParallelRestrictedOverlappingSchwarz<Range,
Domain,
ParallelInfo,
SeqPreconditioner> >
{
typedef DefaultParallelConstructionArgs<SeqPreconditioner,ParallelInfo> Arguments;
typedef ConstructionTraits<SeqPreconditioner> SeqConstructionTraits;
/// \brief Construct a parallel restricted overlapping schwarz preconditioner.
static inline Opm::ParallelRestrictedOverlappingSchwarz<Range,
Domain,
ParallelInfo,
SeqPreconditioner>*
construct(Arguments& args)
{
return new Opm::ParallelRestrictedOverlappingSchwarz
<Range,Domain,ParallelInfo,SeqPreconditioner>(*SeqConstructionTraits
::construct(args),
args.getComm());
}
/// \brief Deconstruct and free a parallel restricted overlapping schwarz preconditioner.
static inline void deconstruct(Opm::ParallelRestrictedOverlappingSchwarz
<Range,Domain,ParallelInfo,SeqPreconditioner>* bp)
{
SeqConstructionTraits
::deconstruct(static_cast<SeqPreconditioner*>(&bp->preconditioner));
delete bp;
}
};
/// \brief Tells AMG how to use Opm::ParallelOverlappingILU0 smoother
/// \tparam Domain The type of the Vector representing the domain.
/// \tparam Range The type of the Vector representing the range.
/// \tparam ParallelInfo The type of the parallel information object
/// used, e.g. Dune::OwnerOverlapCommunication
/// \tparam SeqPreconditioner The underlying sequential preconditioner to use.
template<class Range, class Domain, class ParallelInfo, class SeqPreconditioner>
struct SmootherTraits<Opm::ParallelRestrictedOverlappingSchwarz<Range,
Domain,
ParallelInfo,
SeqPreconditioner> >
{
typedef DefaultSmootherArgs<typename SeqPreconditioner::matrix_type::field_type> Arguments;
};
} // end namespace Amg
} // end namespace Dune
namespace Opm{
/// \brief Block parallel preconditioner.
///
/// This is essentially a wrapper that takes a sequential
/// preconditioner. In each step the sequential preconditioner
/// is applied to the whole subdomain and then all owner data
/// points are updated on all other processes from the processor
/// that knows the complete matrix row for this data point (in dune-istl
/// speak that is the one that owns the data).
///
/// Note that this is different from the usual approach in dune-istl where
/// the application of the sequential preconditioner only takes place on
/// the (owner) partition of the process disregarding any overlap/ghost region.
///
/// For more information see https://www.cs.colorado.edu/~cai/papers/rash.pdf
///
/// \tparam Domain The type of the Vector representing the domain.
/// \tparam Range The type of the Vector representing the range.
/// \tparam ParallelInfo The type of the parallel information object
/// used, e.g. Dune::OwnerOverlapCommunication
/// \tparam SeqPreconditioner The underlying sequential preconditioner to use.
template<class Range, class Domain, class ParallelInfo, class SeqPreconditioner=Dune::Preconditioner<Range,Domain> >
class ParallelRestrictedOverlappingSchwarz
: public Dune::Preconditioner<Range,Domain> {
friend class Dune::Amg
::ConstructionTraits<ParallelRestrictedOverlappingSchwarz<Range,
Domain,
ParallelInfo,
SeqPreconditioner> >;
public:
//! \brief The domain type of the preconditioner.
typedef Domain domain_type;
//! \brief The range type of the preconditioner.
typedef Range range_type;
//! \brief The field type of the preconditioner.
typedef typename Domain::field_type field_type;
//! \brief The type of the communication object.
typedef ParallelInfo communication_type;
// define the category
enum {
//! \brief The category the precondtioner is part of.
category=Dune::SolverCategory::overlapping
};
/*! \brief Constructor.
constructor gets all parameters to operate the prec.
\param p The sequential preconditioner.
\param c The communication object for syncing overlap and copy
data points. (E.~g. OwnerOverlapCommunication )
*/
ParallelRestrictedOverlappingSchwarz (SeqPreconditioner& p, const communication_type& c)
: preconditioner_(p), communication_(c)
{ }
/*!
\brief Prepare the preconditioner.
\copydoc Preconditioner::pre(X&,Y&)
*/
virtual void pre (Domain& x, Range& b)
{
communication_.copyOwnerToAll(x,x); // make dirichlet values consistent
preconditioner_.pre(x,b);
}
/*!
\brief Apply the preconditioner
\copydoc Preconditioner::apply(X&,const Y&)
*/
virtual void apply (Domain& v, const Range& d)
{
apply<true>(v, d);
}
template<bool forward>
void apply (Domain& v, const Range& d)
{
// hack us a mutable d to prevent copying.
Range& md = const_cast<Range&>(d);
communication_.copyOwnerToAll(md,md);
preconditioner_.template apply<forward>(v,d);
communication_.copyOwnerToAll(v,v);
// Make sure that d is the same as at the beginning of apply.
communication_.project(md);
}
/*!
\brief Clean up.
\copydoc Preconditioner::post(X&)
*/
virtual void post (Range& x)
{
preconditioner_.post(x);
}
private:
//! \brief a sequential preconditioner
SeqPreconditioner& preconditioner_;
//! \brief the communication object
const communication_type& communication_;
};
} // end namespace OPM
#endif

2
opm/autodiff/RedistributeDataHandles.hpp
View File

@@ -344,7 +344,7 @@ void distributeGridAndData( Dune::CpGrid& grid,
global_grid.switchToGlobalView();
// distribute the grid and switch to the distributed view
grid.loadBalance(eclipseState);
grid.loadBalance(eclipseState, geology.transmissibility().data());
grid.switchToDistributedView();
std::vector<int> compressedToCartesianIdx;
Opm::createGlobalCellArray(grid, compressedToCartesianIdx);