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

changed: put MILU in separate compile unit

Browse Source
This commit is contained in:
Arne Morten Kvarving
2022-08-16 10:45:32 +02:00
parent 874e4f818c
commit a32a32649d
9 changed files with 440 additions and 312 deletions
Download Patch File Download Diff File Expand all files Collapse all files

306
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/MILU.hpp>
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <dune/common/version.hh>
@@ -44,36 +45,6 @@ namespace Opm
template<class Matrix, class Domain, class Range, class ParallelInfo = Dune::Amg::SequentialInformation>
class ParallelOverlappingILU0;
enum class MILU_VARIANT{
/// \brief Do not perform modified ILU
ILU = 0,
/// \brief \f$U_{ii} = U_{ii} +\f$ sum(dropped entries)
MILU_1 = 1,
/// \brief \f$U_{ii} = U_{ii} + sign(U_{ii}) * \f$ sum(dropped entries)
MILU_2 = 2,
/// \brief \f$U_{ii} = U_{ii} sign(U_{ii}) * \f$ sum(|dropped entries|)
MILU_3 = 3,
/// \brief \f$U_{ii} = U_{ii} + (U_{ii}>0?1:0) * \f$ sum(dropped entries)
MILU_4 = 4
};
inline MILU_VARIANT convertString2Milu(std::string milu)
{
if( 0 == milu.compare("MILU_1") )
{
return MILU_VARIANT::MILU_1;
}
if ( 0 == milu.compare("MILU_2") )
{
return MILU_VARIANT::MILU_2;
}
if ( 0 == milu.compare("MILU_3") )
{
return MILU_VARIANT::MILU_3;
}
return MILU_VARIANT::ILU;
}
template<class F>
class ParallelOverlappingILU0Args
: public Dune::Amg::DefaultSmootherArgs<F>
@@ -165,281 +136,6 @@ namespace Opm
{
namespace detail
{
struct Reorderer
{
virtual std::size_t operator[](std::size_t i) const = 0;
virtual ~Reorderer() {}
};
struct NoReorderer : public Reorderer
{
virtual std::size_t operator[](std::size_t i) const
{
return i;
}
};
struct RealReorderer : public Reorderer
{
RealReorderer(const std::vector<std::size_t>& ordering)
: ordering_(&ordering)
{}
virtual std::size_t operator[](std::size_t i) const
{
return (*ordering_)[i];
}
const std::vector<std::size_t>* ordering_;
};
struct IdentityFunctor
{
template<class T>
T operator()(const T& t)
{
return t;
}
};
struct OneFunctor
{
template<class T>
T operator()(const T&)
{
return 1.0;
}
};
struct SignFunctor
{
template<class T>
double operator()(const T& t)
{
if ( t < 0.0 )
{
return -1;
}
else
{
return 1.0;
}
}
};
struct IsPositiveFunctor
{
template<class T>
double operator()(const T& t)
{
if ( t < 0.0 )
{
return 0;
}
else
{
return 1;
}
}
};
struct AbsFunctor
{
template<class T>
T operator()(const T& t)
{
using std::abs;
return abs(t);
}
};
template<class M, class F1=detail::IdentityFunctor, class F2=detail::OneFunctor >
void milu0_decomposition(M& A, F1 absFunctor = F1(), F2 signFunctor = F2(),
std::vector<typename M::block_type>* diagonal = nullptr)
{
if( diagonal )
{
diagonal->reserve(A.N());
}
for ( auto irow = A.begin(), iend = A.end(); irow != iend; ++irow)
{
auto a_i_end = irow->end();
auto a_ik = irow->begin();
std::array<typename M::field_type, M::block_type::rows> sum_dropped{};
// Eliminate entries in lower triangular matrix
// and store factors for L
for ( ; a_ik.index() < irow.index(); ++a_ik )
{
auto k = a_ik.index();
auto a_kk = A[k].find(k);
// L_ik = A_kk^-1 * A_ik
a_ik->rightmultiply(*a_kk);
// modify the rest of the row, everything right of a_ik
// a_i* -=a_ik * a_k*
auto a_k_end = A[k].end();
auto a_kj = a_kk, a_ij = a_ik;
++a_kj; ++a_ij;
while ( a_kj != a_k_end)
{
auto modifier = *a_kj;
modifier.leftmultiply(*a_ik);
while( a_ij != a_i_end && a_ij.index() < a_kj.index())
{
++a_ij;
}
if ( a_ij != a_i_end && a_ij.index() == a_kj.index() )
{
// Value is not dropped
*a_ij -= modifier;
++a_ij; ++a_kj;
}
else
{
auto entry = sum_dropped.begin();
for( const auto& row: modifier )
{
for( const auto& colEntry: row )
{
*entry += absFunctor(-colEntry);
}
++entry;
}
++a_kj;
}
}
}
if ( a_ik.index() != irow.index() )
OPM_THROW(std::logic_error, "Matrix is missing diagonal for row " << irow.index());
int index = 0;
for(const auto& entry: sum_dropped)
{
auto& bdiag = (*a_ik)[index][index];
bdiag += signFunctor(bdiag) * entry;
++index;
}
if ( diagonal )
{
diagonal->push_back(*a_ik);
}
a_ik->invert(); // compute inverse of diagonal block
}
}
template<class M>
void milu0_decomposition(M& A,
std::vector<typename M::block_type>* diagonal)
{
milu0_decomposition(A, detail::IdentityFunctor(), detail::OneFunctor(),
diagonal);
}
template<class M>
void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU,
Reorderer& ordering, Reorderer& inverseOrdering)
{
using Map = std::map<std::size_t, int>;
auto iluRow = ILU.createbegin();
for(std::size_t i = 0, iend = A.N(); i < iend; ++i)
{
auto& orow = A[inverseOrdering[i]];
Map rowPattern;
for ( auto col = orow.begin(), cend = orow.end(); col != cend; ++col)
{
rowPattern[ordering[col.index()]] = 0;
}
for(auto ik = rowPattern.begin(); ik->first < i; ++ik)
{
if ( ik->second < n )
{
auto& rowk = ILU[ik->first];
for ( auto kj = rowk.find(ik->first), endk = rowk.end();
kj != endk; ++kj)
{
// Assume double and block_type FieldMatrix
// first element is misused to store generation number
int generation = (*kj)[0][0];
if(generation < n)
{
auto ij = rowPattern.find(kj.index());
if ( ij == rowPattern.end() )
{
rowPattern[ordering[kj.index()]] = generation + 1;
}
}
}
}
}
// create the row
for (const auto& entry : rowPattern)
{
iluRow.insert(entry.first);
}
++iluRow;
// write generation to newly created row.
auto generationPair = rowPattern.begin();
for ( auto col = ILU[i].begin(), cend = ILU[i].end(); col != cend;
++col, ++generationPair)
{
assert(col.index() == generationPair->first);
(*col)[0][0] = generationPair->second;
}
}
// copy Entries from A
for(auto iter=A.begin(), iend = A.end(); iter != iend; ++iter)
{
auto& newRow = ILU[ordering[iter.index()]];
// reset stored generation
for ( auto& col: newRow)
{
col = 0;
}
// copy row.
for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
{
newRow[ordering[col.index()]] = *col;
}
}
// call decomposition on pattern
switch ( milu )
{
case MILU_VARIANT::MILU_1:
detail::milu0_decomposition ( ILU);
break;
case MILU_VARIANT::MILU_2:
detail::milu0_decomposition ( ILU, detail::IdentityFunctor(),
detail::SignFunctor() );
break;
case MILU_VARIANT::MILU_3:
detail::milu0_decomposition ( ILU, detail::AbsFunctor(),
detail::SignFunctor() );
break;
case MILU_VARIANT::MILU_4:
detail::milu0_decomposition ( ILU, detail::IdentityFunctor(),
detail::IsPositiveFunctor() );
break;
default:
#if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
bilu0_decomposition( ILU );
#else
Dune::ILU::blockILU0Decomposition( ILU );
#endif
break;
}
}
//! Compute Blocked ILU0 decomposition, when we know junk ghost rows are located at the end of A
template<class M>
void ghost_last_bilu0_decomposition (M& A, size_t interiorSize)