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 #4006 from akva2/milu_refactor

Browse Source 
changed: put MILU in separate compile unit
This commit is contained in:
Markus Blatt 2022-08-16 15:31:23 +02:00 committed by GitHub
commit acff71c3d1
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 440 additions and 312 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -51,6 +51,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/linalg/FlexibleSolver4.cpp
opm/simulators/linalg/FlexibleSolver5.cpp
opm/simulators/linalg/FlexibleSolver6.cpp
opm/simulators/linalg/MILU.cpp
opm/simulators/linalg/PropertyTree.cpp
opm/simulators/linalg/setupPropertyTree.cpp
opm/simulators/utils/PartiallySupportedFlowKeywords.cpp

6
opm/simulators/linalg/FlowLinearSolverParameters.hpp
View File

@ -24,15 +24,11 @@
#ifndef OPM_FLOWLINEARSOLVERPARAMETERS_HEADER_INCLUDED
#define OPM_FLOWLINEARSOLVERPARAMETERS_HEADER_INCLUDED
#include <opm/common/utility/parameters/ParameterGroup.hpp>
#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
#include <opm/simulators/linalg/MILU.hpp>
#include <opm/simulators/linalg/linalgproperties.hh>
#include <opm/models/utils/parametersystem.hh>
#include <array>
#include <memory>
namespace Opm {
template <class TypeTag>
class ISTLSolverEbos;

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

@ -26,8 +26,10 @@
#include <opm/models/utils/propertysystem.hh>
#include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
#include <opm/simulators/linalg/FlexibleSolver.hpp>
#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
#include <opm/simulators/linalg/MatrixBlock.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
#include <opm/simulators/linalg/WellOperators.hpp>
#include <opm/simulators/linalg/WriteSystemMatrixHelper.hpp>
#include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>

277
opm/simulators/linalg/MILU.cpp Normal file
View File

@ -0,0 +1,277 @@
/*
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/>.
*/
#include <config.h>
#include <opm/simulators/linalg/MILU.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>
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, class F1, class F2>
void milu0_decomposition(M& A, F1 absFunctor, F2 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 " << 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, 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;
}
}
#define INSTANCE(F1,F2,...) \
template void milu0_decomposition<__VA_ARGS__, F1, F2> \
(__VA_ARGS__&, F1, F2, \
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(AbsFunctor,SignFunctor,__VA_ARGS__) \
INSTANCE(IdentityFunctor,IsPositiveFunctor,__VA_ARGS__) \
INSTANCE(IdentityFunctor,OneFunctor,__VA_ARGS__) \
INSTANCE(IdentityFunctor,SignFunctor,__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_BLOCK(1)
INSTANCE_BLOCK(2)
INSTANCE_BLOCK(3)
INSTANCE_BLOCK(4)
INSTANCE_BLOCK(5)
INSTANCE_BLOCK(6)
} // end namespace detail
} // end namespace Opm

152
opm/simulators/linalg/MILU.hpp Normal file
View File

@ -0,0 +1,152 @@
/*
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_MILU_HEADER_INCLUDED
#define OPM_MILU_HEADER_INCLUDED
#include <cmath>
#include <cstddef>
#include <string>
#include <vector>
namespace Opm
{
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
};
MILU_VARIANT convertString2Milu(const std::string& milu);
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)
{
return std::abs(t);
}
};
template<class M, class F1=IdentityFunctor, class F2=OneFunctor >
void milu0_decomposition(M& A, F1 absFunctor = F1(), F2 signFunctor = F2(),
std::vector<typename M::block_type>* diagonal = nullptr);
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);
} // end namespace details
} // end namespace Opm
#endif

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)

1
opm/simulators/linalg/WellOperators.hpp
View File

@ -24,7 +24,6 @@
#include <dune/istl/operators.hh>
namespace Opm
{

4
opm/simulators/linalg/setupPropertyTree.cpp
View File

@ -21,6 +21,10 @@
#include <opm/simulators/linalg/setupPropertyTree.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
#include <filesystem>
#include <boost/version.hpp>

3
opm/simulators/linalg/setupPropertyTree.hpp
View File

@ -20,7 +20,6 @@
#ifndef OPM_SETUPPROPERTYTREE_HEADER_INCLUDED
#define OPM_SETUPPROPERTYTREE_HEADER_INCLUDED
#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
#include <opm/simulators/linalg/PropertyTree.hpp>
#include <string>
@ -28,6 +27,8 @@
namespace Opm
{
class FlowLinearSolverParameters;
PropertyTree setupPropertyTree(FlowLinearSolverParameters p,
bool LinearSolverMaxIterSet,
bool CprMaxEllIterSet);