diff --git a/opm/core/linalg/LinearSolverFactory.cpp b/opm/core/linalg/LinearSolverFactory.cpp index bd517b668..207c4ea54 100644 --- a/opm/core/linalg/LinearSolverFactory.cpp +++ b/opm/core/linalg/LinearSolverFactory.cpp @@ -97,9 +97,10 @@ namespace Opm const int* ja, const double* sa, const double* rhs, - double* solution) const + double* solution, + const boost::any& add) const { - return solver_->solve(size, nonzeros, ia, ja, sa, rhs, solution); + return solver_->solve(size, nonzeros, ia, ja, sa, rhs, solution, add); } void LinearSolverFactory::setTolerance(const double tol) diff --git a/opm/core/linalg/LinearSolverFactory.hpp b/opm/core/linalg/LinearSolverFactory.hpp index 6cc2354b6..f27140bc0 100644 --- a/opm/core/linalg/LinearSolverFactory.hpp +++ b/opm/core/linalg/LinearSolverFactory.hpp @@ -74,7 +74,8 @@ namespace Opm const int* ja, const double* sa, const double* rhs, - double* solution) const; + double* solution, + const boost::any& add=boost::any()) const; /// Set tolerance for the linear solver. /// \param[in] tol tolerance value diff --git a/opm/core/linalg/LinearSolverInterface.hpp b/opm/core/linalg/LinearSolverInterface.hpp index 559eaaec1..ba60e7b71 100644 --- a/opm/core/linalg/LinearSolverInterface.hpp +++ b/opm/core/linalg/LinearSolverInterface.hpp @@ -20,6 +20,7 @@ #ifndef OPM_LINEARSOLVERINTERFACE_HEADER_INCLUDED #define OPM_LINEARSOLVERINTERFACE_HEADER_INCLUDED +#include struct CSRMatrix; @@ -69,7 +70,8 @@ namespace Opm const int* ja, const double* sa, const double* rhs, - double* solution) const = 0; + double* solution, + const boost::any& add=boost::any()) const = 0; /// Set tolerance for the linear solver. /// \param[in] tol tolerance value diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index 181b9ed2f..822782904 100644 --- a/opm/core/linalg/LinearSolverIstl.cpp +++ b/opm/core/linalg/LinearSolverIstl.cpp @@ -23,6 +23,7 @@ #endif #include +#include // Silence compatibility warning from DUNE headers since we don't use // the deprecated member anyway (in this compilation unit) @@ -35,7 +36,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -134,7 +137,8 @@ namespace Opm const int* ja, const double* sa, const double* rhs, - double* solution) const + double* solution, + const boost::any& comm) const { // Build Istl structures from input. // System matrix @@ -173,10 +177,26 @@ namespace Opm if (maxit == 0) { maxit = 5000; } - Dune::SeqScalarProduct sp; - Dune::Amg::SequentialInformation comm; - Operator opA(A); - return solveSystem(opA, x, solution, b, sp, comm, maxit); +#if HAVE_MPI + if(comm.type()==typeid(ParallelISTLInformation)) + { + typedef Dune::OwnerOverlapCopyCommunication Comm; + const ParallelISTLInformation& info = boost::any_cast(comm); + Comm istlComm(info.communicator()); + info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices()); + Dune::OverlappingSchwarzOperator + opA(A, istlComm); + Dune::OverlappingSchwarzScalarProduct sp(istlComm); + return solveSystem(opA, x, solution, b, sp, istlComm, maxit); + } + else +#endif + { + Dune::SeqScalarProduct sp; + Dune::Amg::SequentialInformation comm; + Operator opA(A); + return solveSystem(opA, x, solution, b, sp, comm, maxit); + } } template @@ -237,24 +257,44 @@ template namespace { + template + struct SmootherChooser + { + typedef P Type; + }; + +#if HAVE_MPI template + struct SmootherChooser > + { + typedef Dune::OwnerOverlapCopyCommunication Comm; + typedef Dune::BlockPreconditioner + Type; + }; +#endif + + template struct PreconditionerTraits { - typedef std::shared_ptr

PointerType; + typedef typename SmootherChooser::Type SmootherType; + typedef std::shared_ptr PointerType; }; template - typename PreconditionerTraits::PointerType + typename PreconditionerTraits::PointerType makePreconditioner(O& opA, double relax, const C& comm, int iterations=1) { - typename Dune::Amg::SmootherTraits

::Arguments args; - typename Dune::Amg::ConstructionTraits

::Arguments cargs; + typedef typename SmootherChooser::Type SmootherType; + typedef typename PreconditionerTraits::PointerType PointerType; + typename Dune::Amg::SmootherTraits::Arguments args; + typename Dune::Amg::ConstructionTraits::Arguments cargs; cargs.setMatrix(opA.getmat()); args.iterations=iterations; args.relaxationFactor=relax; cargs.setArgs(args); cargs.setComm(comm); - return std::shared_ptr

(Dune::Amg::ConstructionTraits

::construct(cargs)); + return PointerType(Dune::Amg::ConstructionTraits::construct(cargs)); } template @@ -322,19 +362,20 @@ template #endif #if SMOOTHER_ILU - typedef Dune::SeqILU0 Smoother; + typedef Dune::SeqILU0 SeqSmoother; #else - typedef Dune::SeqSOR Smoother; + typedef Dune::SeqSOR SeqSmoother; #endif + typedef typename SmootherChooser::Type Smoother; typedef Dune::Amg::CoarsenCriterion Criterion; - typedef Dune::Amg::AMG Precond; + typedef Dune::Amg::AMG Precond; // Construct preconditioner. Criterion criterion; - Precond::SmootherArgs smootherArgs; + typename Precond::SmootherArgs smootherArgs; setUpCriterion(criterion, linsolver_prolongate_factor, verbosity, linsolver_smooth_steps); - Precond precond(opA, criterion, smootherArgs); + Precond precond(opA, criterion, smootherArgs, comm); // Construct linear solver. Dune::CGSolver linsolve(opA, sp, precond, tolerance, maxit, verbosity); @@ -359,6 +400,7 @@ template double linsolver_prolongate_factor, int linsolver_smooth_steps) { // Solve with AMG solver. + Dune::MatrixAdapter sOpA(opA.getmat()); #if FIRST_DIAGONAL typedef Dune::Amg::FirstDiagonal CouplingMetric; @@ -385,10 +427,10 @@ template Criterion criterion; setUpCriterion(criterion, linsolver_prolongate_factor, verbosity, linsolver_smooth_steps); - Precond precond(opA, criterion, smootherArgs); + Precond precond(sOpA, criterion, smootherArgs); // Construct linear solver. - Dune::GeneralizedPCGSolver linsolve(opA, sp, precond, tolerance, maxit, verbosity); + Dune::GeneralizedPCGSolver linsolve(sOpA, precond, tolerance, maxit, verbosity); // Solve system. Dune::InverseOperatorResult result; @@ -408,6 +450,8 @@ template double linsolver_prolongate_factor) { // Solve with AMG solver. + typedef Dune::MatrixAdapter Operator; + Operator sOpA(opA.getmat()); #if FIRST_DIAGONAL typedef Dune::Amg::FirstDiagonal CouplingMetric; @@ -433,10 +477,10 @@ template parms.setNoPreSmoothSteps(smooth_steps); parms.setNoPostSmoothSteps(smooth_steps); parms.setProlongationDampingFactor(linsolver_prolongate_factor); - Precond precond(opA, criterion, parms); + Precond precond(sOpA, criterion, parms); // Construct linear solver. - Dune::GeneralizedPCGSolver linsolve(opA, sp, precond, tolerance, maxit, verbosity); + Dune::GeneralizedPCGSolver linsolve(sOpA, precond, tolerance, maxit, verbosity); // Solve system. Dune::InverseOperatorResult result; diff --git a/opm/core/linalg/LinearSolverIstl.hpp b/opm/core/linalg/LinearSolverIstl.hpp index df7b681e9..de32783f3 100644 --- a/opm/core/linalg/LinearSolverIstl.hpp +++ b/opm/core/linalg/LinearSolverIstl.hpp @@ -24,12 +24,11 @@ #include #include #include - +#include namespace Opm { - /// Concrete class encapsulating some dune-istl linear solvers. class LinearSolverIstl : public LinearSolverInterface { @@ -75,7 +74,8 @@ namespace Opm const int* ja, const double* sa, const double* rhs, - double* solution) const; + double* solution, + const boost::any& comm=boost::any()) const; /// Set tolerance for the residual in dune istl linear solver. /// \param[in] tol tolerance value diff --git a/opm/core/linalg/LinearSolverUmfpack.cpp b/opm/core/linalg/LinearSolverUmfpack.cpp index 8980d4b19..e7067b69f 100644 --- a/opm/core/linalg/LinearSolverUmfpack.cpp +++ b/opm/core/linalg/LinearSolverUmfpack.cpp @@ -46,7 +46,8 @@ namespace Opm const int* ja, const double* sa, const double* rhs, - double* solution) const + double* solution, + const boost::any&) const { CSRMatrix A = { (size_t)size, diff --git a/opm/core/linalg/LinearSolverUmfpack.hpp b/opm/core/linalg/LinearSolverUmfpack.hpp index f2562a156..7126bc4a6 100644 --- a/opm/core/linalg/LinearSolverUmfpack.hpp +++ b/opm/core/linalg/LinearSolverUmfpack.hpp @@ -55,7 +55,8 @@ namespace Opm const int* ja, const double* sa, const double* rhs, - double* solution) const; + double* solution, + const boost::any& add=boost::any()) const; /// Set tolerance for the linear solver. /// \param[in] tol tolerance value diff --git a/opm/core/linalg/ParallelIstlInformation.hpp b/opm/core/linalg/ParallelIstlInformation.hpp new file mode 100644 index 000000000..38795774c --- /dev/null +++ b/opm/core/linalg/ParallelIstlInformation.hpp @@ -0,0 +1,160 @@ +/* + Copyright 2014 Dr. Markus 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_PARALLELISTLINFORMTION_HEADER_INCLUDED +#define OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED + +#if HAVE_MPI && HAVE_DUNE_ISTL + +#include "mpi.h" +#include +#include +#include +#include + +#include + +namespace Opm +{ + +/// \brief Class that encapsulates the parallelization information needed by the +/// ISTL solvers. +class ParallelISTLInformation +{ +public: + /// \brief The type of the parallel index set used. + typedef Dune::OwnerOverlapCopyCommunication::ParallelIndexSet ParallelIndexSet; + /// \brief The type of the remote indices information used. + typedef Dune::OwnerOverlapCopyCommunication::RemoteIndices RemoteIndices; + + /// \brief Constructs an empty parallel information object using MPI_COMM_WORLD + ParallelISTLInformation() + : indexSet_(new ParallelIndexSet), + remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, MPI_COMM_WORLD)), + communicator_(MPI_COMM_WORLD) + {} + /// \brief Constructs an empty parallel information object using a communicator. + /// \param communicator The communicator to use. + ParallelISTLInformation(MPI_Comm communicator) + : indexSet_(new ParallelIndexSet), + remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, communicator)), + communicator_(communicator) + {} + /// \brief Constructs a parallel information object from the specified information. + /// \param indexSet The parallel index set to use. + /// \param remoteIndices The remote indices information to use. + /// \param communicator The communicator to use. + ParallelISTLInformation(const std::shared_ptr& indexSet, + const std::shared_ptr& remoteIndices, + MPI_Comm communicator) + : indexSet_(indexSet), remoteIndices_(remoteIndices), communicator_(communicator) + {} + /// \brief Copy constructor. + /// + /// The information will be shared by the the two objects. + ParallelISTLInformation(const ParallelISTLInformation& other) + : indexSet_(other.indexSet_), remoteIndices_(other.remoteIndices_), + communicator_(other.communicator_) + {} + /// \brief Get a pointer to the underlying index set. + std::shared_ptr indexSet() const + { + return indexSet_; + } + /// \brief Get a pointer to the remote indices information. + std::shared_ptr remoteIndices() const + { + return remoteIndices_; + } + /// \brief Get the MPI communicator that we use. + MPI_Comm communicator() const + { + return communicator_; + } + /// \brief Copy the information stored to the specified objects. + /// \param[out] indexSet The object to store the index set in. + /// \param[out] remoteIndices The object to store the remote indices information in. + void copyValuesTo(ParallelIndexSet& indexSet, RemoteIndices& remoteIndices) const + { + indexSet.beginResize(); + IndexSetInserter inserter(indexSet); + std::for_each(indexSet_->begin(), indexSet_->end(), inserter); + indexSet.endResize(); + remoteIndices.rebuild(); + } + /// \brief Communcate the dofs owned by us to the other process. + /// + /// Afterwards all associated dofs will contain the same data. + template + void copyOwnerToAll (const T& source, T& dest) const + { + typedef Dune::Combine,Dune::EnumItem,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> OwnerOverlapSet; + typedef Dune::Combine,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> AllSet; + OwnerOverlapSet sourceFlags; + AllSet destFlags; + Dune::Interface interface(communicator_); + if(!remoteIndices_->isSynced()) + remoteIndices_->rebuild(); + interface.build(*remoteIndices_,sourceFlags,destFlags); + Dune::BufferedCommunicator communicator; + communicator.template build(interface); + communicator.template forward >(source,dest); + communicator.free(); + } +private: + /** \brief gather/scatter callback for communcation */ + template + struct CopyGatherScatter + { + typedef typename Dune::CommPolicy::IndexedType V; + + static V gather(const T& a, std::size_t i) + { + return a[i]; + } + + static void scatter(T& a, V v, std::size_t i) + { + a[i] = v; + } + }; + template + class IndexSetInserter + { + public: + typedef T ParallelIndexSet; + + IndexSetInserter(ParallelIndexSet& indexSet) + : indexSet_(&indexSet) + {} + void operator()(const typename ParallelIndexSet::IndexPair& pair) + { + indexSet_->add(pair.global(), pair.local()); + } + + private: + ParallelIndexSet* indexSet_; + }; + + std::shared_ptr indexSet_; + std::shared_ptr remoteIndices_; + MPI_Comm communicator_; +}; +} // end namespace Opm +#endif +#endif diff --git a/tests/test_parallel_linearsolver.cpp b/tests/test_parallel_linearsolver.cpp new file mode 100644 index 000000000..aa670d276 --- /dev/null +++ b/tests/test_parallel_linearsolver.cpp @@ -0,0 +1,254 @@ +/* + Copyright 2014 Dr. Markus 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 . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif +#define NVERBOSE // to suppress our messages when throwing + +#define BOOST_TEST_MODULE OPM-ParallelIterativeSolverTest +#include + +// MPI header +#if HAVE_MPI +#include +#include +#include +#include +#include +#include +#include +#include +#else +#error "This file needs to compiled with MPI support!" +#endif +#include +#include + +#include +#include +#include + +struct MPIFixture { + MPIFixture() + { + int m_argc = boost::unit_test::framework::master_test_suite().argc; + char** m_argv = boost::unit_test::framework::master_test_suite().argv; + MPI_Init(&m_argc, &m_argv); + } + ~MPIFixture() + { + MPI_Finalize(); + } +}; + +BOOST_GLOBAL_FIXTURE(MPIFixture); + +struct MyMatrix +{ + MyMatrix(std::size_t rows, std::size_t nnz) + : data(nnz, 0.0), rowStart(rows+1, -1), + colIndex(nnz, -1) + {} + MyMatrix() + : data(), rowStart(), colIndex() + {} + + std::vector data; + std::vector rowStart; + std::vector colIndex; +}; + +typedef int LocalId; +typedef int GlobalId; +typedef Dune::OwnerOverlapCopyCommunication Communication; +typedef Dune::OwnerOverlapCopyAttributeSet GridAttributes; +typedef GridAttributes::AttributeSet GridFlag; +typedef Dune::ParallelLocalIndex LocalIndex; + +/// \brief Sets up a paralle Laplacian. +/// +/// The process stores the unknowns with indices in the range [start, end). +/// As we use an overlapping domain decomposition, the process owns the indices +/// in the range [istart, iend]. If we would only used the indices in this range then +/// they form a partitioning of the whole index set. +/// \tparam I The type of the parallel index set (for convenience) +/// \param indexset The parallel index set for marking owner and copy region. +/// \param N The global number of unknowns of the system. +/// \param start The first index stored on this process +/// \param end One past the last index stored on this process +/// \param istart The first index that the process owns. +/// \param iend One past the last index the process owns. +template +std::shared_ptr create1DLaplacian(I& indexset, int N, int start, int end, + int istart, int iend) +{ + indexset.beginResize(); + MyMatrix* mm=new MyMatrix(end-start, (end-start)*3); + int nnz=0; + mm->rowStart[0]=0; + assert(start==0||start=iend) + { + // We are in the overlap region of the grid + // therefore we setup the system such that + // right hand side will equal the left hand side + // of the linear system. + mm->colIndex[nnz]=localRow; + mm->data[nnz++]=1.0; + indexset.add(row, LocalIndex(localRow, GridAttributes::copy, true)); + mm->rowStart[localRow+1]=nnz; + continue; + } + + double dval=0; + if(row>0) + { + mm->colIndex[nnz]=localRow-1; + mm->data[nnz++]=-1; + dval+=1; + } + mm->colIndex[nnz]=localRow; + mm->data[nnz++]=2;//dval+(rowcolIndex[nnz]=localRow+1; + mm->data[nnz++]=-1; + dval+=1; + } + mm->rowStart[localRow+1]=nnz; + indexset.add(row, LocalIndex(localRow, GridAttributes::owner, true)); + } + mm->data.resize(nnz); + mm->colIndex.resize(nnz); + indexset.endResize(); + return std::shared_ptr(mm); +} + +template +void createRandomVectors(O& pinfo, int NN, std::vector& x, std::vector& b, + const MyMatrix& mat) +{ + x.resize(NN); + for(auto entry=x.begin(), end =x.end(); entry!=end; ++entry) + *entry=((double) (rand()%100))/10.0; + + pinfo.copyOwnerToAll(x,x); + + b.resize(NN); + + // Construct the right hand side as b=A*x + std::fill(b.begin(), b.end(), 0.0); + for(std::size_t row=0; row0) + start = istart - 1; + else + start = istart; + + if(iend x(end-start), b(end-start); + createRandomVectors(comm, end-start, x, b, *mat); + std::vector exact(x); + std::fill(x.begin(), x.end(), 0.0); + Opm::LinearSolverFactory ls(param); + boost::any anyComm(comm); + ls.solve(N, mat->data.size(), &(mat->rowStart[0]), + &(mat->colIndex[0]), &(mat->data[0]), &(b[0]), + &(x[0]), anyComm); +} + +#ifdef HAVE_DUNE_ISTL +BOOST_AUTO_TEST_CASE(CGAMGTest) +{ + Opm::parameter::ParameterGroup param; + param.insertParameter(std::string("linsolver"), std::string("istl")); + param.insertParameter(std::string("linsolver_type"), std::string("1")); + param.insertParameter(std::string("linsolver_max_iterations"), std::string("200")); + run_test(param); +} + +BOOST_AUTO_TEST_CASE(CGILUTest) +{ + Opm::parameter::ParameterGroup param; + param.insertParameter(std::string("linsolver"), std::string("istl")); + param.insertParameter(std::string("linsolver_type"), std::string("0")); + param.insertParameter(std::string("linsolver_max_iterations"), std::string("200")); + run_test(param); +} + +BOOST_AUTO_TEST_CASE(BiCGILUTest) +{ + Opm::parameter::ParameterGroup param; + param.insertParameter(std::string("linsolver"), std::string("istl")); + param.insertParameter(std::string("linsolver_type"), std::string("2")); + param.insertParameter(std::string("linsolver_max_iterations"), std::string("200")); + run_test(param); +} +#endif +