mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-07-07 04:53:03 -05:00
301 lines
9.0 KiB
C++
301 lines
9.0 KiB
C++
/*
|
|
Copyright 2015, 2022 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/>.
|
|
*/
|
|
|
|
#include <config.h>
|
|
#include <opm/simulators/linalg/MILU.hpp>
|
|
#include <opm/common/TimingMacros.hpp>
|
|
#include <dune/common/version.hh>
|
|
#include <dune/common/fmatrix.hh>
|
|
#include <dune/istl/bcrsmatrix.hh>
|
|
#include <dune/istl/ilu.hh>
|
|
|
|
#include <opm/common/ErrorMacros.hpp>
|
|
|
|
#include <opm/simulators/linalg/matrixblock.hh>
|
|
|
|
#include <array>
|
|
|
|
template <typename T>
|
|
T Opm::detail::identityFunctor(const T& t){return t;};
|
|
|
|
template <typename T>
|
|
T Opm::detail::oneFunctor(const T&){return 1.0;};
|
|
|
|
template <typename T>
|
|
T Opm::detail::signFunctor(const T& t){if (t< 0){return -1;} else{return 1;}};
|
|
|
|
template <typename T>
|
|
T Opm::detail::isPositiveFunctor(const T& t){if (t<0){return 0;} else{return 1;}};
|
|
|
|
template <typename T>
|
|
T Opm::detail::absFunctor(const T& t){return std::abs(t);};
|
|
|
|
namespace Opm
|
|
{
|
|
|
|
MILU_VARIANT convertString2Milu(const 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;
|
|
}
|
|
|
|
namespace detail
|
|
{
|
|
|
|
template<class M>
|
|
void milu0_decomposition(M& A, FieldFunct<M> absFunctor, FieldFunct<M> signFunctor,
|
|
std::vector<typename M::block_type>* diagonal)
|
|
{
|
|
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 " + std::to_string(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 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, identityFunctor<typename M::field_type>,
|
|
signFunctor<typename M::field_type>);
|
|
|
|
break;
|
|
case MILU_VARIANT::MILU_3:
|
|
detail::milu0_decomposition ( ILU, absFunctor<typename M::field_type>,
|
|
signFunctor<typename M::field_type>);
|
|
break;
|
|
case MILU_VARIANT::MILU_4:
|
|
detail::milu0_decomposition ( ILU, identityFunctor<typename M::field_type> ,
|
|
isPositiveFunctor<typename M::field_type>);
|
|
break;
|
|
default:
|
|
#if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
|
|
bilu0_decomposition( ILU );
|
|
#else
|
|
Dune::ILU::blockILU0Decomposition( ILU );
|
|
#endif
|
|
break;
|
|
}
|
|
}
|
|
|
|
|
|
template double Opm::detail::identityFunctor(const double&);
|
|
template double Opm::detail::oneFunctor(const double&);
|
|
template double Opm::detail::signFunctor(const double&);
|
|
template double Opm::detail::isPositiveFunctor(const double&);
|
|
template double Opm::detail::absFunctor(const double&);
|
|
|
|
#define INSTANCE(...) \
|
|
template void milu0_decomposition<__VA_ARGS__> \
|
|
(__VA_ARGS__&,std::function<double(const double&)>, std::function<double(const double&)>, \
|
|
std::vector<typename __VA_ARGS__::block_type>*);
|
|
|
|
#define INSTANCE_ILUN(...) \
|
|
template void milun_decomposition(const __VA_ARGS__&, int, MILU_VARIANT, \
|
|
__VA_ARGS__&,Reorderer&,Reorderer&);
|
|
|
|
#define INSTANCE_FULL(...) \
|
|
INSTANCE(__VA_ARGS__) \
|
|
INSTANCE_ILUN(__VA_ARGS__)
|
|
|
|
#define INSTANCE_BLOCK(Dim) \
|
|
INSTANCE_FULL(Dune::BCRSMatrix<MatrixBlock<double,Dim,Dim>>)
|
|
|
|
#define INSTANCE_FM(Dim) \
|
|
INSTANCE_FULL(Dune::BCRSMatrix<Dune::FieldMatrix<double,Dim,Dim>>)
|
|
|
|
INSTANCE_FM(1)
|
|
INSTANCE_FM(2)
|
|
INSTANCE_FM(3)
|
|
INSTANCE_FM(4)
|
|
INSTANCE_FM(5)
|
|
INSTANCE_FM(6)
|
|
|
|
INSTANCE_BLOCK(1)
|
|
INSTANCE_BLOCK(2)
|
|
INSTANCE_BLOCK(3)
|
|
INSTANCE_BLOCK(4)
|
|
INSTANCE_BLOCK(5)
|
|
INSTANCE_BLOCK(6)
|
|
|
|
} // end namespace detail
|
|
|
|
} // end namespace Opm
|