From a32a32649dcca30b73bd4246fa31c638f5be4339 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 16 Aug 2022 10:45:32 +0200 Subject: [PATCH] changed: put MILU in separate compile unit --- CMakeLists_files.cmake | 1 + .../linalg/FlowLinearSolverParameters.hpp | 6 +- opm/simulators/linalg/ISTLSolverEbos.hpp | 2 + opm/simulators/linalg/MILU.cpp | 277 ++++++++++++++++ opm/simulators/linalg/MILU.hpp | 152 +++++++++ .../linalg/ParallelOverlappingILU0.hpp | 306 +----------------- opm/simulators/linalg/WellOperators.hpp | 1 - opm/simulators/linalg/setupPropertyTree.cpp | 4 + opm/simulators/linalg/setupPropertyTree.hpp | 3 +- 9 files changed, 440 insertions(+), 312 deletions(-) create mode 100644 opm/simulators/linalg/MILU.cpp create mode 100644 opm/simulators/linalg/MILU.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index e2bd95250..fc364b370 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -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 diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.hpp b/opm/simulators/linalg/FlowLinearSolverParameters.hpp index 4fdfa6df5..8a20b732f 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.hpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.hpp @@ -24,15 +24,11 @@ #ifndef OPM_FLOWLINEARSOLVERPARAMETERS_HEADER_INCLUDED #define OPM_FLOWLINEARSOLVERPARAMETERS_HEADER_INCLUDED -#include -#include +#include #include #include -#include -#include - namespace Opm { template class ISTLSolverEbos; diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 07248cce9..3b8acc9a0 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -26,8 +26,10 @@ #include #include #include +#include #include #include +#include #include #include #include diff --git a/opm/simulators/linalg/MILU.cpp b/opm/simulators/linalg/MILU.cpp new file mode 100644 index 000000000..2a6269ce6 --- /dev/null +++ b/opm/simulators/linalg/MILU.cpp @@ -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 . +*/ + +#include +#include + +#include +#include +#include +#include + +#include + +#include + +#include + +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 +void milu0_decomposition(M& A, F1 absFunctor, F2 signFunctor, + std::vector* 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 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 +void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU, + Reorderer& ordering, Reorderer& inverseOrdering) +{ + using Map = std::map; + + 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*); + +#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>) + +#define INSTANCE_FM(Dim) \ + INSTANCE_FULL(Dune::BCRSMatrix>) + +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 diff --git a/opm/simulators/linalg/MILU.hpp b/opm/simulators/linalg/MILU.hpp new file mode 100644 index 000000000..6491986d8 --- /dev/null +++ b/opm/simulators/linalg/MILU.hpp @@ -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 . +*/ +#ifndef OPM_MILU_HEADER_INCLUDED +#define OPM_MILU_HEADER_INCLUDED + +#include +#include +#include +#include + +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& ordering) + : ordering_(&ordering) + {} + virtual std::size_t operator[](std::size_t i) const + { + return (*ordering_)[i]; + } + const std::vector* ordering_; +}; + +struct IdentityFunctor +{ + template + T operator()(const T& t) + { + return t; + } +}; + +struct OneFunctor +{ + template + T operator()(const T&) + { + return 1.0; + } +}; +struct SignFunctor +{ + template + double operator()(const T& t) + { + if (t < 0.0) + { + return -1; + } + else + { + return 1.0; + } + } +}; + +struct IsPositiveFunctor +{ + template + double operator()(const T& t) + { + if (t < 0.0) + { + return 0; + } + else + { + return 1; + } + } +}; +struct AbsFunctor +{ + template + T operator()(const T& t) + { + return std::abs(t); + } +}; + +template +void milu0_decomposition(M& A, F1 absFunctor = F1(), F2 signFunctor = F2(), + std::vector* diagonal = nullptr); + +template +void milu0_decomposition(M& A, + std::vector* diagonal) +{ + milu0_decomposition(A, detail::IdentityFunctor(), detail::OneFunctor(), + diagonal); +} + +template +void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU, + Reorderer& ordering, Reorderer& inverseOrdering); + +} // end namespace details + +} // end namespace Opm + +#endif diff --git a/opm/simulators/linalg/ParallelOverlappingILU0.hpp b/opm/simulators/linalg/ParallelOverlappingILU0.hpp index 53e4290aa..0dfde2151 100644 --- a/opm/simulators/linalg/ParallelOverlappingILU0.hpp +++ b/opm/simulators/linalg/ParallelOverlappingILU0.hpp @@ -21,6 +21,7 @@ #define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED #include +#include #include #include #include @@ -44,36 +45,6 @@ namespace Opm template 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 ParallelOverlappingILU0Args : public Dune::Amg::DefaultSmootherArgs @@ -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& ordering) - : ordering_(&ordering) - {} - virtual std::size_t operator[](std::size_t i) const - { - return (*ordering_)[i]; - } - const std::vector* ordering_; - }; - - struct IdentityFunctor - { - template - T operator()(const T& t) - { - return t; - } - }; - - struct OneFunctor - { - template - T operator()(const T&) - { - return 1.0; - } - }; - struct SignFunctor - { - template - double operator()(const T& t) - { - if ( t < 0.0 ) - { - return -1; - } - else - { - return 1.0; - } - } - }; - - struct IsPositiveFunctor - { - template - double operator()(const T& t) - { - if ( t < 0.0 ) - { - return 0; - } - else - { - return 1; - } - } - }; - struct AbsFunctor - { - template - T operator()(const T& t) - { - using std::abs; - return abs(t); - } - }; - - template - void milu0_decomposition(M& A, F1 absFunctor = F1(), F2 signFunctor = F2(), - std::vector* 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 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 - void milu0_decomposition(M& A, - std::vector* diagonal) - { - milu0_decomposition(A, detail::IdentityFunctor(), detail::OneFunctor(), - diagonal); - } - - template - void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU, - Reorderer& ordering, Reorderer& inverseOrdering) - { - using Map = std::map; - - 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 void ghost_last_bilu0_decomposition (M& A, size_t interiorSize) diff --git a/opm/simulators/linalg/WellOperators.hpp b/opm/simulators/linalg/WellOperators.hpp index 06afa260a..692442c30 100644 --- a/opm/simulators/linalg/WellOperators.hpp +++ b/opm/simulators/linalg/WellOperators.hpp @@ -24,7 +24,6 @@ #include - namespace Opm { diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 6b5f1d077..764eb35d6 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -21,6 +21,10 @@ #include +#include + +#include + #include #include diff --git a/opm/simulators/linalg/setupPropertyTree.hpp b/opm/simulators/linalg/setupPropertyTree.hpp index 9bfd5b7a7..6c92db417 100644 --- a/opm/simulators/linalg/setupPropertyTree.hpp +++ b/opm/simulators/linalg/setupPropertyTree.hpp @@ -20,7 +20,6 @@ #ifndef OPM_SETUPPROPERTYTREE_HEADER_INCLUDED #define OPM_SETUPPROPERTYTREE_HEADER_INCLUDED -#include #include #include @@ -28,6 +27,8 @@ namespace Opm { +class FlowLinearSolverParameters; + PropertyTree setupPropertyTree(FlowLinearSolverParameters p, bool LinearSolverMaxIterSet, bool CprMaxEllIterSet);