From f5d81513da1678ef7b32d332b096cb5e15efe94b Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Thu, 7 Sep 2017 09:15:39 +0200 Subject: [PATCH] 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. --- CMakeLists_files.cmake | 2 + opm/autodiff/BlackoilAmg.hpp | 727 ++++++++++++++++++ opm/autodiff/CPRPreconditioner.hpp | 40 +- opm/autodiff/ISTLSolver.hpp | 47 +- .../NewtonIterationBlackoilInterleaved.hpp | 3 + tests/test_blackoil_amg.cpp | 327 ++++++++ 6 files changed, 1124 insertions(+), 22 deletions(-) create mode 100644 opm/autodiff/BlackoilAmg.hpp create mode 100644 tests/test_blackoil_amg.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index c62b18c85..5f53b3021 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -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 diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp new file mode 100644 index 000000000..0f2c1ab72 --- /dev/null +++ b/opm/autodiff/BlackoilAmg.hpp @@ -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 . +*/ +#ifndef OPM_AMG_HEADER_INCLUDED +#define OPM_AMG_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +namespace Dune +{ +namespace Amg +{ +template +class UnSymmetricCriterion; +} +} + +namespace Dune +{ + +template +class MatrixBlock; + +} + +namespace Opm +{ + +namespace Detail +{ + +template +struct ScalarType +{ +}; + +template +struct ScalarType > +{ + typedef Dune::FieldVector value; +}; + +template +struct ScalarType > +{ + typedef Dune::FieldMatrix value; +}; + +template +struct ScalarType > +{ + typedef Dune::MatrixBlock value; +}; + +template +struct ScalarType > +{ + using ScalarBlock = typename ScalarType::value; + using ScalarAllocator = typename Allocator::template rebind::other; + typedef Dune::BCRSMatrix value; +}; + +template +struct ScalarType > +{ + using ScalarBlock = typename ScalarType::value; + using ScalarAllocator = typename Allocator::template rebind::other; + typedef Dune::BlockVector value; +}; + +template +struct ScalarType > +{ + typedef Dune::SeqScalarProduct::value> value; +}; + +#define ComposeScalarTypeForSeqPrecond(PREC) \ + template \ + struct ScalarType > \ + { \ + typedef PREC::value, \ + typename ScalarType::value, \ + typename ScalarType::value, \ + l> value; \ + }; + +ComposeScalarTypeForSeqPrecond(Dune::SeqJac); +ComposeScalarTypeForSeqPrecond(Dune::SeqSOR); +ComposeScalarTypeForSeqPrecond(Dune::SeqSSOR); +ComposeScalarTypeForSeqPrecond(Dune::SeqGS); +ComposeScalarTypeForSeqPrecond(Dune::SeqILU0); +ComposeScalarTypeForSeqPrecond(Dune::SeqILUn); + +#undef ComposeScalarTypeForSeqPrecond + +template +struct ScalarType > +{ + typedef Dune::Richardson::value, + typename ScalarType::value> value; +}; + +template +struct ScalarType > +{ + typedef Dune::OverlappingSchwarzOperator::value, + typename ScalarType::value, + typename ScalarType::value, + C> value; +}; + +template +struct ScalarType > +{ + typedef Dune::MatrixAdapter::value, + typename ScalarType::value, + typename ScalarType::value> value; +}; + +template +struct ScalarType > +{ + typedef Dune::OverlappingSchwarzScalarProduct::value, + C> value; +}; + +template +struct ScalarType > +{ + typedef Dune::NonoverlappingSchwarzScalarProduct::value, + C> value; +}; + +template +struct ScalarType > +{ + typedef Dune::BlockPreconditioner::value, + typename ScalarType::value, + C, + typename ScalarType::value> value; +}; + +template +struct ScalarType > +{ + typedef ParallelOverlappingILU0::value, + typename ScalarType::value, + typename ScalarType::value, + C> value; +}; + +template +struct ScalarType,N> > > +{ + using value = Dune::Amg::CoarsenCriterion::value>, Dune::Amg::FirstDiagonal> >; +}; + +template +struct ScalarType,N> > > +{ + using value = Dune::Amg::CoarsenCriterion::value>, Dune::Amg::FirstDiagonal> >; +}; + +template +struct OneComponentCriterionType +{}; + +template +struct OneComponentCriterionType,N> >,COMPONENT_INDEX> +{ + using value = Dune::Amg::CoarsenCriterion, Dune::Amg::Diagonal > >; +}; + +template +struct OneComponentCriterionType,N> >,COMPONENT_INDEX> +{ + using value = Dune::Amg::CoarsenCriterion, Dune::Amg::Diagonal > >; +}; + +template +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 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::Arguments SmootherArgs; + /** @brief The type of the AMG construct on the coarse level.*/ + typedef Dune::Amg::AMG 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 + { + 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; + auto sp = Chooser::construct(comm_); + Dune::BiCGSTABSolver solver(const_cast(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 + CoarseLevelSolver* createCoarseLevelSolver(LTP& transferPolicy) + { + coarseOperator_=transferPolicy.getCoarseLevelOperator(); + const LevelTransferPolicy& transfer = + reinterpret_cast(transferPolicy); + AMGInverseOperator* inv = new AMGInverseOperator(*coarseOperator_, + criterion_, + smootherArgs_, + transfer.getCoarseLevelCommunication()); + + return inv; //std::shared_ptr >(inv); + + } + +private: + /** @brief The coarse level operator. */ + std::shared_ptr coarseOperator_; + /** @brief The arguments used to construct the smoother. */ + SmootherArgs smootherArgs_; + /** @brief The coarsening criterion. */ + Criterion criterion_; +}; + +template +Smoother* constructSmoother(const Operator& op, + const typename Dune::Amg::SmootherTraits::Arguments& smargs, + const Communication& comm) +{ + typename Dune::Amg::ConstructionTraits::Arguments args; + args.setMatrix(op.getmat()); + args.setComm(comm); + args.setArgs(smargs); + return Dune::Amg::ConstructionTraits::construct(args); +} + +template +const Dune::Amg::OverlapVertex* +buildOverlapVertices(const G& graph, const C& pinfo, + Dune::Amg::AggregatesMap& 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; + + 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 +void buildCoarseSparseMatrix(M& coarseMatrix, G& fineGraph, + const V& visitedMap, + const C& pinfo, + Dune::Amg::AggregatesMap& aggregates, + const S& overlap) +{ + using OverlapVertex = Dune::Amg ::OverlapVertex; + 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::UNAGGREGATED; +#endif + const auto ISOLATED = Dune::Amg::AggregatesMap::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 sparsityBuilder(coarseMatrix); + + Dune::Amg::ConnectivityConstructor::examine(fineGraph, visitedMap, pinfo, + aggregates, overlap, + overlapVertices, + overlapVertices+count, + sparsityBuilder); + delete[] overlapVertices; +} + +template +void buildCoarseSparseMatrix(M& coarseMatrix, G& fineGraph, const V& visitedMap, + const Dune::Amg::SequentialInformation& pinfo, + Dune::Amg::AggregatesMap& 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::UNAGGREGATED; +#endif + const auto ISOLATED = Dune::Amg::AggregatesMap::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 sparsityBuilder(coarseMatrix); + + Dune::Amg::ConnectivityConstructor + ::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 OneComponentAggregationLevelTransferPolicy + : public Dune::Amg::LevelTransferPolicy::value> +{ + typedef Dune::Amg::AggregatesMap AggregatesMap; +public: + using CoarseOperator = typename Detail::ScalarType::value; + typedef Dune::Amg::LevelTransferPolicy FatherType; + typedef Communication ParallelInformation; + +public: + OneComponentAggregationLevelTransferPolicy(const Criterion& crit, const Communication& comm) + : criterion_(crit), communication_(&const_cast(comm)) + {} + + void createCoarseLevelSystem(const Operator& fineOperator) + { + prolongDamp_ = 1; + typedef Dune::Amg::PropertiesGraphCreator GraphCreator; + typedef typename GraphCreator::PropertiesGraph PropertiesGraph; + typedef typename GraphCreator::GraphTuple GraphTuple; + + typedef typename PropertiesGraph::VertexDescriptor Vertex; + + std::vector excluded(fineOperator.getmat().N(), false); + + using OverlapFlags = Dune::NegateSet; + 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="<. +*/ +#include +#if HAVE_DYNAMIC_BOOST_TEST && HAVE_MPI +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE BlackoilAmgTest +#define BOOST_TEST_NO_MAIN +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +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:"< LocalIndex; + +template +void setupPattern(int N, M& mat, Dune::ParallelIndexSet& indices, int overlapStart, int overlapEnd, + int start, int end); + +template +void fillValues(int N, M& mat, int overlapStart, int overlapEnd, int start, int end); + + +template +M setupAnisotropic2d(int N, Dune::ParallelIndexSet& indices, + const Dune::CollectiveCommunication& p, int *nout, typename M::block_type::value_type eps=1.0); + + +template +void setupPattern(int N, M& mat, Dune::ParallelIndexSet& 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 0) || (i>= end && i < N-1)) + flag=GridAttributes::copy; + + if(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 +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= 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 +void setBoundary(V& lhs, V& rhs, const G& n, Dune::ParallelIndexSet& indices) +{ + typedef typename Dune::ParallelIndexSet::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 +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 +M setupAnisotropic2d(int N, Dune::ParallelIndexSet& indices, const Dune::CollectiveCommunication& 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(rank0) + overlapStart = start - 1; + else + overlapStart = start; + + if(end MatrixBlock; + typedef Dune::BCRSMatrix BCRSMat; + typedef Dune::FieldVector VectorBlock; + typedef Dune::BlockVector Vector; + typedef int GlobalId; + typedef Dune::OwnerOverlapCopyCommunication Communication; + typedef Dune::OverlappingSchwarzOperator Operator; + int argc; + char** argv; + const auto& ccomm = Dune::MPIHelper::instance(argc, argv).getCollectiveCommunication(); + + Communication comm(ccomm); + int n=0; + BCRSMat mat = setupAnisotropic2d(N, comm.indexSet(), comm.communicator(), &n, 1); + + + comm.remoteIndices().template rebuild(); + + 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 > + Criterion; + typedef Dune::SeqSSOR Smoother; + //typedef Dune::SeqJac Smoother; + //typedef Dune::SeqILU0 Smoother; + //typedef Dune::SeqILUn Smoother; + typedef Dune::BlockPreconditioner ParSmoother; + typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; + + Dune::OverlappingSchwarzScalarProduct sp(comm); + + Dune::InverseOperatorResult r; + SmootherArgs smootherArgs; + Criterion criterion; + + smootherArgs.iterations = 1; + + Opm::BlackoilAmg amg(fop, criterion, + smootherArgs, + comm); + Dune::CGSolver 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