mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
First version of a AMG for the Blackoil equations.
The approach is inspired by Geiger's system-amg but we use dune-istl aggregation AMG for it. On the fine level all unknowns attached to a cell form a matrix block and are treated fully coupled. To form the first coarse level system we use only the pressure component to guide the aggregation and neglect all other unknowns on the fine level. All other level are formed in the usual way by scalar aggregation. Currently,it has to be requested for flow_ebos manually by passing "linear_solver_use_amg=true amg_blackoil_system=true" to it.
This commit is contained in:
Markus Blatt
parent
a2419e285f
commit
f5d81513da
@ -137,6 +137,7 @@ list (APPEND MAIN_SOURCE_FILES
|
||||
list (APPEND TEST_SOURCE_FILES
|
||||
tests/test_autodiffhelpers.cpp
|
||||
tests/test_autodiffmatrix.cpp
|
||||
tests/test_blackoil_amg.cpp
|
||||
tests/test_block.cpp
|
||||
tests/test_boprops_ad.cpp
|
||||
tests/test_rateconverter.cpp
|
||||
@ -256,6 +257,7 @@ list (APPEND PUBLIC_HEADER_FILES
|
||||
opm/autodiff/AutoDiffHelpers.hpp
|
||||
opm/autodiff/AutoDiffMatrix.hpp
|
||||
opm/autodiff/AutoDiff.hpp
|
||||
opm/autodiff/BlackoilAmg.hpp
|
||||
opm/autodiff/BlackoilDetails.hpp
|
||||
opm/autodiff/BlackoilLegacyDetails.hpp
|
||||
opm/autodiff/BlackoilModel.hpp
|
||||
|
||||
727
opm/autodiff/BlackoilAmg.hpp
Normal file
727
opm/autodiff/BlackoilAmg.hpp
Normal file
@ -0,0 +1,727 @@
|
||||
/*
|
||||
Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
|
||||
|
||||
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_AMG_HEADER_INCLUDED
|
||||
#define OPM_AMG_HEADER_INCLUDED
|
||||
|
||||
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
|
||||
#include <opm/autodiff/CPRPreconditioner.hpp>
|
||||
#include <dune/istl/paamg/twolevelmethod.hh>
|
||||
#include <dune/istl/paamg/aggregates.hh>
|
||||
#include <dune/istl/bvector.hh>
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/istl/preconditioners.hh>
|
||||
#include <dune/istl/schwarz.hh>
|
||||
#include <dune/istl/scalarproducts.hh>
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
namespace Dune
|
||||
{
|
||||
namespace Amg
|
||||
{
|
||||
template<class M, class Norm>
|
||||
class UnSymmetricCriterion;
|
||||
}
|
||||
}
|
||||
|
||||
namespace Dune
|
||||
{
|
||||
|
||||
template <class Scalar, int n, int m>
|
||||
class MatrixBlock;
|
||||
|
||||
}
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
namespace Detail
|
||||
{
|
||||
|
||||
template<typename NonScalarType>
|
||||
struct ScalarType
|
||||
{
|
||||
};
|
||||
|
||||
template<typename FieldType, int SIZE>
|
||||
struct ScalarType<Dune::FieldVector<FieldType, SIZE> >
|
||||
{
|
||||
typedef Dune::FieldVector<FieldType, 1> value;
|
||||
};
|
||||
|
||||
template<typename FieldType, int ROWS, int COLS>
|
||||
struct ScalarType<Dune::FieldMatrix<FieldType, ROWS, COLS> >
|
||||
{
|
||||
typedef Dune::FieldMatrix<FieldType, 1, 1> value;
|
||||
};
|
||||
|
||||
template<typename FieldType, int ROWS, int COLS>
|
||||
struct ScalarType<Dune::MatrixBlock<FieldType, ROWS, COLS> >
|
||||
{
|
||||
typedef Dune::MatrixBlock<FieldType, 1, 1> value;
|
||||
};
|
||||
|
||||
template<typename BlockType, typename Allocator>
|
||||
struct ScalarType<Dune::BCRSMatrix<BlockType, Allocator> >
|
||||
{
|
||||
using ScalarBlock = typename ScalarType<BlockType>::value;
|
||||
using ScalarAllocator = typename Allocator::template rebind<ScalarBlock>::other;
|
||||
typedef Dune::BCRSMatrix<ScalarBlock,ScalarAllocator> value;
|
||||
};
|
||||
|
||||
template<typename BlockType, typename Allocator>
|
||||
struct ScalarType<Dune::BlockVector<BlockType, Allocator> >
|
||||
{
|
||||
using ScalarBlock = typename ScalarType<BlockType>::value;
|
||||
using ScalarAllocator = typename Allocator::template rebind<ScalarBlock>::other;
|
||||
typedef Dune::BlockVector<ScalarBlock,ScalarAllocator> value;
|
||||
};
|
||||
|
||||
template<typename X>
|
||||
struct ScalarType<Dune::SeqScalarProduct<X> >
|
||||
{
|
||||
typedef Dune::SeqScalarProduct<typename ScalarType<X>::value> value;
|
||||
};
|
||||
|
||||
#define ComposeScalarTypeForSeqPrecond(PREC) \
|
||||
template<typename M, typename X, typename Y, int l> \
|
||||
struct ScalarType<PREC<M,X,Y,l> > \
|
||||
{ \
|
||||
typedef PREC<typename ScalarType<M>::value, \
|
||||
typename ScalarType<X>::value, \
|
||||
typename ScalarType<Y>::value, \
|
||||
l> value; \
|
||||
};
|
||||
|
||||
ComposeScalarTypeForSeqPrecond(Dune::SeqJac);
|
||||
ComposeScalarTypeForSeqPrecond(Dune::SeqSOR);
|
||||
ComposeScalarTypeForSeqPrecond(Dune::SeqSSOR);
|
||||
ComposeScalarTypeForSeqPrecond(Dune::SeqGS);
|
||||
ComposeScalarTypeForSeqPrecond(Dune::SeqILU0);
|
||||
ComposeScalarTypeForSeqPrecond(Dune::SeqILUn);
|
||||
|
||||
#undef ComposeScalarTypeForSeqPrecond
|
||||
|
||||
template<typename X, typename Y>
|
||||
struct ScalarType<Dune::Richardson<X,Y> >
|
||||
{
|
||||
typedef Dune::Richardson<typename ScalarType<X>::value,
|
||||
typename ScalarType<Y>::value> value;
|
||||
};
|
||||
|
||||
template<class M, class X, class Y, class C>
|
||||
struct ScalarType<Dune::OverlappingSchwarzOperator<M,X,Y,C> >
|
||||
{
|
||||
typedef Dune::OverlappingSchwarzOperator<typename ScalarType<M>::value,
|
||||
typename ScalarType<X>::value,
|
||||
typename ScalarType<Y>::value,
|
||||
C> value;
|
||||
};
|
||||
|
||||
template<class M, class X, class Y>
|
||||
struct ScalarType<Dune::MatrixAdapter<M,X,Y> >
|
||||
{
|
||||
typedef Dune::MatrixAdapter<typename ScalarType<M>::value,
|
||||
typename ScalarType<X>::value,
|
||||
typename ScalarType<Y>::value> value;
|
||||
};
|
||||
|
||||
template<class X, class C>
|
||||
struct ScalarType<Dune::OverlappingSchwarzScalarProduct<X,C> >
|
||||
{
|
||||
typedef Dune::OverlappingSchwarzScalarProduct<typename ScalarType<X>::value,
|
||||
C> value;
|
||||
};
|
||||
|
||||
template<class X, class C>
|
||||
struct ScalarType<Dune::NonoverlappingSchwarzScalarProduct<X,C> >
|
||||
{
|
||||
typedef Dune::NonoverlappingSchwarzScalarProduct<typename ScalarType<X>::value,
|
||||
C> value;
|
||||
};
|
||||
|
||||
template<class X, class Y, class C, class T>
|
||||
struct ScalarType<Dune::BlockPreconditioner<X,Y,C,T> >
|
||||
{
|
||||
typedef Dune::BlockPreconditioner<typename ScalarType<X>::value,
|
||||
typename ScalarType<Y>::value,
|
||||
C,
|
||||
typename ScalarType<T>::value> value;
|
||||
};
|
||||
|
||||
template<class M, class X, class Y, class C>
|
||||
struct ScalarType<ParallelOverlappingILU0<M,X,Y,C> >
|
||||
{
|
||||
typedef ParallelOverlappingILU0<typename ScalarType<M>::value,
|
||||
typename ScalarType<X>::value,
|
||||
typename ScalarType<Y>::value,
|
||||
C> value;
|
||||
};
|
||||
|
||||
template<class B, class N>
|
||||
struct ScalarType<Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<B>,N> > >
|
||||
{
|
||||
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<typename ScalarType<B>::value>, Dune::Amg::FirstDiagonal> >;
|
||||
};
|
||||
|
||||
template<class B, class N>
|
||||
struct ScalarType<Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<B>,N> > >
|
||||
{
|
||||
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<typename ScalarType<B>::value>, Dune::Amg::FirstDiagonal> >;
|
||||
};
|
||||
|
||||
template<class C, std::size_t COMPONENT_INDEX>
|
||||
struct OneComponentCriterionType
|
||||
{};
|
||||
|
||||
template<class B, class N, std::size_t COMPONENT_INDEX>
|
||||
struct OneComponentCriterionType<Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<B>,N> >,COMPONENT_INDEX>
|
||||
{
|
||||
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<Dune::BCRSMatrix<B>, Dune::Amg::Diagonal<COMPONENT_INDEX> > >;
|
||||
};
|
||||
|
||||
template<class B, class N, std::size_t COMPONENT_INDEX>
|
||||
struct OneComponentCriterionType<Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<B>,N> >,COMPONENT_INDEX>
|
||||
{
|
||||
using value = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Dune::BCRSMatrix<B>, Dune::Amg::Diagonal<COMPONENT_INDEX> > >;
|
||||
};
|
||||
|
||||
template<class Operator, class Criterion, class Communication, std::size_t COMPONENT_INDEX>
|
||||
class OneComponentAggregationLevelTransferPolicy;
|
||||
|
||||
|
||||
/**
|
||||
* @brief A policy class for solving the coarse level system using one step of AMG.
|
||||
* @tparam O The type of the linear operator used.
|
||||
* @tparam S The type of the smoother used in AMG.
|
||||
* @tparam C The type of the crition used for the aggregation within AMG.
|
||||
* @tparam C1 The type of the information about the communication. Either
|
||||
* Dune::OwnerOverlapCopyCommunication or Dune::SequentialInformation.
|
||||
*/
|
||||
template<class O, class S, class C, class P>
|
||||
class OneStepAMGCoarseSolverPolicy
|
||||
{
|
||||
public:
|
||||
typedef P LevelTransferPolicy;
|
||||
/** @brief The type of the linear operator used. */
|
||||
typedef O Operator;
|
||||
/** @brief The type of the range and domain of the operator. */
|
||||
typedef typename O::range_type X;
|
||||
/** @brief The type of the crition used for the aggregation within AMG.*/
|
||||
typedef C Criterion;
|
||||
/** @brief The type of the communication used for AMG.*/
|
||||
typedef typename P::ParallelInformation Communication;
|
||||
/** @brief The type of the smoother used in AMG. */
|
||||
typedef S Smoother;
|
||||
/** @brief The type of the arguments used for constructing the smoother. */
|
||||
typedef typename Dune::Amg::SmootherTraits<S>::Arguments SmootherArgs;
|
||||
/** @brief The type of the AMG construct on the coarse level.*/
|
||||
typedef Dune::Amg::AMG<Operator,X,Smoother,Communication> AMGType;
|
||||
/**
|
||||
* @brief Constructs the coarse solver policy.
|
||||
* @param args The arguments used for constructing the smoother.
|
||||
* @param c The crition used for the aggregation within AMG.
|
||||
*/
|
||||
OneStepAMGCoarseSolverPolicy(const SmootherArgs& args, const Criterion& c)
|
||||
: smootherArgs_(args), criterion_(c)
|
||||
{}
|
||||
/** @brief Copy constructor. */
|
||||
OneStepAMGCoarseSolverPolicy(const OneStepAMGCoarseSolverPolicy& other)
|
||||
: coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_),
|
||||
criterion_(other.criterion_)
|
||||
{}
|
||||
private:
|
||||
/**
|
||||
* @brief A wrapper that makes an inverse operator out of AMG.
|
||||
*
|
||||
* The operator will use one step of AMG to approximately solve
|
||||
* the coarse level system.
|
||||
*/
|
||||
struct AMGInverseOperator : public Dune::InverseOperator<X,X>
|
||||
{
|
||||
AMGInverseOperator(const typename AMGType::Operator& op,
|
||||
const Criterion& crit,
|
||||
const typename AMGType::SmootherArgs& args,
|
||||
const Communication& comm)
|
||||
: amg_(op, crit,args, comm), op_(op), comm_(comm), first_(true)
|
||||
{}
|
||||
|
||||
void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res)
|
||||
{
|
||||
DUNE_UNUSED_PARAMETER(reduction);
|
||||
DUNE_UNUSED_PARAMETER(res);
|
||||
/*
|
||||
if(first_)
|
||||
{
|
||||
amg_.pre(x,b);
|
||||
first_=false;
|
||||
x_=x;
|
||||
}
|
||||
amg_.apply(x,b);
|
||||
*/
|
||||
using Chooser = Dune::ScalarProductChooser<X,Communication,AMGType::category>;
|
||||
auto sp = Chooser::construct(comm_);
|
||||
Dune::BiCGSTABSolver<X> solver(const_cast<typename AMGType::Operator&>(op_), *sp, amg_, 1e-2, 25, 0);
|
||||
solver.apply(x,b,res);
|
||||
delete sp;
|
||||
}
|
||||
|
||||
void apply(X& x, X& b, Dune::InverseOperatorResult& res)
|
||||
{
|
||||
return apply(x,b,1e-8,res);
|
||||
}
|
||||
|
||||
~AMGInverseOperator()
|
||||
{
|
||||
/*
|
||||
if(!first_)
|
||||
amg_.post(x_);
|
||||
*/
|
||||
}
|
||||
AMGInverseOperator(const AMGInverseOperator& other)
|
||||
: x_(other.x_), amg_(other.amg_), first_(other.first_)
|
||||
{
|
||||
}
|
||||
private:
|
||||
X x_;
|
||||
AMGType amg_;
|
||||
const typename AMGType::Operator& op_;
|
||||
const Communication& comm_;
|
||||
bool first_;
|
||||
};
|
||||
|
||||
public:
|
||||
/** @brief The type of solver constructed for the coarse level. */
|
||||
typedef AMGInverseOperator CoarseLevelSolver;
|
||||
|
||||
/**
|
||||
* @brief Constructs a coarse level solver.
|
||||
*
|
||||
* @param transferPolicy The policy describing the transfer between levels.
|
||||
* @return A pointer to the constructed coarse level solver.
|
||||
*/
|
||||
template<class LTP>
|
||||
CoarseLevelSolver* createCoarseLevelSolver(LTP& transferPolicy)
|
||||
{
|
||||
coarseOperator_=transferPolicy.getCoarseLevelOperator();
|
||||
const LevelTransferPolicy& transfer =
|
||||
reinterpret_cast<const LevelTransferPolicy&>(transferPolicy);
|
||||
AMGInverseOperator* inv = new AMGInverseOperator(*coarseOperator_,
|
||||
criterion_,
|
||||
smootherArgs_,
|
||||
transfer.getCoarseLevelCommunication());
|
||||
|
||||
return inv; //std::shared_ptr<InverseOperator<X,X> >(inv);
|
||||
|
||||
}
|
||||
|
||||
private:
|
||||
/** @brief The coarse level operator. */
|
||||
std::shared_ptr<Operator> coarseOperator_;
|
||||
/** @brief The arguments used to construct the smoother. */
|
||||
SmootherArgs smootherArgs_;
|
||||
/** @brief The coarsening criterion. */
|
||||
Criterion criterion_;
|
||||
};
|
||||
|
||||
template<class Smoother, class Operator, class Communication>
|
||||
Smoother* constructSmoother(const Operator& op,
|
||||
const typename Dune::Amg::SmootherTraits<Smoother>::Arguments& smargs,
|
||||
const Communication& comm)
|
||||
{
|
||||
typename Dune::Amg::ConstructionTraits<Smoother>::Arguments args;
|
||||
args.setMatrix(op.getmat());
|
||||
args.setComm(comm);
|
||||
args.setArgs(smargs);
|
||||
return Dune::Amg::ConstructionTraits<Smoother>::construct(args);
|
||||
}
|
||||
|
||||
template<class G, class C, class S>
|
||||
const Dune::Amg::OverlapVertex<typename G::VertexDescriptor>*
|
||||
buildOverlapVertices(const G& graph, const C& pinfo,
|
||||
Dune::Amg::AggregatesMap<typename G::VertexDescriptor>& aggregates,
|
||||
const S& overlap,
|
||||
std::size_t& overlapCount)
|
||||
{
|
||||
// count the overlap vertices.
|
||||
auto end = graph.end();
|
||||
overlapCount = 0;
|
||||
|
||||
const auto& lookup=pinfo.globalLookup();
|
||||
|
||||
for(auto vertex=graph.begin(); vertex != end; ++vertex) {
|
||||
const auto* pair = lookup.pair(*vertex);
|
||||
|
||||
if(pair!=0 && overlap.contains(pair->local().attribute()))
|
||||
++overlapCount;
|
||||
}
|
||||
// Allocate space
|
||||
using Vertex = typename G::VertexDescriptor;
|
||||
using OverlapVertex = Dune::Amg::OverlapVertex<Vertex>;
|
||||
|
||||
auto* overlapVertices = new OverlapVertex[overlapCount=0 ? 1 : overlapCount];
|
||||
if(overlapCount==0)
|
||||
return overlapVertices;
|
||||
|
||||
// Initialize them
|
||||
overlapCount=0;
|
||||
for(auto vertex=graph.begin(); vertex != end; ++vertex) {
|
||||
const auto* pair = lookup.pair(*vertex);
|
||||
|
||||
if(pair!=0 && overlap.contains(pair->local().attribute())) {
|
||||
overlapVertices[overlapCount].aggregate = &aggregates[pair->local()];
|
||||
overlapVertices[overlapCount].vertex = pair->local();
|
||||
++overlapCount;
|
||||
}
|
||||
}
|
||||
std::sort(overlapVertices, overlapVertices+overlapCount,
|
||||
[](const OverlapVertex& v1, const OverlapVertex& v2)
|
||||
{
|
||||
return *v1.aggregate < *v2.aggregate;
|
||||
});
|
||||
// due to the sorting the isolated aggregates (to be skipped) are at the end.
|
||||
|
||||
return overlapVertices;
|
||||
}
|
||||
|
||||
template<class M, class G, class V, class C, class S>
|
||||
void buildCoarseSparseMatrix(M& coarseMatrix, G& fineGraph,
|
||||
const V& visitedMap,
|
||||
const C& pinfo,
|
||||
Dune::Amg::AggregatesMap<typename G::VertexDescriptor>& aggregates,
|
||||
const S& overlap)
|
||||
{
|
||||
using OverlapVertex = Dune::Amg ::OverlapVertex<typename G::VertexDescriptor>;
|
||||
std::size_t count;
|
||||
|
||||
const OverlapVertex* overlapVertices = buildOverlapVertices(fineGraph,
|
||||
pinfo,
|
||||
aggregates,
|
||||
overlap,
|
||||
count);
|
||||
|
||||
// Reset the visited flags of all vertices.
|
||||
// As the isolated nodes will be skipped we simply mark them as visited
|
||||
#ifndef NDEBUG
|
||||
const auto UNAGGREGATED = Dune::Amg::AggregatesMap<typename G::VertexDescriptor>::UNAGGREGATED;
|
||||
#endif
|
||||
const auto ISOLATED = Dune::Amg::AggregatesMap<typename G::VertexDescriptor>::ISOLATED;
|
||||
auto vend = fineGraph.end();
|
||||
|
||||
for(auto vertex = fineGraph.begin(); vertex != vend; ++vertex) {
|
||||
assert(aggregates[*vertex] != UNAGGREGATED);
|
||||
put(visitedMap, *vertex, aggregates[*vertex]==ISOLATED);
|
||||
}
|
||||
|
||||
Dune::Amg::SparsityBuilder<M> sparsityBuilder(coarseMatrix);
|
||||
|
||||
Dune::Amg::ConnectivityConstructor<G,C>::examine(fineGraph, visitedMap, pinfo,
|
||||
aggregates, overlap,
|
||||
overlapVertices,
|
||||
overlapVertices+count,
|
||||
sparsityBuilder);
|
||||
delete[] overlapVertices;
|
||||
}
|
||||
|
||||
template<class M, class G, class V, class S>
|
||||
void buildCoarseSparseMatrix(M& coarseMatrix, G& fineGraph, const V& visitedMap,
|
||||
const Dune::Amg::SequentialInformation& pinfo,
|
||||
Dune::Amg::AggregatesMap<typename G::VertexDescriptor>& aggregates,
|
||||
const S&)
|
||||
{
|
||||
// Reset the visited flags of all vertices.
|
||||
// As the isolated nodes will be skipped we simply mark them as visited
|
||||
#ifndef NDEBUG
|
||||
const auto UNAGGREGATED = Dune::Amg::AggregatesMap<typename G::VertexDescriptor>::UNAGGREGATED;
|
||||
#endif
|
||||
const auto ISOLATED = Dune::Amg::AggregatesMap<typename G::VertexDescriptor>::ISOLATED;
|
||||
using Vertex = typename G::VertexIterator;
|
||||
Vertex vend = fineGraph.end();
|
||||
for(Vertex vertex = fineGraph.begin(); vertex != vend; ++vertex) {
|
||||
assert(aggregates[*vertex] != UNAGGREGATED);
|
||||
put(visitedMap, *vertex, aggregates[*vertex]==ISOLATED);
|
||||
}
|
||||
|
||||
Dune::Amg::SparsityBuilder<M> sparsityBuilder(coarseMatrix);
|
||||
|
||||
Dune::Amg::ConnectivityConstructor<G,Dune::Amg::SequentialInformation>
|
||||
::examine(fineGraph, visitedMap, pinfo, aggregates, sparsityBuilder);
|
||||
}
|
||||
|
||||
} // end namespace Detail
|
||||
|
||||
/**
|
||||
* @brief A LevelTransferPolicy that uses aggregation to construct the coarse level system.
|
||||
* @tparam Operator The type of the fine level operator.
|
||||
* @tparam Criterion The criterion that describes the aggregation procedure.
|
||||
* @tparam Communication The class that describes the communication pattern.
|
||||
*/
|
||||
template<class Operator, class Criterion, class Communication, std::size_t COMPONENT_INDEX>
|
||||
class OneComponentAggregationLevelTransferPolicy
|
||||
: public Dune::Amg::LevelTransferPolicy<Operator, typename Detail::ScalarType<Operator>::value>
|
||||
{
|
||||
typedef Dune::Amg::AggregatesMap<typename Operator::matrix_type::size_type> AggregatesMap;
|
||||
public:
|
||||
using CoarseOperator = typename Detail::ScalarType<Operator>::value;
|
||||
typedef Dune::Amg::LevelTransferPolicy<Operator,CoarseOperator> FatherType;
|
||||
typedef Communication ParallelInformation;
|
||||
|
||||
public:
|
||||
OneComponentAggregationLevelTransferPolicy(const Criterion& crit, const Communication& comm)
|
||||
: criterion_(crit), communication_(&const_cast<Communication&>(comm))
|
||||
{}
|
||||
|
||||
void createCoarseLevelSystem(const Operator& fineOperator)
|
||||
{
|
||||
prolongDamp_ = 1;
|
||||
typedef Dune::Amg::PropertiesGraphCreator<Operator> GraphCreator;
|
||||
typedef typename GraphCreator::PropertiesGraph PropertiesGraph;
|
||||
typedef typename GraphCreator::GraphTuple GraphTuple;
|
||||
|
||||
typedef typename PropertiesGraph::VertexDescriptor Vertex;
|
||||
|
||||
std::vector<bool> excluded(fineOperator.getmat().N(), false);
|
||||
|
||||
using OverlapFlags = Dune::NegateSet<typename ParallelInformation::OwnerSet>;
|
||||
GraphTuple graphs = GraphCreator::create(fineOperator, excluded,
|
||||
*communication_, OverlapFlags());
|
||||
|
||||
aggregatesMap_.reset(new AggregatesMap(std::get<1>(graphs)->maxVertex()+1));
|
||||
|
||||
int noAggregates, isoAggregates, oneAggregates, skippedAggregates;
|
||||
using std::get;
|
||||
std::tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) =
|
||||
aggregatesMap_->buildAggregates(fineOperator.getmat(), *(get<1>(graphs)),
|
||||
criterion_, true);
|
||||
std::cout<<"no aggregates="<<noAggregates<<" iso="<<isoAggregates<<" one="<<oneAggregates<<" skipped="<<skippedAggregates<<std::endl;
|
||||
|
||||
using CommunicationArgs = typename Dune::Amg::ConstructionTraits<Communication>::Arguments;
|
||||
CommunicationArgs commArgs(communication_->communicator(), communication_->getSolverCategory());
|
||||
coarseLevelCommunication_.reset(Dune::Amg::ConstructionTraits<Communication>::construct(commArgs));
|
||||
using Iterator = typename std::vector<bool>::iterator;
|
||||
auto visitedMap = get(Dune::Amg::VertexVisitedTag(), *(get<1>(graphs)));
|
||||
communication_->buildGlobalLookup(fineOperator.getmat().N());
|
||||
std::size_t aggregates =
|
||||
Dune::Amg::IndicesCoarsener<ParallelInformation,OverlapFlags>
|
||||
::coarsen(*communication_, *get<1>(graphs), visitedMap,
|
||||
*aggregatesMap_, *coarseLevelCommunication_,
|
||||
noAggregates);
|
||||
GraphCreator::free(graphs);
|
||||
coarseLevelCommunication_->buildGlobalLookup(aggregates);
|
||||
Dune::Amg::AggregatesPublisher<Vertex,OverlapFlags,ParallelInformation>
|
||||
::publish(*aggregatesMap_,
|
||||
*communication_,
|
||||
coarseLevelCommunication_->globalLookup());
|
||||
std::vector<bool>& visited=excluded;
|
||||
|
||||
std::fill(visited.begin(), visited.end(), false);
|
||||
|
||||
Dune::IteratorPropertyMap<Iterator, Dune::IdentityMap>
|
||||
visitedMap2(visited.begin(), Dune::IdentityMap());
|
||||
using CoarseMatrix = typename CoarseOperator::matrix_type;
|
||||
coarseLevelMatrix_.reset(new CoarseMatrix(aggregates, aggregates,
|
||||
CoarseMatrix::row_wise));
|
||||
Detail::buildCoarseSparseMatrix(*coarseLevelMatrix_, *get<0>(graphs), visitedMap2,
|
||||
*communication_,
|
||||
*aggregatesMap_,
|
||||
OverlapFlags());
|
||||
delete get<0>(graphs);
|
||||
communication_->freeGlobalLookup();
|
||||
if( static_cast<int>(this->coarseLevelMatrix_->N())
|
||||
< criterion_.coarsenTarget())
|
||||
{
|
||||
coarseLevelCommunication_->freeGlobalLookup();
|
||||
}
|
||||
calculateCoarseEntries(fineOperator.getmat());
|
||||
this->lhs_.resize(this->coarseLevelMatrix_->M());
|
||||
this->rhs_.resize(this->coarseLevelMatrix_->N());
|
||||
using OperatorArgs = typename Dune::Amg::ConstructionTraits<CoarseOperator>::Arguments;
|
||||
OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_);
|
||||
this->operator_.reset(Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs));
|
||||
}
|
||||
|
||||
template<class M>
|
||||
void calculateCoarseEntries(const M& fineMatrix)
|
||||
{
|
||||
*coarseLevelMatrix_ = 0;
|
||||
for(auto row = fineMatrix.begin(), rowEnd = fineMatrix.end();
|
||||
row != rowEnd; ++row)
|
||||
{
|
||||
const auto& i = (*aggregatesMap_)[row.index()];
|
||||
if(i != AggregatesMap::ISOLATED)
|
||||
{
|
||||
for(auto entry = row->begin(), entryEnd = row->end();
|
||||
entry != entryEnd; ++entry)
|
||||
{
|
||||
const auto& j = (*aggregatesMap_)[entry.index()];
|
||||
if ( j != AggregatesMap::ISOLATED )
|
||||
{
|
||||
(*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][COMPONENT_INDEX];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void moveToCoarseLevel(const typename FatherType::FineRangeType& fine)
|
||||
{
|
||||
// Set coarse vector to zero
|
||||
this->rhs_=0;
|
||||
|
||||
auto end = fine.end(), begin=fine.begin();
|
||||
for(auto block=begin; block != end; ++block)
|
||||
{
|
||||
const auto& vertex = (*aggregatesMap_)[block-begin];
|
||||
if(vertex != AggregatesMap::ISOLATED)
|
||||
{
|
||||
this->rhs_[vertex] += (*block)[COMPONENT_INDEX];
|
||||
}
|
||||
}
|
||||
this->lhs_=0;
|
||||
}
|
||||
|
||||
void moveToFineLevel(typename FatherType::FineDomainType& fine)
|
||||
{
|
||||
//for(auto& entry: this->lhs_)
|
||||
// entry*=prolongDamp_;
|
||||
this->lhs_ *= prolongDamp_;
|
||||
auto end=fine.end(), begin=fine.begin();
|
||||
|
||||
for(auto block=begin; block != end; ++block)
|
||||
{
|
||||
const auto& vertex = (*aggregatesMap_)[block-begin];
|
||||
if(vertex != AggregatesMap::ISOLATED)
|
||||
(*block)[COMPONENT_INDEX] += this->lhs_[vertex];
|
||||
}
|
||||
communication_->copyOwnerToAll(fine,fine);
|
||||
}
|
||||
|
||||
OneComponentAggregationLevelTransferPolicy* clone() const
|
||||
{
|
||||
return new OneComponentAggregationLevelTransferPolicy(*this);
|
||||
}
|
||||
|
||||
const Communication& getCoarseLevelCommunication() const
|
||||
{
|
||||
return *coarseLevelCommunication_;
|
||||
}
|
||||
private:
|
||||
typename Operator::matrix_type::field_type prolongDamp_;
|
||||
std::shared_ptr<AggregatesMap> aggregatesMap_;
|
||||
Criterion criterion_;
|
||||
Communication* communication_;
|
||||
std::shared_ptr<Communication> coarseLevelCommunication_;
|
||||
std::shared_ptr<typename CoarseOperator::matrix_type> coarseLevelMatrix_;
|
||||
};
|
||||
|
||||
template<typename O, typename S, typename C,
|
||||
typename P, std::size_t COMPONENT_INDEX>
|
||||
class BlackoilAmg
|
||||
: public Dune::Preconditioner<typename O::domain_type, typename O::range_type>
|
||||
{
|
||||
public:
|
||||
using Operator = O;
|
||||
using Criterion = C;
|
||||
using Communication = P;
|
||||
using Smoother = S;
|
||||
|
||||
protected:
|
||||
using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
|
||||
using CoarseOperator = typename Detail::ScalarType<Operator>::value;
|
||||
using CoarseSmoother = typename Detail::ScalarType<Smoother>::value;
|
||||
using FineCriterion =
|
||||
typename Detail::OneComponentCriterionType<Criterion,COMPONENT_INDEX>::value;
|
||||
using CoarseCriterion = typename Detail::ScalarType<Criterion>::value;
|
||||
using LevelTransferPolicy =
|
||||
OneComponentAggregationLevelTransferPolicy<Operator,
|
||||
FineCriterion,
|
||||
Communication,
|
||||
COMPONENT_INDEX>;
|
||||
using CoarseSolverPolicy =
|
||||
Detail::OneStepAMGCoarseSolverPolicy<CoarseOperator,
|
||||
CoarseSmoother,
|
||||
CoarseCriterion,
|
||||
LevelTransferPolicy>;
|
||||
using TwoLevelMethod =
|
||||
Dune::Amg::TwoLevelMethod<Operator,
|
||||
CoarseSolverPolicy,
|
||||
Smoother>;
|
||||
public:
|
||||
// define the category
|
||||
enum {
|
||||
//! \brief The category the precondtioner is part of.
|
||||
category = Operator::category
|
||||
};
|
||||
BlackoilAmg(const Operator& fineOperator, const Criterion& criterion,
|
||||
const SmootherArgs& smargs, const Communication& comm)
|
||||
: smoother_(Detail::constructSmoother<Smoother>(fineOperator,smargs,comm)),
|
||||
levelTransferPolicy_(criterion, comm),
|
||||
coarseSolverPolicy_(smargs, criterion),
|
||||
twoLevelMethod_(fineOperator, smoother_,
|
||||
levelTransferPolicy_,
|
||||
coarseSolverPolicy_, 0, 1)
|
||||
{}
|
||||
|
||||
void pre(typename TwoLevelMethod::FineDomainType& x,
|
||||
typename TwoLevelMethod::FineRangeType& b)
|
||||
{
|
||||
twoLevelMethod_.pre(x,b);
|
||||
}
|
||||
|
||||
void post(typename TwoLevelMethod::FineDomainType& x)
|
||||
{
|
||||
twoLevelMethod_.post(x);
|
||||
}
|
||||
|
||||
void apply(typename TwoLevelMethod::FineDomainType& v,
|
||||
const typename TwoLevelMethod::FineRangeType& d)
|
||||
{
|
||||
twoLevelMethod_.apply(v, d);
|
||||
}
|
||||
private:
|
||||
std::shared_ptr<Smoother> smoother_;
|
||||
LevelTransferPolicy levelTransferPolicy_;
|
||||
CoarseSolverPolicy coarseSolverPolicy_;
|
||||
TwoLevelMethod twoLevelMethod_;
|
||||
};
|
||||
|
||||
namespace ISTLUtility
|
||||
{
|
||||
///
|
||||
/// \brief A traits class for selecting the types of the preconditioner.
|
||||
///
|
||||
/// \tparam M The type of the matrix.
|
||||
/// \tparam X The type of the domain of the linear problem.
|
||||
/// \tparam Y The type of the range of the linear problem.
|
||||
/// \tparam P The type of the parallel information.
|
||||
/// \tparam C The type of the coarsening criterion to use.
|
||||
/// \tparam index The pressure index.
|
||||
////
|
||||
template<class M, class X, class Y, class P, class C, std::size_t index>
|
||||
struct BlackoilAmgSelector
|
||||
{
|
||||
using Criterion = C;
|
||||
using Selector = CPRSelector<M,X,Y,P>;
|
||||
using ParallelInformation = typename Selector::ParallelInformation;
|
||||
using Operator = typename Selector::Operator;
|
||||
using Smoother = typename Selector::EllipticPreconditioner;
|
||||
using AMG = BlackoilAmg<Operator,Smoother,Criterion,ParallelInformation,index>;
|
||||
};
|
||||
} // end namespace ISTLUtility
|
||||
} // end namespace Opm
|
||||
#endif
|
||||
@ -271,6 +271,29 @@ createEllipticPreconditionerPointer(const M& Ae, double relax,
|
||||
|
||||
#endif // #if HAVE_MPI
|
||||
|
||||
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 )
|
||||
{
|
||||
// TODO: revise choice of parameters
|
||||
int coarsenTarget=1200;
|
||||
using Criterion = C;
|
||||
Criterion criterion(15, coarsenTarget);
|
||||
criterion.setDebugLevel( 0 ); // no debug information, 1 for printing hierarchy information
|
||||
criterion.setDefaultValuesIsotropic(2);
|
||||
criterion.setNoPostSmoothSteps( 1 );
|
||||
criterion.setNoPreSmoothSteps( 1 );
|
||||
|
||||
// for DUNE 2.2 we also need to pass the smoother args
|
||||
typedef typename AMG::Smoother Smoother;
|
||||
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
|
||||
SmootherArgs smootherArgs;
|
||||
smootherArgs.iterations = 1;
|
||||
smootherArgs.relaxationFactor = relax;
|
||||
|
||||
amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) );
|
||||
}
|
||||
|
||||
/// \brief Creates the elliptic preconditioner (ILU0)
|
||||
/// \param opA The operator representing the matrix of the system.
|
||||
/// \param relax The relaxation parameter for ILU0.
|
||||
@ -292,22 +315,7 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
|
||||
// The coarsening criterion used in the AMG
|
||||
typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
|
||||
|
||||
// TODO: revise choice of parameters
|
||||
int coarsenTarget=1200;
|
||||
Criterion criterion(15,coarsenTarget);
|
||||
criterion.setDebugLevel( 0 ); // no debug information, 1 for printing hierarchy information
|
||||
criterion.setDefaultValuesIsotropic(2);
|
||||
criterion.setNoPostSmoothSteps( 1 );
|
||||
criterion.setNoPreSmoothSteps( 1 );
|
||||
|
||||
// for DUNE 2.2 we also need to pass the smoother args
|
||||
typedef typename AMG::Smoother Smoother;
|
||||
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
|
||||
SmootherArgs smootherArgs;
|
||||
smootherArgs.iterations = 1;
|
||||
smootherArgs.relaxationFactor = relax;
|
||||
|
||||
amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) );
|
||||
createAMGPreconditionerPointer<Criterion>(opA, relax, comm, amgPtr);
|
||||
}
|
||||
|
||||
} // end namespace ISTLUtility
|
||||
|
||||
@ -21,6 +21,7 @@
|
||||
#define OPM_ISTLSOLVER_HEADER_INCLUDED
|
||||
|
||||
#include <opm/autodiff/AdditionalObjectDeleter.hpp>
|
||||
#include <opm/autodiff/BlackoilAmg.hpp>
|
||||
#include <opm/autodiff/CPRPreconditioner.hpp>
|
||||
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
|
||||
#include <opm/autodiff/NewtonIterationUtilities.hpp>
|
||||
@ -407,10 +408,8 @@ namespace Opm
|
||||
if( parameters_.linear_solver_use_amg_ )
|
||||
{
|
||||
typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType;
|
||||
typedef typename CPRSelectorType::AMG AMG;
|
||||
typedef typename CPRSelectorType::Operator MatrixOperator;
|
||||
|
||||
std::unique_ptr< AMG > amg;
|
||||
std::unique_ptr< MatrixOperator > opA;
|
||||
|
||||
if( ! std::is_same< LinearOperator, MatrixOperator > :: value )
|
||||
@ -420,12 +419,34 @@ namespace Opm
|
||||
}
|
||||
|
||||
const double relax = 1.0;
|
||||
if ( ! parameters_.amg_blackoil_system_ )
|
||||
{
|
||||
typedef typename CPRSelectorType::AMG AMG;
|
||||
std::unique_ptr< AMG > amg;
|
||||
|
||||
// Construct preconditioner.
|
||||
constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
|
||||
// Construct preconditioner.
|
||||
constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
|
||||
|
||||
// Solve.
|
||||
solve(linearOperator, x, istlb, *sp, *amg, result);
|
||||
// Solve.
|
||||
solve(linearOperator, x, istlb, *sp, *amg, result);
|
||||
}
|
||||
else
|
||||
{
|
||||
using Matrix = typename MatrixOperator::matrix_type;
|
||||
using CouplingMetric = Dune::Amg::Diagonal<pressureIndex>;
|
||||
using CritBase = Dune::Amg::SymmetricCriterion<Matrix, CouplingMetric>;
|
||||
using Criterion = Dune::Amg::CoarsenCriterion<CritBase>;
|
||||
using AMG = typename ISTLUtility
|
||||
::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureIndex >::AMG;
|
||||
|
||||
std::unique_ptr< AMG > amg;
|
||||
// Construct preconditioner.
|
||||
Criterion crit(15, 2000);
|
||||
constructAMGPrecond<Criterion>( linearOperator, parallelInformation_arg, amg, opA, relax );
|
||||
|
||||
// Solve.
|
||||
solve(linearOperator, x, istlb, *sp, *amg, result);
|
||||
}
|
||||
}
|
||||
else
|
||||
#endif
|
||||
@ -495,6 +516,20 @@ namespace Opm
|
||||
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax, 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
|
||||
{
|
||||
ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax, comm, amg );
|
||||
}
|
||||
|
||||
|
||||
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
|
||||
{
|
||||
ISTLUtility::template createAMGPreconditionerPointer<C>( opA, relax, comm, amg );
|
||||
}
|
||||
/// \brief Solve the system using the given preconditioner and scalar product.
|
||||
template <class Operator, class ScalarProd, class Precond>
|
||||
void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const
|
||||
|
||||
@ -46,6 +46,7 @@ namespace Opm
|
||||
bool require_full_sparsity_pattern_;
|
||||
bool ignoreConvergenceFailure_;
|
||||
bool linear_solver_use_amg_;
|
||||
bool amg_blackoil_system_;
|
||||
|
||||
NewtonIterationBlackoilInterleavedParameters() { reset(); }
|
||||
// read values from parameter class
|
||||
@ -63,6 +64,7 @@ namespace Opm
|
||||
require_full_sparsity_pattern_ = param.getDefault("require_full_sparsity_pattern", require_full_sparsity_pattern_);
|
||||
ignoreConvergenceFailure_ = param.getDefault("linear_solver_ignoreconvergencefailure", ignoreConvergenceFailure_);
|
||||
linear_solver_use_amg_ = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ );
|
||||
amg_blackoil_system_ = param.getDefault("amg_blackoil_system", amg_blackoil_system_);
|
||||
ilu_relaxation_ = param.getDefault("ilu_relaxation", ilu_relaxation_ );
|
||||
ilu_fillin_level_ = param.getDefault("ilu_fillin_level", ilu_fillin_level_ );
|
||||
}
|
||||
@ -70,6 +72,7 @@ namespace Opm
|
||||
// set default values
|
||||
void reset()
|
||||
{
|
||||
amg_blackoil_system_ = false;
|
||||
newton_use_gmres_ = false;
|
||||
linear_solver_reduction_ = 1e-2;
|
||||
linear_solver_maxiter_ = 150;
|
||||
|
||||
327
tests/test_blackoil_amg.cpp
Normal file
327
tests/test_blackoil_amg.cpp
Normal file
@ -0,0 +1,327 @@
|
||||
/*
|
||||
Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
|
||||
|
||||
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>
|
||||
#if HAVE_DYNAMIC_BOOST_TEST && HAVE_MPI
|
||||
#define BOOST_TEST_DYN_LINK
|
||||
#define BOOST_TEST_MODULE BlackoilAmgTest
|
||||
#define BOOST_TEST_NO_MAIN
|
||||
#include <boost/test/unit_test.hpp>
|
||||
#include <opm/autodiff/BlackoilAmg.hpp>
|
||||
|
||||
#include <dune/common/parallel/mpihelper.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
#include <dune/common/unused.hh>
|
||||
#include <dune/common/parallel/indexset.hh>
|
||||
#include <dune/common/parallel/plocalindex.hh>
|
||||
#include <dune/common/parallel/collectivecommunication.hh>
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/istl/owneroverlapcopy.hh>
|
||||
|
||||
class MPIError {
|
||||
public:
|
||||
/** @brief Constructor. */
|
||||
MPIError(std::string s, int e) : errorstring(s), errorcode(e){}
|
||||
/** @brief The error string. */
|
||||
std::string errorstring;
|
||||
/** @brief The mpi error code. */
|
||||
int errorcode;
|
||||
};
|
||||
|
||||
void MPI_err_handler(MPI_Comm *, int *err_code, ...){
|
||||
char *err_string=new char[MPI_MAX_ERROR_STRING];
|
||||
int err_length;
|
||||
MPI_Error_string(*err_code, err_string, &err_length);
|
||||
std::string s(err_string, err_length);
|
||||
std::cerr << "An MPI Error ocurred:"<<std::endl<<s<<std::endl;
|
||||
delete[] err_string;
|
||||
throw MPIError(s, *err_code);
|
||||
}
|
||||
|
||||
typedef Dune::OwnerOverlapCopyAttributeSet GridAttributes;
|
||||
typedef GridAttributes::AttributeSet GridFlag;
|
||||
typedef Dune::ParallelLocalIndex<GridFlag> LocalIndex;
|
||||
|
||||
template<class M, class G, class L, int n>
|
||||
void setupPattern(int N, M& mat, Dune::ParallelIndexSet<G,L,n>& indices, int overlapStart, int overlapEnd,
|
||||
int start, int end);
|
||||
|
||||
template<class M>
|
||||
void fillValues(int N, M& mat, int overlapStart, int overlapEnd, int start, int end);
|
||||
|
||||
|
||||
template<class M, class G, class L, class C, int n>
|
||||
M setupAnisotropic2d(int N, Dune::ParallelIndexSet<G,L,n>& indices,
|
||||
const Dune::CollectiveCommunication<C>& p, int *nout, typename M::block_type::value_type eps=1.0);
|
||||
|
||||
|
||||
template<class M, class G, class L, int s>
|
||||
void setupPattern(int N, M& mat, Dune::ParallelIndexSet<G,L,s>& indices, int overlapStart, int overlapEnd,
|
||||
int start, int end)
|
||||
{
|
||||
int n = overlapEnd - overlapStart;
|
||||
|
||||
typename M::CreateIterator iter = mat.createbegin();
|
||||
indices.beginResize();
|
||||
|
||||
for(int j=0; j < N; j++)
|
||||
for(int i=overlapStart; i < overlapEnd; i++, ++iter) {
|
||||
int global = j*N+i;
|
||||
GridFlag flag = GridAttributes::owner;
|
||||
bool isPublic = false;
|
||||
|
||||
if((i<start && i > 0) || (i>= end && i < N-1))
|
||||
flag=GridAttributes::copy;
|
||||
|
||||
if(i<start+1 || i>= end-1) {
|
||||
isPublic = true;
|
||||
indices.add(global, LocalIndex(iter.index(), flag, isPublic));
|
||||
}
|
||||
|
||||
|
||||
iter.insert(iter.index());
|
||||
|
||||
// i direction
|
||||
if(i > overlapStart )
|
||||
// We have a left neighbour
|
||||
iter.insert(iter.index()-1);
|
||||
|
||||
if(i < overlapEnd-1)
|
||||
// We have a rigt neighbour
|
||||
iter.insert(iter.index()+1);
|
||||
|
||||
// j direction
|
||||
// Overlap is a dirichlet border, discard neighbours
|
||||
if(flag != GridAttributes::copy) {
|
||||
if(j>0)
|
||||
// lower neighbour
|
||||
iter.insert(iter.index()-n);
|
||||
if(j<N-1)
|
||||
// upper neighbour
|
||||
iter.insert(iter.index()+n);
|
||||
}
|
||||
}
|
||||
indices.endResize();
|
||||
}
|
||||
|
||||
template<class M, class T>
|
||||
void fillValues(int N, M& mat, int overlapStart, int overlapEnd, int start, int end, T eps)
|
||||
{
|
||||
DUNE_UNUSED_PARAMETER(N);
|
||||
typedef typename M::block_type Block;
|
||||
Block dval(0), bone(0), bmone(0), beps(0);
|
||||
|
||||
for(typename Block::RowIterator b = dval.begin(); b != dval.end(); ++b)
|
||||
b->operator[](b.index())=2.0+2.0*eps;
|
||||
|
||||
for(typename Block::RowIterator b=bone.begin(); b != bone.end(); ++b)
|
||||
b->operator[](b.index())=1.0;
|
||||
|
||||
for(typename Block::RowIterator b=bmone.begin(); b != bmone.end(); ++b)
|
||||
b->operator[](b.index())=-1.0;
|
||||
|
||||
for(typename Block::RowIterator b=beps.begin(); b != beps.end(); ++b)
|
||||
b->operator[](b.index())=-eps;
|
||||
|
||||
int n = overlapEnd-overlapStart;
|
||||
typedef typename M::ColIterator ColIterator;
|
||||
typedef typename M::RowIterator RowIterator;
|
||||
|
||||
for (RowIterator i = mat.begin(); i != mat.end(); ++i) {
|
||||
// calculate coordinate
|
||||
int y = i.index() / n;
|
||||
int x = overlapStart + i.index() - y * n;
|
||||
|
||||
ColIterator endi = (*i).end();
|
||||
|
||||
if(x<start || x >= end) { // || x==0 || x==N-1 || y==0 || y==N-1){
|
||||
// overlap node is dirichlet border
|
||||
ColIterator j = (*i).begin();
|
||||
|
||||
for(; j.index() < i.index(); ++j)
|
||||
*j=0;
|
||||
|
||||
*j = bone;
|
||||
|
||||
for(++j; j != endi; ++j)
|
||||
*j=0;
|
||||
}else{
|
||||
for(ColIterator j = (*i).begin(); j != endi; ++j)
|
||||
if(j.index() == i.index())
|
||||
*j=dval;
|
||||
else if(j.index()+1==i.index() || j.index()-1==i.index())
|
||||
*j=beps;
|
||||
else
|
||||
*j=bmone;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<class V, class G, class L, int s>
|
||||
void setBoundary(V& lhs, V& rhs, const G& n, Dune::ParallelIndexSet<G,L,s>& indices)
|
||||
{
|
||||
typedef typename Dune::ParallelIndexSet<G,L,s>::const_iterator Iter;
|
||||
for(Iter i=indices.begin(); i != indices.end(); ++i) {
|
||||
G x = i->global()/n;
|
||||
G y = i->global()%n;
|
||||
|
||||
if(x==0 || y ==0 || x==n-1 || y==n-1) {
|
||||
//double h = 1.0 / ((double) (n-1));
|
||||
lhs[i->local()]=rhs[i->local()]=0; //((double)x)*((double)y)*h*h;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<class V, class G>
|
||||
void setBoundary(V& lhs, V& rhs, const G& N)
|
||||
{
|
||||
typedef typename V::block_type Block;
|
||||
typedef typename Block::value_type T;
|
||||
for(int j=0; j < N; ++j)
|
||||
for(int i=0; i < N; i++)
|
||||
if(i==0 || j ==0 || i==N-1 || j==N-1) {
|
||||
T h = 1.0 / ((double) (N-1));
|
||||
T x, y;
|
||||
if(i==N-1)
|
||||
x=1;
|
||||
else
|
||||
x = ((T) i)*h;
|
||||
if(j==N-1)
|
||||
y = 1;
|
||||
else
|
||||
y = ((T) j)*h;
|
||||
|
||||
lhs[j*N+i]=rhs[j*N+i]=0; //x*y;
|
||||
}
|
||||
}
|
||||
|
||||
template<class M, class G, class L, class C, int s>
|
||||
M setupAnisotropic2d(int N, Dune::ParallelIndexSet<G,L,s>& indices, const Dune::CollectiveCommunication<C>& p, int *nout, typename M::block_type::value_type eps)
|
||||
{
|
||||
int procs=p.size(), rank=p.rank();
|
||||
|
||||
typedef M BCRSMat;
|
||||
|
||||
// calculate size of local matrix in the distributed direction
|
||||
int start, end, overlapStart, overlapEnd;
|
||||
|
||||
int n = N/procs; // number of unknowns per process
|
||||
int bigger = N%procs; // number of process with n+1 unknows
|
||||
|
||||
// Compute owner region
|
||||
if(rank<bigger) {
|
||||
start = rank*(n+1);
|
||||
end = start+(n+1);
|
||||
}else{
|
||||
start = bigger*(n+1) + (rank-bigger) * n;
|
||||
end = start+n;
|
||||
}
|
||||
|
||||
// Compute overlap region
|
||||
if(start>0)
|
||||
overlapStart = start - 1;
|
||||
else
|
||||
overlapStart = start;
|
||||
|
||||
if(end<N)
|
||||
overlapEnd = end + 1;
|
||||
else
|
||||
overlapEnd = end;
|
||||
|
||||
int noKnown = overlapEnd-overlapStart;
|
||||
|
||||
*nout = noKnown;
|
||||
|
||||
BCRSMat mat(noKnown*N, noKnown*N, noKnown*N*5, BCRSMat::row_wise);
|
||||
|
||||
setupPattern(N, mat, indices, overlapStart, overlapEnd, start, end);
|
||||
fillValues(N, mat, overlapStart, overlapEnd, start, end, eps);
|
||||
|
||||
// Dune::printmatrix(std::cout,mat,"aniso 2d","row",9,1);
|
||||
|
||||
return mat;
|
||||
}
|
||||
|
||||
//BOOST_AUTO_TEST_CASE(runBlackoilAmgLaplace)
|
||||
void runBlackoilAmgLaplace()
|
||||
{
|
||||
const int BS=2, N=100;
|
||||
typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
|
||||
typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
|
||||
typedef Dune::FieldVector<double,BS> VectorBlock;
|
||||
typedef Dune::BlockVector<VectorBlock> Vector;
|
||||
typedef int GlobalId;
|
||||
typedef Dune::OwnerOverlapCopyCommunication<GlobalId> Communication;
|
||||
typedef Dune::OverlappingSchwarzOperator<BCRSMat,Vector,Vector,Communication> Operator;
|
||||
int argc;
|
||||
char** argv;
|
||||
const auto& ccomm = Dune::MPIHelper::instance(argc, argv).getCollectiveCommunication();
|
||||
|
||||
Communication comm(ccomm);
|
||||
int n=0;
|
||||
BCRSMat mat = setupAnisotropic2d<BCRSMat>(N, comm.indexSet(), comm.communicator(), &n, 1);
|
||||
|
||||
|
||||
comm.remoteIndices().template rebuild<false>();
|
||||
|
||||
Vector b(mat.N()), x(mat.M());
|
||||
|
||||
b=0;
|
||||
x=100;
|
||||
setBoundary(x, b, N, comm.indexSet());
|
||||
|
||||
Operator fop(mat, comm);
|
||||
typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
|
||||
Criterion;
|
||||
typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother;
|
||||
//typedef Dune::SeqJac<BCRSMat,Vector,Vector> Smoother;
|
||||
//typedef Dune::SeqILU0<BCRSMat,Vector,Vector> Smoother;
|
||||
//typedef Dune::SeqILUn<BCRSMat,Vector,Vector> Smoother;
|
||||
typedef Dune::BlockPreconditioner<Vector,Vector,Communication,Smoother> ParSmoother;
|
||||
typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
|
||||
|
||||
Dune::OverlappingSchwarzScalarProduct<Vector,Communication> sp(comm);
|
||||
|
||||
Dune::InverseOperatorResult r;
|
||||
SmootherArgs smootherArgs;
|
||||
Criterion criterion;
|
||||
|
||||
smootherArgs.iterations = 1;
|
||||
|
||||
Opm::BlackoilAmg<Operator,ParSmoother,Criterion,Communication,0> amg(fop, criterion,
|
||||
smootherArgs,
|
||||
comm);
|
||||
Dune::CGSolver<Vector> amgCG(fop, sp, amg, 10e-8, 300, (ccomm.rank()==0) ? 2 : 0);
|
||||
amgCG.apply(x,b,r);
|
||||
|
||||
}
|
||||
|
||||
bool init_unit_test_func()
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
const auto& helper = Dune::MPIHelper::instance(argc, argv);
|
||||
boost::unit_test::unit_test_main(&init_unit_test_func,
|
||||
argc, argv);
|
||||
}
|
||||
#else
|
||||
int main () { return 0; }
|
||||
#endif // #if HAVE_DYNAMIC_BOOST_TEST
|
||||
Loading…
Reference in New Issue
Block a user