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 #1495 from blattms/red-black-ilu

Browse Source 
Added support for red black ilu
This commit is contained in:
Bård Skaflestad 2018-08-01 18:52:07 +02:00 committed by GitHub
commit b087cf54af
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 1199 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists_files.cmake
View File

@ -142,6 +142,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_blackoil_amg.cpp
tests/test_block.cpp
tests/test_boprops_ad.cpp
tests/test_graphcoloring.cpp
tests/test_rateconverter.cpp
tests/test_span.cpp
tests/test_syntax.cpp
@ -150,6 +151,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_welldensitysegmented.cpp
tests/test_vfpproperties.cpp
tests/test_singlecellsolves.cpp
tests/test_milu.cpp
tests/test_multmatrixtransposed.cpp
tests/test_multiphaseupwind.cpp
tests/test_wellmodel.cpp
@ -290,6 +292,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/FlowMainEbos.hpp
opm/autodiff/FlowMainSequential.hpp
opm/autodiff/GeoProps.hpp
opm/autodiff/GraphColoring.hpp
opm/autodiff/GridHelpers.hpp
opm/autodiff/GridInit.hpp
opm/autodiff/ImpesTPFAAD.hpp

96
opm/autodiff/CPRPreconditioner.hpp
View File

@ -29,6 +29,7 @@
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <opm/common/utility/parameters/ParameterGroup.hpp>
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
@ -64,6 +65,9 @@ struct CPRParameter
double cpr_relax_;
double cpr_solver_tol_;
int cpr_ilu_n_;
MILU_VARIANT cpr_ilu_milu_;
bool cpr_ilu_redblack_;
bool cpr_ilu_reorder_sphere_;
int cpr_max_ell_iter_;
bool cpr_use_amg_;
bool cpr_use_bicgstab_;
@ -77,25 +81,33 @@ struct CPRParameter
// reset values to default
reset();
cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_);
cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_);
cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_);
cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_);
cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_);
cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_);
cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_);
cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_);
cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_);
cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_);
cpr_ilu_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_);
cpr_ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_);
cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_);
cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_);
cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_);
cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_);
cpr_pressure_aggregation_ = param.getDefault("cpr_pressure_aggregation", cpr_pressure_aggregation_);
std::string milu("ILU");
cpr_ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu));
}
void reset()
{
cpr_relax_ = 1.0;
cpr_solver_tol_ = 1e-2;
cpr_ilu_n_ = 0;
cpr_max_ell_iter_ = 25;
cpr_use_amg_ = true;
cpr_use_bicgstab_ = true;
cpr_solver_verbose_ = false;
cpr_relax_ = 1.0;
cpr_solver_tol_ = 1e-2;
cpr_ilu_n_ = 0;
cpr_ilu_milu_ = MILU_VARIANT::ILU;
cpr_ilu_redblack_ = false;
cpr_ilu_reorder_sphere_ = true;
cpr_max_ell_iter_ = 25;
cpr_use_amg_ = true;
cpr_use_bicgstab_ = true;
cpr_solver_verbose_ = false;
cpr_pressure_aggregation_ = false;
}
};
@ -103,6 +115,31 @@ struct CPRParameter
namespace ISTLUtility
{
template<class T>
void setILUParameters(Opm::ParallelOverlappingILU0Args<T>& args,
const CPRParameter& params)
{
args.setN(params.cpr_ilu_n_);
args.setMilu(params.cpr_ilu_milu_);
}
template<class T>
void setILUParameters(Opm::ParallelOverlappingILU0Args<T>& args,
MILU_VARIANT milu, int n=0)
{
args.setN(n);
args.setMilu(milu);
}
template<class S, class P>
void setILUParameters(S&, const P&)
{}
template<class S, class P>
void setILUParameters(S&, bool, int)
{}
///
/// \brief A traits class for selecting the types of the preconditioner.
///
@ -167,9 +204,9 @@ struct CPRSelector<M,X,Y,Dune::Amg::SequentialInformation>
//! \param relax The relaxation factor to use.
template<class M, class X, class C>
std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
createILU0Ptr(const M& A, const C& comm, double relax)
createILU0Ptr(const M& A, const C& comm, double relax, MILU_VARIANT milu)
{
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >(A, comm, relax);
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >(A, comm, relax, milu);
}
//! \brief Creates and initializes a shared pointer to an ILUn preconditioner.
//! \param A The matrix of the linear system to solve.
@ -177,25 +214,28 @@ createILU0Ptr(const M& A, const C& comm, double relax)
//! \param relax The relaxation factor to use.
template<class M, class X, class C>
std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
createILUnPtr(const M& A, const C& comm, int ilu_n, double relax)
createILUnPtr(const M& A, const C& comm, int ilu_n, double relax, MILU_VARIANT milu)
{
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >( A, comm, ilu_n, relax);
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >( A, comm, ilu_n, relax, milu );
}
/// \brief Creates the elliptic preconditioner (ILU0)
/// \param Ae The matrix of the elliptic system.
/// \param relax The relaxation parameter for ILU0.
/// \param milu If true, the modified ilu approach is used. Dropped elements
/// will get added to the diagonal of U to preserve the row sum
/// for constant vectors (Ae = LUe).
/// \param comm The object describing the parallelization information and communication.
template<class M, class X=typename M::range_type, class P>
typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
createEllipticPreconditionerPointer(const M& Ae, double relax,
const P& comm)
MILU_VARIANT milu, const P& comm)
{
typedef typename CPRSelector<M,X,X,P >
::EllipticPreconditioner ParallelPreconditioner;
typedef typename CPRSelector<M,X,X,P>
::EllipticPreconditionerPointer EllipticPreconditionerPointer;
return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax));
return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax, milu));
}
template < class C, class Op, class P, class S, std::size_t index >
@ -220,13 +260,14 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm,
SmootherArgs smootherArgs;
smootherArgs.iterations = 1;
smootherArgs.relaxationFactor = relax;
setILUParameters(smootherArgs, params);
amgPtr.reset( new AMG( params, opA, criterion, smootherArgs, comm ) );
}
template < class C, class Op, class P, class AMG >
inline void
createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::unique_ptr< AMG >& amgPtr)
createAMGPreconditionerPointer(Op& opA, const double relax, const MILU_VARIANT milu, const P& comm, std::unique_ptr< AMG >& amgPtr)
{
// TODO: revise choice of parameters
int coarsenTarget=1200;
@ -243,6 +284,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::
SmootherArgs smootherArgs;
smootherArgs.iterations = 1;
smootherArgs.relaxationFactor = relax;
setILUParameters(smootherArgs, milu);
amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) );
}
@ -254,7 +296,7 @@ createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm, std::
// \param amgPtr The unique_ptr to be filled (return)
template < int pressureIndex=0, class Op, class P, class AMG >
inline void
createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std::unique_ptr< AMG >& amgPtr )
createAMGPreconditionerPointer( Op& opA, const double relax, const MILU_VARIANT milu, const P& comm, std::unique_ptr< AMG >& amgPtr )
{
// type of matrix
typedef typename Op::matrix_type M;
@ -268,7 +310,7 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
// The coarsening criterion used in the AMG
typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
createAMGPreconditionerPointer<Criterion>(opA, relax, comm, amgPtr);
createAMGPreconditionerPointer<Criterion>(opA, relax, milu, comm, amgPtr);
}
} // end namespace ISTLUtility
@ -370,10 +412,10 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
// create the preconditioner for the whole system.
if( param_.cpr_ilu_n_ == 0 ) {
pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, comm, param_.cpr_relax_ );
pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, comm, param_.cpr_relax_, param_.cpr_ilu_milu_ );
}
else {
pre_ = ISTLUtility::createILUnPtr<M,X>( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_);
pre_ = ISTLUtility::createILUnPtr<M,X>( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_, param_.cpr_ilu_milu_);
}
}
@ -531,11 +573,11 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
{
if( amg )
{
ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, comm, amg_ );
ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, param_.cpr_ilu_milu_, comm, amg_ );
}
else
{
precond_ = ISTLUtility::createEllipticPreconditionerPointer<M,X>( Ae_, param_.cpr_relax_, comm);
precond_ = ISTLUtility::createEllipticPreconditionerPointer<M,X>( Ae_, param_.cpr_relax_, param_.cpr_ilu_milu_, comm);
}
}
};

245
opm/autodiff/GraphColoring.hpp Normal file
View File

@ -0,0 +1,245 @@
/*
Copyright 2018 Equinor
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_GRAPHCOLORING_HEADER_INCLUDED
#define OPM_GRAPHCOLORING_HEADER_INCLUDED
#include <vector>
#include <deque>
#include <tuple>
#include <algorithm>
#include <numeric>
#include <queue>
#include <cstddef>
namespace Opm
{
namespace Detail
{
template<class Graph>
std::size_t colorGraphWelshPowell(const Graph& graph,
std::deque<typename Graph::VertexDescriptor>& orderedVertices,
std::vector<int>& colors,
int color, int noVertices)
{
std::vector<int> forbidden(noVertices, false);
std::size_t noColored = 0;
for(auto vertex = orderedVertices.begin(),
vertexEnd = orderedVertices.end();
vertex != vertexEnd; ++vertex)
{
// Skip forbidden vertices
while(vertex != vertexEnd && forbidden[*vertex])
++vertex;
if ( vertex == vertexEnd )
{
break;
}
// Color Vertex
colors[*vertex] = color;
++noColored;
// Forbid neighors
for(auto edge = graph.beginEdges(*vertex), endEdge = graph.endEdges(*vertex);
edge != endEdge; ++edge)
{
forbidden[edge.target()] = true;
}
}
// forbidden vertices will be colored next for coloring
using Vertex = typename Graph::VertexDescriptor;
auto newEnd = std::remove_if(orderedVertices.begin(), orderedVertices.end(),
[&forbidden](const Vertex& vertex)
{
return !forbidden[vertex];
});
orderedVertices.resize(newEnd-orderedVertices.begin());
return noColored;
}
template<class Graph, class Functor>
std::size_t breadthFirstSearch(const Graph& graph, typename Graph::VertexDescriptor root,
Functor functor)
{
std::vector<int> visited(graph.maxVertex() + 1, false);
using Vertex = typename Graph::VertexDescriptor;
std::queue<Vertex> nextVertices;
std::size_t noVisited = 0;
nextVertices.push(root);
visited[root] = true; // We do not visit root.
while( !nextVertices.empty() )
{
auto current = nextVertices.front();
for(auto edge = graph.beginEdges(current),
endEdge = graph.endEdges(current);
edge != endEdge; ++edge)
{
if ( ! visited[edge.target()] )
{
visited[edge.target()] = true;
nextVertices.push(edge.target());
functor(edge.target());
++noVisited;
}
}
nextVertices.pop();
}
return noVisited;
}
} // end namespace Detail
/// \brief Color the vertices of graph.
///
/// It uses the algorithm of Welsh and Powell for this.
/// \param graph The graph to color. Must adhere to the graph interface of dune-istl.
/// \return A pair of a vector with the colors of the vertices and the number of colors
/// assigned
template<class Graph>
std::tuple<std::vector<int>, int, std::vector<std::size_t> >
colorVerticesWelshPowell(const Graph& graph)
{
using Vertex = typename Graph::VertexDescriptor;
std::deque<Vertex> orderedVertices;
auto noVertices = graph.maxVertex()+1;
std::vector<int> degrees(noVertices, 0);
int maxDegree = 0;
std::ptrdiff_t firstDegreeChange = 0;
// populate deque
for( auto vertex = graph.begin(), endVertex = graph.end();
vertex != endVertex; ++vertex)
{
auto currentVertex = *vertex;
auto& degree = degrees[currentVertex];
for(auto edge = graph.beginEdges(currentVertex),
endEdge = graph.endEdges(currentVertex);
edge != endEdge; ++edge)
{
++degree;
}
if( degree >= maxDegree )
{
orderedVertices.emplace_front(currentVertex);
++firstDegreeChange;
if(degree > maxDegree)
{
firstDegreeChange = 1;
maxDegree = degree;
}
}
else
{
orderedVertices.emplace_back(currentVertex);
}
}
// order deque by descending degree
std::stable_sort(orderedVertices.begin() + firstDegreeChange,
orderedVertices.end(),
[&degrees](const Vertex& v1, const Vertex& v2)
{
return degrees[v1] > degrees[v2];
});
// Overwrite degree with color
auto& colors = degrees;
std::fill(colors.begin(), colors.end(), -1);
int color = 0;
std::vector<std::size_t> verticesPerColor;
verticesPerColor.reserve(10);
while(!orderedVertices.empty())
{
verticesPerColor
.push_back(Detail::colorGraphWelshPowell(graph, orderedVertices, colors,
color++, noVertices));
}
return std::make_tuple(colors, color, verticesPerColor);
}
/// \! Reorder colored graph preserving order of vertices with the same color.
template<class Graph>
std::vector<std::size_t>
reorderVerticesPreserving(const std::vector<int>& colors, int noColors,
const std::vector<std::size_t>& verticesPerColor,
const Graph& graph)
{
std::vector<std::size_t> colorIndex(noColors, 0);
std::vector<std::size_t> indices(graph.maxVertex() + 1);
std::partial_sum(verticesPerColor.begin(),
verticesPerColor.begin()+verticesPerColor.size() - 1,
colorIndex.begin() + 1);
for(const auto& vertex: graph)
{
indices[vertex] = colorIndex[colors[vertex]]++;
}
return indices;
}
/// \! Reorder Vetrices in spheres
template<class Graph>
std::vector<std::size_t>
reorderVerticesSpheres(const std::vector<int>& colors, int noColors,
const std::vector<std::size_t>& verticesPerColor,
const Graph& graph,
typename Graph::VertexDescriptor root)
{
std::vector<std::size_t> colorIndex(noColors, 0);
const auto notVisitedTag = std::numeric_limits<std::size_t>::max();
std::vector<std::size_t> indices(graph.maxVertex() + 1, notVisitedTag);
using Vertex = typename Graph::VertexDescriptor;
std::partial_sum(verticesPerColor.begin(),
verticesPerColor.begin()+verticesPerColor.size() - 1,
colorIndex.begin() + 1);
std::size_t noVisited = 0;
auto numberer = [&colorIndex, &colors, &indices](Vertex vertex)
{
indices[vertex] = colorIndex[colors[vertex]]++;
};
while ( noVisited < graph.maxVertex() + 1 )
{
numberer(root);
++noVisited; //root node already visited and not visited in BFS
noVisited += Detail::breadthFirstSearch(graph, root, numberer);
if ( noVisited < graph.maxVertex() + 1 )
{
// Graph is disconnected search for not yet visited node
for(auto vertex: graph)
{
if ( indices[vertex] == notVisitedTag )
{
// \todo make sure that this is a peripheral node!
root = vertex;
break;
}
}
}
}
return indices;
}
} // end namespace Opm
#endif

36
opm/autodiff/ISTLSolver.hpp
View File

@ -440,6 +440,7 @@ namespace Detail
}
const double relax = parameters_.ilu_relaxation_;
const MILU_VARIANT ilu_milu = parameters_.ilu_milu_;
if ( parameters_.use_cpr_ )
{
using Matrix = typename MatrixOperator::matrix_type;
@ -452,7 +453,7 @@ namespace Detail
std::unique_ptr< AMG > amg;
// Construct preconditioner.
Criterion crit(15, 2000);
constructAMGPrecond<Criterion>( linearOperator, parallelInformation_arg, amg, opA, relax );
constructAMGPrecond<Criterion>( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu );
// Solve.
solve(linearOperator, x, istlb, *sp, *amg, result);
@ -463,7 +464,7 @@ namespace Detail
std::unique_ptr< AMG > amg;
// Construct preconditioner.
constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu );
// Solve.
solve(linearOperator, x, istlb, *sp, *amg, result);
@ -495,7 +496,10 @@ namespace Detail
{
const double relax = parameters_.ilu_relaxation_;
const int ilu_fillin = parameters_.ilu_fillin_level_;
std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax));
const MILU_VARIANT ilu_milu = parameters_.ilu_milu_;
const bool ilu_redblack = parameters_.ilu_redblack_;
const bool ilu_reorder_spheres = parameters_.ilu_reorder_sphere_;
std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax, ilu_milu, ilu_redblack, ilu_reorder_spheres));
return precond;
}
@ -517,38 +521,46 @@ namespace Detail
{
typedef std::unique_ptr<ParPreconditioner> Pointer;
const double relax = parameters_.ilu_relaxation_;
return Pointer(new ParPreconditioner(opA.getmat(), comm, relax));
const MILU_VARIANT ilu_milu = parameters_.ilu_milu_;
const bool ilu_redblack = parameters_.ilu_redblack_;
const bool ilu_reorder_spheres = parameters_.ilu_reorder_sphere_;
return Pointer(new ParPreconditioner(opA.getmat(), comm, relax, ilu_milu, ilu_redblack, ilu_reorder_spheres));
}
#endif
template <class LinearOperator, class MatrixOperator, class POrComm, class AMG >
void
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax, const MILU_VARIANT milu) const
{
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, comm, amg );
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( *opA, relax, milu, comm, amg );
}
template <class MatrixOperator, class POrComm, class AMG >
void
constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const
constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax,
const MILU_VARIANT milu) const
{
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax, comm, amg );
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax,
milu, comm, amg );
}
template <class C, class LinearOperator, class MatrixOperator, class POrComm, class AMG >
void
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax,
const MILU_VARIANT milu ) const
{
ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax, comm, amg, parameters_ );
ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax,
comm, amg, parameters_ );
}
template <class C, class MatrixOperator, class POrComm, class AMG >
void
constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const
constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax, const MILU_VARIANT milu ) const
{
ISTLUtility::template createAMGPreconditionerPointer<C>( opA, relax, comm, amg, parameters_ );
ISTLUtility::template createAMGPreconditionerPointer<C>( opA, relax, milu,
comm, amg, parameters_ );
}
/// \brief Solve the system using the given preconditioner and scalar product.
template <class Operator, class ScalarProd, class Precond>

11
opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
View File

@ -28,6 +28,7 @@
#include <opm/autodiff/CPRPreconditioner.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/common/utility/parameters/ParameterGroup.hpp>
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
#include <array>
#include <memory>
@ -44,6 +45,9 @@ namespace Opm
int linear_solver_restart_;
int linear_solver_verbosity_;
int ilu_fillin_level_;
Opm::MILU_VARIANT ilu_milu_;
bool ilu_redblack_;
bool ilu_reorder_sphere_;
bool newton_use_gmres_;
bool require_full_sparsity_pattern_;
bool ignoreConvergenceFailure_;
@ -69,6 +73,10 @@ namespace Opm
linear_solver_use_amg_ = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ );
ilu_relaxation_ = param.getDefault("ilu_relaxation", ilu_relaxation_ );
ilu_fillin_level_ = param.getDefault("ilu_fillin_level", ilu_fillin_level_ );
ilu_redblack_ = param.getDefault("ilu_redblack", cpr_ilu_redblack_);
ilu_reorder_sphere_ = param.getDefault("ilu_reorder_sphere", cpr_ilu_reorder_sphere_);
std::string milu("ILU");
ilu_milu_ = convertString2Milu(param.getDefault("ilu_milu", milu));
// Check whether to use cpr approach
const std::string cprSolver = "cpr";
@ -89,6 +97,9 @@ namespace Opm
linear_solver_use_amg_ = false;
ilu_fillin_level_ = 0;
ilu_relaxation_ = 0.9;
ilu_milu_ = MILU_VARIANT::ILU;
ilu_redblack_ = false;
ilu_reorder_sphere_ = true;
}
};

580
opm/autodiff/ParallelOverlappingILU0.hpp
View File

@ -20,14 +20,20 @@
#ifndef OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
#define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
#include <opm/autodiff/GraphColoring.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <dune/common/version.hh>
#include <dune/istl/preconditioner.hh>
#include <dune/istl/paamg/smoother.hh>
#include <dune/istl/paamg/graph.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <type_traits>
#include <numeric>
#include <limits>
#include <cstddef>
#include <string>
namespace Opm
{
@ -37,6 +43,64 @@ 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>
{
public:
ParallelOverlappingILU0Args(MILU_VARIANT milu = MILU_VARIANT::ILU )
: milu_(milu)
{}
void setMilu(MILU_VARIANT milu)
{
milu_ = milu;
}
MILU_VARIANT getMilu() const
{
return milu_;
}
void setN(int n)
{
n_ = n;
}
int getN() const
{
return n_;
}
private:
MILU_VARIANT milu_;
int n_;
};
} // end namespace Opm
namespace Dune
@ -45,6 +109,13 @@ namespace Dune
namespace Amg
{
template<class M, class X, class Y, class C>
struct SmootherTraits<Opm::ParallelOverlappingILU0<M,X,Y,C> >
{
using Arguments = Opm::ParallelOverlappingILU0Args<typename M::field_type>;
};
/// \brief Tells AMG how to construct the Opm::ParallelOverlappingILU0 smoother
/// \tparam Matrix The type of the Matrix.
/// \tparam Domain The type of the Vector representing the domain.
@ -54,17 +125,18 @@ namespace Amg
template<class Matrix, class Domain, class Range, class ParallelInfo>
struct ConstructionTraits<Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> >
{
typedef Dune::SeqILU0<Matrix,Domain,Range> T;
typedef Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> T;
typedef DefaultParallelConstructionArgs<T,ParallelInfo> Arguments;
typedef ConstructionTraits<T> SeqConstructionTraits;
static inline Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* construct(Arguments& args)
{
return new Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>(args.getMatrix(),
args.getComm(),
args.getArgs().relaxationFactor);
return new T(args.getMatrix(),
args.getComm(),
args.getArgs().getN(),
args.getArgs().relaxationFactor,
args.getArgs().getMilu());
}
static inline void deconstruct(Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* bp)
static inline void deconstruct(T* bp)
{
delete bp;
}
@ -81,6 +153,278 @@ namespace Opm
{
namespace detail
{
struct Reorderer
{
virtual std::size_t operator[](std::size_t i) const = 0;
};
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];
if(entry<0)
std::cout << entry << std::endl;
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:
bilu0_decomposition( ILU );
break;
}
}
//! compute ILU decomposition of A. A is overwritten by its decomposition
template<class M, class CRS, class InvVector>
void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv )
@ -154,6 +498,7 @@ namespace Opm
}
} // end namespace detail
/// \brief A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as
///
/// This preconditioner differs from a ParallelRestrictedOverlappingSchwarz with
@ -257,10 +602,18 @@ public:
\param A The matrix to operate on.
\param n ILU fill in level (for testing). This does not work in parallel.
\param w The relaxation factor.
\param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
\param redblack Whether to use a red-black ordering.
\param reorder_sphere If true, we start the reordering at a root node.
The vertices on each layer aound it (same distance) are
ordered consecutivly. If false, we preserver the order of
the vertices with the same color.
*/
template<class BlockType, class Alloc>
ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const int n, const field_type w )
const int n, const field_type w,
MILU_VARIANT milu, bool redblack=false,
bool reorder_sphere=true)
: lower_(),
upper_(),
inv_(),
@ -269,7 +622,8 @@ public:
{
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), n );
init( reinterpret_cast<const Matrix&>(A), n, milu, redblack,
reorder_sphere );
}
/*! \brief Constructor gets all parameters to operate the prec.
@ -277,10 +631,18 @@ public:
\param comm communication object, e.g. Dune::OwnerOverlapCopyCommunication
\param n ILU fill in level (for testing). This does not work in parallel.
\param w The relaxation factor.
\param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
\param redblack Whether to use a red-black ordering.
\param reorder_sphere If true, we start the reordering at a root node.
The vertices on each layer aound it (same distance) are
ordered consecutivly. If false, we preserver the order of
the vertices with the same color.
*/
template<class BlockType, class Alloc>
ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const ParallelInfo& comm, const int n, const field_type w )
const ParallelInfo& comm, const int n, const field_type w,
MILU_VARIANT milu, bool redblack=false,
bool reorder_sphere=true)
: lower_(),
upper_(),
inv_(),
@ -289,7 +651,8 @@ public:
{
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), n );
init( reinterpret_cast<const Matrix&>(A), n, milu, redblack,
reorder_sphere );
}
/*! \brief Constructor.
@ -297,11 +660,18 @@ public:
Constructor gets all parameters to operate the prec.
\param A The matrix to operate on.
\param w The relaxation factor.
\param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
\param redblack Whether to use a red-black ordering.
\param reorder_sphere If true, we start the reordering at a root node.
The vertices on each layer aound it (same distance) are
ordered consecutivly. If false, we preserver the order of
the vertices with the same color.
*/
template<class BlockType, class Alloc>
ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const field_type w)
: ParallelOverlappingILU0( A, 0, w )
const field_type w, MILU_VARIANT milu, bool redblack=false,
bool reorder_sphere=true)
: ParallelOverlappingILU0( A, 0, w, milu, redblack, reorder_sphere )
{
}
@ -311,10 +681,18 @@ public:
\param A The matrix to operate on.
\param comm communication object, e.g. Dune::OwnerOverlapCopyCommunication
\param w The relaxation factor.
\param milu The modified ILU variant to use. 0 means traditional ILU. \see MILU_VARIANT.
\param redblack Whether to use a red-black ordering.
\param reorder_sphere If true, we start the reordering at a root node.
The vertices on each layer aound it (same distance) are
ordered consecutivly. If false, we preserver the order of
the vertices with the same color.
*/
template<class BlockType, class Alloc>
ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
const ParallelInfo& comm, const field_type w)
const ParallelInfo& comm, const field_type w,
MILU_VARIANT milu, bool redblack=false,
bool reorder_sphere=true)
: lower_(),
upper_(),
inv_(),
@ -323,7 +701,8 @@ public:
{
// BlockMatrix is a Subclass of FieldMatrix that just adds
// methods. Therefore this cast should be safe.
init( reinterpret_cast<const Matrix&>(A), 0 );
init( reinterpret_cast<const Matrix&>(A), 0, milu, redblack,
reorder_sphere );
}
/*!
@ -344,7 +723,8 @@ public:
*/
virtual void apply (Domain& v, const Range& d)
{
Range& md = const_cast<Range&>(d);
Range& md = reorderD(d);
Domain& mv = reorderV(v);
copyOwnerToAll( md );
// iterator types
@ -361,41 +741,42 @@ public:
// lower triangular solve
for( size_type i=0; i<iEnd; ++ i )
{
dblock rhs( d[ i ] );
dblock rhs( md[ i ] );
const size_type rowI = lower_.rows_[ i ];
const size_type rowINext = lower_.rows_[ i+1 ];
for( size_type col = rowI; col < rowINext; ++ col )
{
lower_.values_[ col ].mmv( v[ lower_.cols_[ col ] ], rhs );
lower_.values_[ col ].mmv( mv[ lower_.cols_[ col ] ], rhs );
}
v[ i ] = rhs; // Lii = I
mv[ i ] = rhs; // Lii = I
}
copyOwnerToAll( v );
copyOwnerToAll( mv );
for( size_type i=0; i<iEnd; ++ i )
{
vblock& vBlock = v[ lastRow - i ];
vblock& vBlock = mv[ lastRow - i ];
vblock rhs ( vBlock );
const size_type rowI = upper_.rows_[ i ];
const size_type rowINext = upper_.rows_[ i+1 ];
for( size_type col = rowI; col < rowINext; ++ col )
{
upper_.values_[ col ].mmv( v[ upper_.cols_[ col ] ], rhs );
upper_.values_[ col ].mmv( mv[ upper_.cols_[ col ] ], rhs );
}
// apply inverse and store result
inv_[ i ].mv( rhs, vBlock);
}
copyOwnerToAll( v );
copyOwnerToAll( mv );
if( relaxation_ ) {
v *= w_;
mv *= w_;
}
reorderBack(mv, v);
}
template <class V>
@ -417,7 +798,7 @@ public:
}
protected:
void init( const Matrix& A, const int iluIteration )
void init( const Matrix& A, const int iluIteration, MILU_VARIANT milu, bool redBlack, bool reorderSpheres )
{
// (For older DUNE versions the communicator might be
// invalid if redistribution in AMG happened on the coarset level.
@ -441,17 +822,103 @@ protected:
std::unique_ptr< Matrix > ILU;
if ( redBlack )
{
using Graph = Dune::Amg::MatrixGraph<const Matrix>;
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 )
{
ordering_ = reorderVerticesSpheres(colors, noColors, verticesPerColor,
graph, 0);
}
else
{
ordering_ = reorderVerticesPreserving(colors, noColors, verticesPerColor,
graph);
}
}
std::vector<std::size_t> inverseOrdering(ordering_.size());
std::size_t index = 0;
for( auto newIndex: ordering_)
{
inverseOrdering[newIndex] = index++;
}
try
{
if( iluIteration == 0 ) {
// create ILU-0 decomposition
ILU.reset( new Matrix( A ) );
bilu0_decomposition( *ILU );
if ( ordering_.empty() )
{
ILU.reset( new Matrix( A ) );
}
else
{
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()]];
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)
{
auto& newRow = newA[ordering_[iter.index()]];
for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
{
newRow[ordering_[col.index()]] = *col;
}
}
}
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:
bilu0_decomposition( *ILU );
break;
}
}
else {
// create ILU-n decomposition
ILU.reset( new Matrix( A.N(), A.M(), Matrix::row_wise) );
bilu_decomposition( A, iluIteration, *ILU );
std::unique_ptr<detail::Reorderer> reorderer, inverseReorderer;
if ( ordering_.empty() )
{
reorderer.reset(new detail::NoReorderer());
inverseReorderer.reset(new detail::NoReorderer());
}
else
{
reorderer.reset(new detail::RealReorderer(ordering_));
inverseReorderer.reset(new detail::RealReorderer(inverseOrdering));
}
milun_decomposition( A, iluIteration, milu, *ILU, *reorderer, *inverseReorderer );
}
}
catch ( Dune::MatrixBlockError error )
@ -475,11 +942,70 @@ protected:
detail::convertToCRS( *ILU, lower_, upper_, inv_ );
}
/// \brief Reorder D if needed and return a reference to it.
Range& reorderD(const Range& d)
{
if ( ordering_.empty())
{
// As d is non-const in the apply method of the
// solver casting away constness in this particular
// setting is not undefined. It is ugly though but due
// to the preconditioner interface of dune-istl.
return const_cast<Range&>(d);
}
else
{
reorderedD_.resize(d.size());
std::size_t i = 0;
for(auto index: ordering_)
{
reorderedD_[index]=d[i++];
}
return reorderedD_;
}
}
/// \brief Reorder V if needed and return a reference to it.
Domain& reorderV(Domain& v)
{
if ( ordering_.empty())
{
return v;
}
else
{
reorderedV_.resize(v.size());
std::size_t i = 0;
for(auto index: ordering_)
{
reorderedV_[index]=v[i++];
}
return reorderedV_;
}
}
void reorderBack(const Range& reorderedV, Range& v)
{
if ( !ordering_.empty() )
{
std::size_t i = 0;
for(auto index: ordering_)
{
v[i++] = reorderedV[index];
}
}
}
protected:
//! \brief The ILU0 decomposition of the matrix.
CRS lower_;
CRS upper_;
std::vector< block_type > inv_;
//! \brief the reordering of the unknowns
std::vector< std::size_t > ordering_;
//! \brief The reordered right hand side
Range reorderedD_;
//! \brief The reordered left hand side.
Domain reorderedV_;
const ParallelInfo* comm_;
//! \brief The relaxation factor to use.

101
tests/test_graphcoloring.cpp Normal file
View File

@ -0,0 +1,101 @@
#include <config.h>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/paamg/graph.hh>
#include <opm/autodiff/GraphColoring.hpp>
#define BOOST_TEST_MODULE GraphColoringTest
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>
///! \brief check that all indices are represented in the new ordering.
void checkAllIndices(const std::vector<std::size_t>& ordering)
{
std::vector<int> counters(ordering.size(), 0);
for(auto index: ordering)
{
++counters[index];
}
for(auto count: counters)
{
BOOST_CHECK(count==1);
}
}
BOOST_AUTO_TEST_CASE(TestWelschPowell)
{
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
using Graph = Dune::Amg::MatrixGraph<Matrix>;
int N = 10;
Matrix matrix(N*N, N*N, 5, 0.4, Matrix::implicit);
for( int j = 0; j < N; j++)
{
for(int i = 0; i < N; i++)
{
auto index = j*10+i;
matrix.entry(index,index) = 1;
if ( i > 0 )
{
matrix.entry(index,index-1) = 1;
}
if ( i < N - 1)
{
matrix.entry(index,index+1) = 1;
}
if ( j > 0 )
{
matrix.entry(index,index-N) = 1;
}
if ( j < N - 1)
{
matrix.entry(index,index+N) = 1;
}
}
}
matrix.compress();
Graph graph(matrix);
auto colorsTuple = Opm::colorVerticesWelshPowell(graph);
const auto& colors = std::get<0>(colorsTuple);
const auto& verticesPerColor = std::get<2>(colorsTuple);
auto noColors = std::get<1>(colorsTuple);
auto firstCornerColor = colors[0];
BOOST_CHECK(noColors == 2);
// Check for checkerboard coloring
for( int j = 0, index = 0; j < N; j++)
{
auto expectedColor = firstCornerColor;
for(int i = 0; i < N; i++)
{
BOOST_CHECK(colors[index]==expectedColor);
index++;
expectedColor = (expectedColor + 1) % 2;
}
firstCornerColor=(firstCornerColor + 1) % 2;
}
auto newOrder = Opm::reorderVerticesPreserving(colors, noColors, verticesPerColor,
graph);
std::vector<std::size_t> colorIndex(noColors, 0);
std::partial_sum(verticesPerColor.begin(),
verticesPerColor.begin()+verticesPerColor.size()-1,
colorIndex.begin()+1);
for (auto vertex : graph)
{
BOOST_CHECK(colorIndex[colors[vertex]]++ == newOrder[vertex]);
}
checkAllIndices(newOrder);
newOrder = Opm::reorderVerticesSpheres(colors, noColors, verticesPerColor,
graph, 0);
checkAllIndices(newOrder);
}

193
tests/test_milu.cpp Normal file
View File

@ -0,0 +1,193 @@
#include<config.h>
#define BOOST_TEST_MODULE MILU0Test
#include<vector>
#include<memory>
#include<dune/istl/bcrsmatrix.hh>
#include<dune/istl/bvector.hh>
#include<dune/common/fmatrix.hh>
#include<dune/common/fvector.hh>
#include<opm/autodiff/ParallelOverlappingILU0.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
template<class M>
void test_milu0(M& A)
{
auto ILU = A;
std::shared_ptr<std::vector<typename M::block_type>> diagonal = nullptr;
diagonal.reset(new std::vector<typename M::block_type>());
Opm::detail::milu0_decomposition(ILU, diagonal.get());
#ifdef DEBUG
if ( A.N() < 11)
{
Dune::printmatrix(std::cout, ILU, "ILU", "row");
std::cout << "Diagonal: ";
for (const auto& d : *diagonal)
std::cout << d << " ";
std::cout<<std::endl;
}
#endif
Dune::BlockVector<Dune::FieldVector<typename M::field_type, M::block_type::rows>> e(A.N()), x1(A.N()), x2(A.N()), t(A.N());
e = 1;
A.mv(e, x1);
// Compute L*U*e;
// t=Ue
for ( auto irow = ILU.begin(), iend = ILU.end(); irow != iend; ++irow)
{
auto i = irow.index();
(*diagonal)[i].mv(e[0], t[i]);
auto col = irow->find(i);
auto colend = irow->end();
if ( col == colend )
{
OPM_THROW(std::logic_error, "Missing diagonal entry for row " << irow.index());
}
for (++col ;col != colend; ++col)
{
col->umv(e[0], t[i]);
}
}
for ( auto irow = ILU.begin(), iend = ILU.end(); irow != iend; ++irow)
{
auto i = irow.index();
x2[i] = 0;
for (auto col = irow->begin(); col.index() < irow.index(); ++col)
{
col->umv(t[col.index()], x2[i]);
}
x2[i] += t[i];
}
auto diff = x2;
diff-=x1;
for ( std::size_t i = 0, end = A.N(); i < end; ++i)
{
auto point_difference = diff[i].two_norm();
#ifdef DEBUG
std::cout<<"index "<<i<<" size "<<diff.size()<<" difference"<<point_difference<<std::endl;
#endif
BOOST_CHECK(point_difference < 1e-12);
}
// Test that (LU)^-1Ae=e
A.mv(e, x1);
bilu_backsolve(ILU, x2, x1);
diff = x2;
diff -= e;
for ( std::size_t i = 0, end = A.N(); i < end; ++i)
{
#ifdef DEBUG
auto point_difference = diff[i].two_norm();
std::cout<<"index "<<i<<" size "<<diff.size()<<" point_difference "<<point_difference<<std::endl;
#endif
BOOST_CHECK_CLOSE(x2[i].two_norm(), e[i].two_norm(), 1e-12);
}
}
template<class B, class Alloc>
void setupSparsityPattern(Dune::BCRSMatrix<B,Alloc>& A, int N)
{
typedef typename Dune::BCRSMatrix<B,Alloc> Matrix;
A.setSize(N*N, N*N, N*N*5);
A.setBuildMode(Matrix::row_wise);
for (typename Dune::BCRSMatrix<B,Alloc>::CreateIterator i = A.createbegin(); i != A.createend(); ++i) {
int x = i.index()%N; // x coordinate in the 2d field
int y = i.index()/N; // y coordinate in the 2d field
if(y>0)
// insert lower neighbour
i.insert(i.index()-N);
if(x>0)
// insert left neighbour
i.insert(i.index()-1);
// insert diagonal value
i.insert(i.index());
if(x<N-1)
//insert right neighbour
i.insert(i.index()+1);
if(y<N-1)
// insert upper neighbour
i.insert(i.index()+N);
}
}
template<class B, class Alloc>
void setupLaplacian(Dune::BCRSMatrix<B,Alloc>& A, int N)
{
typedef typename Dune::BCRSMatrix<B,Alloc>::field_type FieldType;
setupSparsityPattern(A,N);
B diagonal(static_cast<FieldType>(0)), bone(static_cast<FieldType>(0)),
beps(static_cast<FieldType>(0));
for(typename B::RowIterator b = diagonal.begin(); b != diagonal.end(); ++b)
b->operator[](b.index())=4;
for(typename B::RowIterator b=bone.begin(); b != bone.end(); ++b)
b->operator[](b.index())=-1.0;
for (typename Dune::BCRSMatrix<B,Alloc>::RowIterator i = A.begin(); i != A.end(); ++i) {
int x = i.index()%N; // x coordinate in the 2d field
int y = i.index()/N; // y coordinate in the 2d field
i->operator[](i.index())=diagonal;
if(y>0)
i->operator[](i.index()-N)=bone;
if(y<N-1)
i->operator[](i.index()+N)=bone;
if(x>0)
i->operator[](i.index()-1)=bone;
if(x < N-1)
i->operator[](i.index()+1)=bone;
}
}
template<int bsize>
void test()
{
std::size_t N = 32;
Dune::BCRSMatrix<Dune::FieldMatrix<double, bsize, bsize> > A;
setupLaplacian(A, N);
test_milu0(A);
#ifdef DEBUG
std::cout<< "Tested block size "<< bsize<<std::endl;
#endif
}
BOOST_AUTO_TEST_CASE(MILULaplace1)
{
test<1>();
}
BOOST_AUTO_TEST_CASE(MILULaplace2)
{
test<2>();
}
BOOST_AUTO_TEST_CASE(MILULaplace3)
{
test<3>();
}
BOOST_AUTO_TEST_CASE(MILULaplace4)
{
test<4>();
}