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 #2491 from blattms/flexible-fixes

Browse Source 
Fixes ILU0 when using flexible solvers.
This commit is contained in:
Atgeirr Flø Rasmussen 2020-03-24 17:39:27 +01:00 committed by GitHub
commit 8ef508cc08
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 66 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -24,6 +24,7 @@
#include <opm/simulators/linalg/MatrixBlock.hpp>
#include <opm/simulators/linalg/BlackoilAmg.hpp>
#include <opm/simulators/linalg/CPRPreconditioner.hpp>
#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
#include <opm/simulators/linalg/ParallelRestrictedAdditiveSchwarz.hpp>
#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
#include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
@ -990,29 +991,7 @@ protected:
Vector getQuasiImpesWeights()
{
Matrix& A = *matrix_;
Vector weights(rhs_->size());
BlockVector rhs(0.0);
rhs[pressureVarIndex] = 1;
const auto endi = A.end();
for (auto i = A.begin(); i!=endi; ++i) {
const auto endj = (*i).end();
MatrixBlockType diag_block(0.0);
for (auto j=(*i).begin(); j!=endj; ++j) {
if (i.index() == j.index()) {
diag_block = (*j);
break;
}
}
BlockVector bweights;
auto diag_block_transpose = Opm::transposeDenseMatrix(diag_block);
diag_block_transpose.solve(bweights, rhs);
double abs_max =
*std::max_element(bweights.begin(), bweights.end(), [](double a, double b){ return std::abs(a) < std::abs(b); } );
bweights /= std::abs(abs_max);
weights[i.index()] = bweights;
}
return weights;
return Amg::getQuasiImpesWeights<Matrix,Vector>(*matrix_, pressureVarIndex, /* transpose=*/ true);
}
Vector getSimpleWeights(const BlockVector& rhs)

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

@ -21,6 +21,7 @@
#define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
#include <opm/simulators/linalg/GraphColoring.hpp>
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <dune/common/version.hh>
@ -506,6 +507,9 @@ namespace Opm
typedef typename M :: size_type size_type;
lower.clear();
upper.clear();
inv.clear();
lower.resize( A.N() );
upper.resize( A.N() );
inv.resize( A.N() );
@ -599,7 +603,7 @@ namespace Opm
/// used, e.g. Dune::OwnerOverlapCommunication
template<class Matrix, class Domain, class Range, class ParallelInfoT>
class ParallelOverlappingILU0
: public Dune::Preconditioner<Domain,Range>
: public Dune::PreconditionerWithUpdate<Domain,Range>
{
typedef ParallelInfoT ParallelInfo;
@ -656,6 +660,14 @@ protected:
cols_.push_back( index );
}
void clear()
{
rows_.clear();
values_.clear();
cols_.clear();
nRows_= 0;
}
std::vector< size_type > rows_;
std::vector< block_type > values_;
std::vector< size_type > cols_;
@ -691,13 +703,14 @@ public:
upper_(),
inv_(),
comm_(nullptr), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 )
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{
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,
reorder_sphere );
update();
}
/*! \brief Constructor gets all parameters to operate the prec.
@ -721,13 +734,14 @@ public:
upper_(),
inv_(),
comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 )
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{
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,
reorder_sphere );
update();
}
/*! \brief Constructor.
@ -772,13 +786,14 @@ public:
upper_(),
inv_(),
comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 )
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{
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 );
update();
}
/*! \brief Constructor.
@ -806,12 +821,13 @@ public:
inv_(),
comm_(&comm), w_(w),
relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
interiorSize_(interiorSize)
interiorSize_(interiorSize),
A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
{
// 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 );
update( );
}
/*!
@ -905,15 +921,14 @@ public:
DUNE_UNUSED_PARAMETER(x);
}
protected:
void init( const Matrix& A, const int iluIteration, MILU_VARIANT milu, bool redBlack, bool reorderSpheres )
virtual void update() override
{
// (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 )
if ( A_->N() > 0 )
{
OPM_THROW(std::logic_error, "Expected a matrix with zero rows for an invalid communicator.");
}
@ -930,15 +945,15 @@ protected:
std::unique_ptr< Matrix > ILU;
if ( redBlack )
if ( redBlack_ )
{
using Graph = Dune::Amg::MatrixGraph<const Matrix>;
Graph graph(A);
Graph graph(*A_);
auto colorsTuple = colorVerticesWelshPowell(graph);
const auto& colors = std::get<0>(colorsTuple);
const auto& verticesPerColor = std::get<2>(colorsTuple);
auto noColors = std::get<1>(colorsTuple);
if ( reorderSpheres )
if ( reorderSphere_ )
{
ordering_ = reorderVerticesSpheres(colors, noColors, verticesPerColor,
graph, 0);
@ -959,27 +974,27 @@ protected:
try
{
if( iluIteration == 0 ) {
if( iluIteration_ == 0 ) {
// create ILU-0 decomposition
if ( ordering_.empty() )
{
ILU.reset( new Matrix( A ) );
ILU.reset( new Matrix( *A_ ) );
}
else
{
ILU.reset( new Matrix(A.N(), A.M(), A.nonzeroes(), Matrix::row_wise));
ILU.reset( new Matrix(A_->N(), A_->M(), A_->nonzeroes(), Matrix::row_wise));
auto& newA = *ILU;
// Create sparsity pattern
for(auto iter=newA.createbegin(), iend = newA.createend(); iter != iend; ++iter)
{
const auto& row = A[inverseOrdering[iter.index()]];
const auto& row = (*A_)[inverseOrdering[iter.index()]];
for(auto col = row.begin(), cend = row.end(); col != cend; ++col)
{
iter.insert(ordering_[col.index()]);
}
}
// Copy values
for(auto iter = A.begin(), iend = A.end(); iter != iend; ++iter)
for(auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter)
{
auto& newRow = newA[ordering_[iter.index()]];
for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
@ -989,7 +1004,7 @@ protected:
}
}
switch ( milu )
switch ( milu_ )
{
case MILU_VARIANT::MILU_1:
detail::milu0_decomposition ( *ILU);
@ -1007,7 +1022,7 @@ protected:
detail::IsPositiveFunctor() );
break;
default:
if (interiorSize_ == A.N())
if (interiorSize_ == A_->N())
bilu0_decomposition( *ILU );
else
detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
@ -1016,7 +1031,7 @@ protected:
}
else {
// create ILU-n decomposition
ILU.reset( new Matrix( A.N(), A.M(), Matrix::row_wise) );
ILU.reset( new Matrix( A_->N(), A_->M(), Matrix::row_wise) );
std::unique_ptr<detail::Reorderer> reorderer, inverseReorderer;
if ( ordering_.empty() )
{
@ -1029,7 +1044,7 @@ protected:
inverseReorderer.reset(new detail::RealReorderer(inverseOrdering));
}
milun_decomposition( A, iluIteration, milu, *ILU, *reorderer, *inverseReorderer );
milun_decomposition( *A_, iluIteration_, milu_, *ILU, *reorderer, *inverseReorderer );
}
}
catch (const Dune::MatrixBlockError& error)
@ -1053,6 +1068,7 @@ protected:
detail::convertToCRS( *ILU, lower_, upper_, inv_ );
}
protected:
/// \brief Reorder D if needed and return a reference to it.
Range& reorderD(const Range& d)
{
@ -1123,6 +1139,11 @@ protected:
const field_type w_;
const bool relaxation_;
size_type interiorSize_;
const Matrix* A_;
int iluIteration_;
MILU_VARIANT milu_;
bool redBlack_;
bool reorderSphere_;
};
} // end namespace Opm

33
opm/simulators/linalg/PreconditionerFactory.hpp
View File

@ -156,26 +156,20 @@ private:
using C = Comm;
doAddCreator("ILU0", [](const O& op, const P& prm, const C& comm) {
const double w = prm.get<double>("relaxation");
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqILU<M, V, V>>>(comm, op.getmat(), w);
#else
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqILU0<M, V, V>>>(comm, op.getmat(), w);
#endif
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const C& comm) {
const double w = prm.get<double>("relaxation");
// Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner.
return wrapPreconditioner<Opm::ParallelOverlappingILU0<M, V, V, C>>(
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("ILUn", [](const O& op, const P& prm, const C& comm) {
const int n = prm.get<int>("ilulevel");
const double w = prm.get<double>("relaxation");
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqILU<M, V, V>>>(comm, op.getmat(), n, w);
#else
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqILUn<M, V, V>>>(comm, op.getmat(), n, w);
#endif
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
op.getmat(), comm, n, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("Jac", [](const O& op, const P& prm, const C& comm) {
const int n = prm.get<int>("repeats");
@ -229,24 +223,19 @@ private:
using P = boost::property_tree::ptree;
doAddCreator("ILU0", [](const O& op, const P& prm) {
const double w = prm.get<double>("relaxation");
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
return wrapPreconditioner<SeqILU<M, V, V>>(op.getmat(), w);
#else
return wrapPreconditioner<SeqILU0<M, V, V>>(op.getmat(), w);
#endif
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
op.getmat(), 0, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("ParOverILU0", [](const O& op, const P& prm) {
const double w = prm.get<double>("relaxation");
return wrapPreconditioner<Opm::ParallelOverlappingILU0<M, V, V>>(op.getmat(), 0, w, Opm::MILU_VARIANT::ILU);
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
op.getmat(), 0, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("ILUn", [](const O& op, const P& prm) {
const int n = prm.get<int>("ilulevel");
const double w = prm.get<double>("relaxation");
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
return wrapPreconditioner<SeqILU<M, V, V>>(op.getmat(), n, w);
#else
return wrapPreconditioner<SeqILUn<M, V, V>>(op.getmat(), n, w);
#endif
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
op.getmat(), n, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("Jac", [](const O& op, const P& prm) {
const int n = prm.get<int>("repeats");

8
opm/simulators/linalg/PressureTransferPolicy.hpp
View File

@ -85,12 +85,12 @@ public:
assert(entry.index() == entryCoarse.index());
double matrix_el = 0;
if (transpose) {
auto bw = weights_[entry.index()];
const auto& bw = weights_[entry.index()];
for (size_t i = 0; i < bw.size(); ++i) {
matrix_el += (*entry)[pressure_var_index_][i] * bw[i];
}
} else {
auto bw = weights_[row.index()];
const auto& bw = weights_[row.index()];
for (size_t i = 0; i < bw.size(); ++i) {
matrix_el += (*entry)[i][pressure_var_index_] * bw[i];
}
@ -109,7 +109,7 @@ public:
auto end = fine.end(), begin = fine.begin();
for (auto block = begin; block != end; ++block) {
auto bw = weights_[block.index()];
const auto& bw = weights_[block.index()];
double rhs_el = 0.0;
if (transpose) {
rhs_el = (*block)[pressure_var_index_];
@ -130,7 +130,7 @@ public:
for (auto block = begin; block != end; ++block) {
if (transpose) {
auto bw = weights_[block.index()];
const auto& bw = weights_[block.index()];
for (size_t i = 0; i < block->size(); ++i) {
(*block)[i] = this->lhs_[block - begin] * bw[i];
}