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

Remove GhostLast ILU

Browse Source
This commit is contained in:
andrthu
2020-03-08 15:45:27 +01:00
parent c91eb330d3
commit 3d0d0f9624
6 changed files with 73 additions and 475 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
ebos/eclbasevanguard.hh
View File

@@ -78,7 +78,7 @@ NEW_PROP_TAG(EclOutputInterval);
NEW_PROP_TAG(IgnoreKeywords);
NEW_PROP_TAG(EnableExperiments);
NEW_PROP_TAG(EdgeWeightsMethod);
NEW_PROP_TAG(OwnersFirst);
NEW_PROP_TAG(OwnerCellsFirst);
SET_STRING_PROP(EclBaseVanguard, IgnoreKeywords, "");
SET_STRING_PROP(EclBaseVanguard, EclDeckFileName, "");
@@ -87,7 +87,7 @@ SET_BOOL_PROP(EclBaseVanguard, EnableOpmRstFile, false);
SET_BOOL_PROP(EclBaseVanguard, EclStrictParsing, false);
SET_BOOL_PROP(EclBaseVanguard, SchedRestart, true);
SET_INT_PROP(EclBaseVanguard, EdgeWeightsMethod, 1);
SET_BOOL_PROP(EclBaseVanguard, OwnersFirst, false);
SET_BOOL_PROP(EclBaseVanguard, OwnerCellsFirst, false);
END_PROPERTIES
@@ -135,7 +135,7 @@ public:
"When restarting: should we try to initialize wells and groups from historical SCHEDULE section.");
EWOMS_REGISTER_PARAM(TypeTag, int, EdgeWeightsMethod,
"Choose edge-weighing strategy: 0=uniform, 1=trans, 2=log(trans).");
EWOMS_REGISTER_PARAM(TypeTag, bool, OwnersFirst,
EWOMS_REGISTER_PARAM(TypeTag, bool, OwnerCellsFirst,
"Order cells owned by rank before ghost/overlap cells.");
}
@@ -269,7 +269,7 @@ public:
std::string fileName = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName);
edgeWeightsMethod_ = Dune::EdgeWeightMethod(EWOMS_GET_PARAM(TypeTag, int, EdgeWeightsMethod));
ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnersFirst);
ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
// Skip processing of filename if external deck already exists.
if (!externalDeck_)

7
ebos/eclcpgridvanguard.hh
View File

@@ -193,10 +193,6 @@ public:
{
const auto wells = this->schedule().getWellsatEnd();
auto& eclState = static_cast<ParallelEclipseState&>(this->eclState());
const EclipseGrid* eclGrid = nullptr;
try
{
auto& eclState = dynamic_cast<ParallelEclipseState&>(this->eclState());
@@ -209,8 +205,7 @@ public:
PropsCentroidsDataHandle<Dune::CpGrid> handle(*grid_, eclState, eclGrid, this->centroids_,
cartesianIndexMapper());
defunctWellNames_ = std::get<1>(grid_->loadBalance(handle, edgeWeightsMethod, &wells, faceTrans.data()));
//defunctWellNames_ = std::get<1>(grid_->loadBalance(edgeWeightsMethod, ownersFirst, &wells, faceTrans.data()));
defunctWellNames_ = std::get<1>(grid_->loadBalance(handle, edgeWeightsMethod, &wells, faceTrans.data(), ownersFirst));
}
catch(const std::bad_cast& e)
{

5
ebos/eclthresholdpressure.hh
View File

@@ -309,8 +309,9 @@ private:
// store intersection, this might be costly
const auto& intersection = *isIt;
// ignore boundary intersections for now (TODO?)
if (!intersection.neighbor())
if (intersection.boundary())
continue; // ignore boundary intersections for now (TODO?)
else if (!intersection.neighbor()) //processor boundary but not domain boundary
continue;
const auto& inside = intersection.inside();

103
opm/simulators/linalg/ISTLSolverEbos.hpp
View File

@@ -287,85 +287,6 @@ protected:
std::unique_ptr< communication_type > comm_;
};
/*!
\brief parallel ScalarProduct, that assumes ghost cells are ordered last
Implements the Dune ScalarProduct interface
*/
template<class X, class C>
class GhostLastScalarProduct : public Dune::ScalarProduct<X>
{
public:
typedef X domain_type;
typedef typename X::field_type field_type;
typedef typename Dune::FieldTraits<field_type>::real_type real_type;
typedef typename X::size_type size_type;
typedef C communication_type;
enum {
category = Dune::SolverCategory::overlapping
};
GhostLastScalarProduct (const communication_type& com, size_type interiorSize)
: interiorSize_(interiorSize)
{
comm_.reset(new communication_type(com.communicator()) );
}
virtual field_type dot (const X& x, const X& y)
{
field_type result = field_type(0.0);
for (size_type i = 0; i < interiorSize_; ++i)
result += x[i]*(y[i]);
return comm_->communicator().sum(result);
}
virtual real_type norm (const X& x)
{
field_type result = field_type(0.0);
for (size_type i = 0; i < interiorSize_; ++i)
result += x[i].two_norm2();
return static_cast<double>(std::sqrt(comm_->communicator().sum(result)));
}
private:
std::unique_ptr< communication_type > comm_;
size_type interiorSize_;
};
template <class X, class C, int c>
struct GhostLastSPChooser
{
typedef C communication_type;
enum { solverCategory = c};
};
template <class X, class C>
struct GhostLastSPChooser<X,C,Dune::SolverCategory::sequential>
{
typedef Dune::SeqScalarProduct<X> ScalarProduct;
enum { category = Dune::SolverCategory::sequential};
static ScalarProduct* construct(const C& comm, size_t is)
{
return new ScalarProduct();
}
};
template <class X, class C>
struct GhostLastSPChooser<X,C,Dune::SolverCategory::overlapping>
{
typedef GhostLastScalarProduct<X,C> ScalarProduct;
typedef C communication_type;
enum {category = Dune::SolverCategory::overlapping};
static ScalarProduct* construct(const communication_type& comm, size_t is)
{
return new ScalarProduct(comm, is);
}
};
/// 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
@@ -434,24 +355,23 @@ struct GhostLastSPChooser<X,C,Dune::SolverCategory::overlapping>
}
#endif
extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
const auto& gridForConn = simulator_.vanguard().grid();
const auto wellsForConn = simulator_.vanguard().schedule().getWellsatEnd();
const bool useWellConn = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnersFirst);
interiorCellNum_ = detail::numInteriorCells(simulator_.vanguard().grid(), ownersFirst_);
ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), ownersFirst_);
if (!ownersFirst_) {
if (!ownersFirst_ || parameters_.linear_solver_use_amg_ || parameters_.use_cpr_ ) {
detail::setWellConnections(gridForConn, wellsForConn, useWellConn, wellConnectionsGraph_);
detail::findOverlapAndInterior(gridForConn, overlapRows_, interiorRows_);
if (gridForConn.comm().size() > 1) {
noGhostAdjacency();
setGhostsInNoGhost(*noGhostMat_);
}
if (ownersFirst_)
OpmLog::warning("OwnerCellsFirst option is true, but ignored.");
}
}
@@ -532,7 +452,7 @@ struct GhostLastSPChooser<X,C,Dune::SolverCategory::overlapping>
if( isParallel() )
{
if ( ownersFirst_ ) {
if ( ownersFirst_ && (!parameters_.linear_solver_use_amg_ || !parameters_.use_cpr_) ) {
typedef WellModelGhostLastMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator;
Operator opA(*matrix_, *matrix_, wellModel, interiorCellNum_,
parallelInformation_ );
@@ -541,11 +461,13 @@ struct GhostLastSPChooser<X,C,Dune::SolverCategory::overlapping>
solve( opA, x, *rhs_, *(opA.comm()) );
}
else {
typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator;
copyJacToNoGhost(*matrix_, *noGhostMat_);
Operator opA(*noGhostMat_, *noGhostMat_, wellModel,
parallelInformation_ );
assert( opA.comm() );
solve( opA, x, *rhs_, *(opA.comm()) );
@@ -589,13 +511,12 @@ struct GhostLastSPChooser<X,C,Dune::SolverCategory::overlapping>
const POrComm& parallelInformation_arg,
Dune::InverseOperatorResult& result) const
{
// Construct scalar product.
auto sp = Dune::createScalarProduct<Vector,POrComm>(parallelInformation_arg, category);
#if FLOW_SUPPORT_AMG // activate AMG if either flow_ebos is used or UMFPack is not available
if( parameters_.linear_solver_use_amg_ || parameters_.use_cpr_)
{
{
typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType;
typedef typename CPRSelectorType::Operator MatrixOperator;
@@ -665,7 +586,6 @@ struct GhostLastSPChooser<X,C,Dune::SolverCategory::overlapping>
// Solve.
solve(linearOperator, x, istlb, *sp, *precond, result);
#endif
}
}
@@ -695,9 +615,7 @@ struct GhostLastSPChooser<X,C,Dune::SolverCategory::overlapping>
typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
// 3x3 matrix block inversion was unstable from at least 2.3 until and
// including 2.5.0
typedef GhostLastParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
typedef ParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
template <class Operator>
std::unique_ptr<ParPreconditioner>
constructPrecond(Operator& opA, const Comm& comm) const
@@ -1151,7 +1069,6 @@ struct GhostLastSPChooser<X,C,Dune::SolverCategory::overlapping>
Vector *rhs_;
std::unique_ptr<Matrix> matrix_for_preconditioner_;
std::vector<int> overlapRows_;
std::vector<int> interiorRows_;
std::vector<std::set<int>> wellConnectionsGraph_;

410
opm/simulators/linalg/ParallelOverlappingILU0.hpp
View File

@@ -445,8 +445,7 @@ namespace Opm
typedef typename M::block_type block;
// implement left looking variant with stored inverse
size_t row_count = 0;
for (rowiterator i = A.begin(); row_count < interiorSize; ++i, ++row_count)
for (rowiterator i = A.begin(); i.index() < interiorSize; ++i)
{
// coliterator is diagonal after the following loop
coliterator endij=(*i).end(); // end of row i
@@ -694,6 +693,7 @@ public:
comm_(nullptr), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 )
{
interiorSize_ = A.N();
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), n, milu, redblack,
@@ -723,6 +723,7 @@ public:
comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 )
{
interiorSize_ = A.N();
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), n, milu, redblack,
@@ -773,12 +774,46 @@ public:
comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 )
{
interiorSize_ = A.N();
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), 0, milu, redblack,
reorder_sphere );
}
/*! \brief Constructor.
Constructor gets all parameters to operate the prec.
\param A The matrix to operate on.
\param n ILU fill in level (for testing). This does not work in parallel.
\param w The relaxation factor.
\param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
\param interiorSize The number of interior/owner rows in the matrix.
\param redblack Whether to use a red-black ordering.
\param reorder_sphere If true, we start the reordering at a root node.
The vertices on each layer aound it (same distance) are
ordered consecutivly. If false, we preserver the order of
the vertices with the same color.
*/
template<class BlockType, class Alloc>
ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const ParallelInfo& comm,
const field_type w, MILU_VARIANT milu,
size_type interiorSize, bool redblack=false,
bool reorder_sphere=true)
: lower_(),
upper_(),
inv_(),
comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
interiorSize_(interiorSize)
{
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), 0, milu, redblack,
reorder_sphere );
}
/*!
\brief Prepare the preconditioner.
@@ -806,13 +841,15 @@ public:
const size_type iEnd = lower_.rows();
const size_type lastRow = iEnd - 1;
size_type upperLoppStart = iEnd - interiorSize_;
size_type lowerLoopEnd = interiorSize_;
if( iEnd != upper_.rows() )
{
OPM_THROW(std::logic_error,"ILU: number of lower and upper rows must be the same");
}
// lower triangular solve
for( size_type i=0; i<iEnd; ++ i )
for( size_type i=0; i<lowerLoopEnd; ++ i )
{
dblock rhs( md[ i ] );
const size_type rowI = lower_.rows_[ i ];
@@ -826,7 +863,7 @@ public:
mv[ i ] = rhs; // Lii = I
}
for( size_type i=0; i<iEnd; ++ i )
for( size_type i=upperLoppStart; i<iEnd; ++ i )
{
vblock& vBlock = mv[ lastRow - i ];
vblock rhs ( vBlock );
@@ -970,7 +1007,10 @@ protected:
detail::IsPositiveFunctor() );
break;
default:
bilu0_decomposition( *ILU );
if (interiorSize_ == A.N())
bilu0_decomposition( *ILU );
else
detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
break;
}
}
@@ -1082,368 +1122,8 @@ protected:
//! \brief The relaxation factor to use.
const field_type w_;
const bool relaxation_;
};
/// \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 ParallelInfoT>
class GhostLastParallelOverlappingILU0
: public Dune::Preconditioner<Domain,Range>
{
typedef ParallelInfoT ParallelInfo;
public:
//! \brief The matrix type the preconditioner is for.
typedef typename std::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;
typedef typename matrix_type::block_type block_type;
typedef typename matrix_type::size_type size_type;
protected:
struct CRS
{
CRS() : nRows_( 0 ) {}
size_type rows() const { return nRows_; }
size_type nonZeros() const
{
assert( rows_[ rows() ] != size_type(-1) );
return rows_[ rows() ];
}
void resize( const size_type nRows )
{
if( nRows_ != nRows )
{
nRows_ = nRows ;
rows_.resize( nRows_+1, size_type(-1) );
}
}
void reserveAdditional( const size_type nonZeros )
{
const size_type needed = values_.size() + nonZeros ;
if( values_.capacity() < needed )
{
const size_type estimate = needed * 1.1;
values_.reserve( estimate );
cols_.reserve( estimate );
}
}
void push_back( const block_type& value, const size_type index )
{
values_.push_back( value );
cols_.push_back( index );
}
std::vector< size_type > rows_;
std::vector< block_type > values_;
std::vector< size_type > cols_;
size_type nRows_;
};
public:
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
Dune::SolverCategory::Category category() const override
{
return std::is_same<ParallelInfoT, Dune::Amg::SequentialInformation>::value ?
Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
}
#else
// define the category
enum {
//! \brief The category the preconditioner is part of.
category = std::is_same<ParallelInfoT, Dune::Amg::SequentialInformation>::value ?
Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping
};
#endif
/*! \brief Constructor.
Constructor gets all parameters to operate the prec.
\param A The matrix to operate on.
\param n ILU fill in level (for testing). This does not work in parallel.
\param w The relaxation factor.
\param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
\param redblack Whether to use a red-black ordering.
\param reorder_sphere If true, we start the reordering at a root node.
The vertices on each layer aound it (same distance) are
ordered consecutivly. If false, we preserver the order of
the vertices with the same color.
*/
template<class BlockType, class Alloc>
GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const int n, const field_type w,
MILU_VARIANT milu, unsigned interiorSize, bool redblack=false,
bool reorder_sphere=true)
: lower_(),
upper_(),
inv_(),
comm_(nullptr), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
interiorSize_(interiorSize)
{
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), n, milu, redblack,
reorder_sphere );
}
/*! \brief Constructor gets all parameters to operate the prec.
\param A The matrix to operate on.
\param comm communication object, e.g. Dune::OwnerOverlapCopyCommunication
\param n ILU fill in level (for testing). This does not work in parallel.
\param w The relaxation factor.
\param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
\param redblack Whether to use a red-black ordering.
\param reorder_sphere If true, we start the reordering at a root node.
The vertices on each layer aound it (same distance) are
ordered consecutivly. If false, we preserver the order of
the vertices with the same color.
*/
template<class BlockType, class Alloc>
GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const ParallelInfo& comm, const int n, const field_type w,
MILU_VARIANT milu, unsigned interiorSize, bool redblack=false,
bool reorder_sphere=true)
: lower_(),
upper_(),
inv_(),
comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
interiorSize_(interiorSize)
{
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), n, milu, redblack,
reorder_sphere );
}
/*! \brief Constructor.
Constructor gets all parameters to operate the prec.
\param A The matrix to operate on.
\param w The relaxation factor.
\param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
\param redblack Whether to use a red-black ordering.
\param reorder_sphere If true, we start the reordering at a root node.
The vertices on each layer aound it (same distance) are
ordered consecutivly. If false, we preserver the order of
the vertices with the same color.
*/
template<class BlockType, class Alloc>
GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const field_type w, MILU_VARIANT milu, unsigned interiorSize,
bool redblack=false,
bool reorder_sphere=true)
: GhostLastParallelOverlappingILU0( A, 0, w, milu, interiorSize, redblack, reorder_sphere )
{
}
/*! \brief Constructor.
Constructor gets all parameters to operate the prec.
\param A The matrix to operate on.
\param comm communication object, e.g. Dune::OwnerOverlapCopyCommunication
\param w The relaxation factor.
\param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
\param redblack Whether to use a red-black ordering.
\param reorder_sphere If true, we start the reordering at a root node.
The vertices on each layer aound it (same distance) are
ordered consecutivly. If false, we preserver the order of
the vertices with the same color.
*/
template<class BlockType, class Alloc>
GhostLastParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const ParallelInfo& comm, const field_type w,
MILU_VARIANT milu, unsigned interiorSize,
bool redblack=false,
bool reorder_sphere=true)
: lower_(),
upper_(),
inv_(),
comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
interiorSize_(interiorSize)
{
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), 0, milu, redblack,
reorder_sphere );
}
/*!
\brief Prepare the preconditioner.
\copydoc Preconditioner::pre(X&,Y&)
*/
virtual void pre (Domain& x, Range& b) override
{
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) override
{
Range& md = const_cast<Range&>(d);;
// iterator types
typedef typename Range ::block_type dblock;
typedef typename Domain::block_type vblock;
const size_type iEnd = lower_.rows();
const size_type lastRow = iEnd - 1;
size_type interiorStart = iEnd - interiorSize_;
size_type lowerLoopEnd = interiorSize_;
if( iEnd != upper_.rows() )
{
OPM_THROW(std::logic_error,"ILU: number of lower and upper rows must be the same");
}
// lower triangular solve
for( size_type i = 0; i < lowerLoopEnd; ++ i )
{
dblock rhs( md[ i ] );
const size_type rowI = lower_.rows_[ i ];
const size_type rowINext = lower_.rows_[ i+1 ];
for( size_type col = rowI; col < rowINext; ++ col )
{
lower_.values_[ col ].mmv( v[ lower_.cols_[ col ] ], rhs );
}
v[ i ] = rhs; // Lii = I
}
for( size_type i = interiorStart; i < iEnd; ++ i )
{
vblock& vBlock = v[ lastRow - i ];
vblock rhs ( vBlock );
const size_type rowI = upper_.rows_[ i ];
const size_type rowINext = upper_.rows_[ i+1 ];
for( size_type col = rowI; col < rowINext; ++ col )
{
upper_.values_[ col ].mmv( v[ upper_.cols_[ col ] ], rhs );
}
// apply inverse and store result
inv_[ i ].mv( rhs, vBlock);
}
copyOwnerToAll( v );
if( relaxation_ ) {
v *= w_;
}
}
template <class V>
void copyOwnerToAll( V& v ) const
{
if( comm_ ) {
comm_->copyOwnerToAll(v, v);
}
}
/*!
\brief Clean up.
\copydoc Preconditioner::post(X&)
*/
virtual void post (Range& x) override
{
DUNE_UNUSED_PARAMETER(x);
}
protected:
void init( const Matrix& A, const int iluIteration, MILU_VARIANT milu, bool redBlack, bool reorderSpheres )
{
// (For older DUNE versions the communicator might be
// invalid if redistribution in AMG happened on the coarset level.
// Therefore we check for nonzero size
if ( comm_ && comm_->communicator().size()<=0 )
{
if ( A.N() > 0 )
{
OPM_THROW(std::logic_error, "Expected a matrix with zero rows for an invalid communicator.");
}
else
{
// simply set the communicator to null
comm_ = nullptr;
}
}
int ilu_setup_successful = 1;
std::string message;
const int rank = ( comm_ ) ? comm_->communicator().rank() : 0;
std::unique_ptr< Matrix > ILU;
try
{
ILU.reset( new Matrix( A ) );
detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
}
catch (const Dune::MatrixBlockError& error)
{
message = error.what();
std::cerr<<"Exception occured on process " << rank << " during " <<
"setup of ILU0 preconditioner with message: " <<
message<<std::endl;
ilu_setup_successful = 0;
}
// Check whether there was a problem on some process
const bool parallel_failure = comm_ && comm_->communicator().min(ilu_setup_successful) == 0;
const bool local_failure = ilu_setup_successful == 0;
if ( local_failure || parallel_failure )
{
throw Dune::MatrixBlockError();
}
// store ILU in simple CRS format
detail::convertToCRS( *ILU, lower_, upper_, inv_ );
}
protected:
//! \brief The ILU0 decomposition of the matrix.
CRS lower_;
CRS upper_;
std::vector< block_type > inv_;
const ParallelInfo* comm_;
//! \brief The relaxation factor to use.
const field_type w_;
const bool relaxation_;
size_type interiorSize_;
};
} // end namespace Opm
#endif

15
opm/simulators/linalg/findOverlapRowsAndColumns.hpp
View File

@@ -112,11 +112,17 @@ namespace detail
}
}
}
/// \brief If ownerFirst=true, returns the number of interior cells in grid, else just numCells().
///
/// If cells in grid is ordered so that interior/owner cells come before overlap/copy cells, the method
/// returns the number of interior cells numInterior. In the linear solver only the first numInterior rows of
/// the matrix are needed.
template <class Grid>
size_t numInteriorCells(const Grid& grid, bool ownerFirst)
size_t numMatrixRowsToUseInSolver(const Grid& grid, bool ownerFirst)
{
size_t numInterior = 0;
if (!ownerFirst)
if (!ownerFirst || grid.comm().size()==1)
return grid.numCells();
const auto& gridView = grid.leafGridView();
auto elemIt = gridView.template begin<0>();
@@ -124,10 +130,9 @@ namespace detail
// loop over cells in mesh
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
// If cell has partition type not equal to interior save row
if (elem.partitionType() == Dune::InteriorEntity) {
// Count only the interior cells.
if (elemIt->partitionType() == Dune::InteriorEntity) {
numInterior++;
}
}