diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 332a36f58..0a48770a3 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -130,6 +130,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/NonlinearSolver_impl.hpp opm/autodiff/LinearisedBlackoilResidual.hpp opm/autodiff/ParallelDebugOutput.hpp + opm/autodiff/ParallelOverlappingILU0.hpp + opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp opm/autodiff/RateConverter.hpp opm/autodiff/RedistributeDataHandles.hpp opm/autodiff/SimulatorBase.hpp diff --git a/opm/autodiff/CPRPreconditioner.hpp b/opm/autodiff/CPRPreconditioner.hpp index afefca5b5..84afac7f5 100644 --- a/opm/autodiff/CPRPreconditioner.hpp +++ b/opm/autodiff/CPRPreconditioner.hpp @@ -48,6 +48,8 @@ #include #include +#include +#include namespace Opm { namespace @@ -95,7 +97,11 @@ struct CPRSelector > /// \brief The operator type; typedef Dune::OverlappingSchwarzOperator Operator; /// \brief The type of the preconditioner used for the elliptic part. - typedef Dune::BlockPreconditioner > + //typedef Dune::BlockPreconditioner > + //EllipticPreconditioner; + //typedef ParallelRestrictedOverlappingSchwarz > + //EllipticPreconditioner; + typedef ParallelOverlappingILU0 EllipticPreconditioner; /// \brief The type of the unique pointer to the preconditioner of the elliptic part. typedef std::unique_ptr > createEllipticPreconditionerPointer(const M& Ae, double relax, const Dune::OwnerOverlapCopyCommunication& comm) { - typedef Dune::BlockPreconditioner, - Dune::SeqILU0 > - ParallelPreconditioner; + typedef typename CPRSelector > + ::EllipticPreconditioner ParallelPreconditioner; - Dune::SeqILU0* ilu=new Dune::SeqILU0(Ae, relax); + //Dune::SeqILU0* ilu=new Dune::SeqILU0(Ae, relax); typedef typename CPRSelector > ::EllipticPreconditionerPointer EllipticPreconditionerPointer; - return EllipticPreconditionerPointer(new ParallelPreconditioner(*ilu, comm), - createParallelDeleter(*ilu, comm)); + return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax)); } #endif } // end namespace diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index b65194f71..a8088c154 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -27,6 +27,8 @@ #include #include #include +#include +#include #include #include #include @@ -141,43 +143,15 @@ namespace Opm #if HAVE_MPI typedef Dune::OwnerOverlapCopyCommunication Comm; - typedef Dune::BlockPreconditioner ParPreconditioner; - + typedef ParallelOverlappingILU0 ParPreconditioner; template - std::unique_ptr > + std::unique_ptr constructPrecond(Operator& opA, const Comm& comm) const { - typedef AdditionalObjectDeleter Deleter; - typedef std::unique_ptr Pointer; - int ilu_setup_successful = 1; - std::string message; + typedef std::unique_ptr Pointer; const double relax = 1.0; - SeqPreconditioner* seq_precond = nullptr; - try { - seq_precond = new SeqPreconditioner(opA.getmat(), relax); - } - catch ( Dune::MatrixBlockError error ) - { - message = error.what(); - std::cerr<<"Exception occured on process " << - comm.communicator().rank() << " during " << - "setup of ILU0 preconditioner with message: " << - message< #include diff --git a/opm/autodiff/ParallelOverlappingILU0.hpp b/opm/autodiff/ParallelOverlappingILU0.hpp new file mode 100644 index 000000000..e757f6932 --- /dev/null +++ b/opm/autodiff/ParallelOverlappingILU0.hpp @@ -0,0 +1,206 @@ +/* + Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015 Statoil AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#ifndef OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED +#define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED + +#include +#include + +namespace Opm +{ + +template +class ParallelOverlappingILU0; + +} // end namespace Opm + +namespace Dune +{ + +namespace Amg +{ + +/// \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. +/// \tparam Range The type of the Vector representing the range. +/// \tparam ParallelInfo The type of the parallel information object +/// used, e.g. Dune::OwnerOverlapCommunication +template +struct ConstructionTraits > +{ + typedef Dune::SeqILU0 T; + typedef DefaultParallelConstructionArgs Arguments; + typedef ConstructionTraits SeqConstructionTraits; + static inline Opm::ParallelOverlappingILU0* construct(Arguments& args) + { + return new Opm::ParallelOverlappingILU0(args.getMatrix(), + args.getComm(), + args.getArgs().relaxationFactor); + } + + static inline void deconstruct(Opm::ParallelOverlappingILU0* bp) + { + delete bp; + } + +}; + +} // end namespace Amg + +} // end namespace Dune + +namespace Opm +{ + +/// \brief A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as +/// +/// This preconditioner differs from a ParallelRestrictedOverlappingSchwarz with +/// Dune:SeqILU0 in the follwing way: +/// During apply we make sure that the current residual is consistent (i.e. +/// each process knows the same value for each index. The we solve +/// Ly= d for y and make y consistent again. Last we solve Ux = y and +/// make sure that x is consistent. +/// In contrast for ParallelRestrictedOverlappingSchwarz we solve (LU)x = d for x +/// without forcing consistency between the two steps. +/// \tparam Matrix The type of the Matrix. +/// \tparam Domain The type of the Vector representing the domain. +/// \tparam Range The type of the Vector representing the range. +/// \tparam ParallelInfo The type of the parallel information object +/// used, e.g. Dune::OwnerOverlapCommunication +template +class ParallelOverlappingILU0 + : public Dune::Preconditioner { +public: + //! \brief The matrix type the preconditioner is for. + typedef typename Dune::remove_const::type matrix_type; + //! \brief The domain type of the preconditioner. + typedef Domain domain_type; + //! \brief The range type of the preconditioner. + typedef Range range_type; + //! \brief The field type of the preconditioner. + typedef typename Domain::field_type field_type; + + // define the category + enum { + //! \brief The category the preconditioner is part of. + category=Dune::SolverCategory::overlapping + }; + + /*! \brief Constructor. + + Constructor gets all parameters to operate the prec. + \param A The matrix to operate on. + \param w The relaxation factor. + */ + ParallelOverlappingILU0 (const Matrix& A, const ParallelInfo& comm, + field_type w) + : ilu_(A), comm_(comm), w_(w) + { + int ilu_setup_successful = 1; + std::string message; + try + { + bilu0_decomposition(ilu_); + } + catch ( Dune::MatrixBlockError error ) + { + message = error.what(); + std::cerr<<"Exception occured on process " << + comm.communicator().rank() << " during " << + "setup of ILU0 preconditioner with message: " << + message<(d); + comm_.copyOwnerToAll(md,md); + auto endrow=ilu_.end(); + for ( auto row = ilu_.begin(); row != endrow; ++row ) + { + auto rhs(d[row.index()]); + for ( auto col = row->begin(); col.index() < row.index(); ++col ) + { + col->mmv(v[col.index()],rhs); + } + v[row.index()] = rhs; + } + comm_.copyOwnerToAll(v, v); + auto rendrow = ilu_.beforeBegin(); + for( auto row = ilu_.beforeEnd(); row != rendrow; --row) + { + auto rhs(v[row.index()]); + auto col = row->beforeEnd(); + for( ; col.index() > row.index(); --col) + { + col->mmv(v[col.index()], rhs); + } + v[row.index()] = 0; + col->umv(rhs, v[row.index()]); + } + comm_.copyOwnerToAll(v, v); + v *= w_; + } + + /*! + \brief Clean up. + + \copydoc Preconditioner::post(X&) + */ + virtual void post (Range& x) + { + DUNE_UNUSED_PARAMETER(x); + } + +private: + //! \brief The ILU0 decomposition of the matrix. + Matrix ilu_; + const ParallelInfo& comm_; + //! \brief The relaxation factor to use. + field_type w_; + +}; + +} // end namespace Opm +#endif diff --git a/opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp b/opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp new file mode 100644 index 000000000..b5ec98023 --- /dev/null +++ b/opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp @@ -0,0 +1,209 @@ +/* + Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015 Statoil AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#ifndef OPM_PARALLELRESTRICTEDADDITIVESCHWARZ_HEADER_INCLUDED +#define OPM_PARALLELRESTRICTEDADDITIVESCHWARZ_HEADER_INCLUDED + +#include +#include + +namespace Opm +{ + +template +class ParallelRestrictedOverlappingSchwarz; + +} // end namespace Opm + +namespace Dune +{ + +namespace Amg +{ + +/// \brief Tells AMG how to construct the Opm::ParallelOverlappingILU0 smoother +/// \tparam Domain The type of the Vector representing the domain. +/// \tparam Range The type of the Vector representing the range. +/// \tparam ParallelInfo The type of the parallel information object +/// used, e.g. Dune::OwnerOverlapCommunication +/// \tparam SeqPreconditioner The underlying sequential preconditioner to use. +template +struct ConstructionTraits > +{ + typedef DefaultParallelConstructionArgs Arguments; + typedef ConstructionTraits SeqConstructionTraits; + + /// \brief Construct a parallel restricted overlapping schwarz preconditioner. + static inline Opm::ParallelRestrictedOverlappingSchwarz* + construct(Arguments& args) + { + return new Opm::ParallelRestrictedOverlappingSchwarz + (*SeqConstructionTraits + ::construct(args), + args.getComm()); + } + + /// \brief Deconstruct and free a parallel restricted overlapping schwarz preconditioner. + static inline void deconstruct(Opm::ParallelRestrictedOverlappingSchwarz + * bp) + { + SeqConstructionTraits + ::deconstruct(static_cast(&bp->preconditioner)); + delete bp; + } + +}; + +/// \brief Tells AMG how to use Opm::ParallelOverlappingILU0 smoother +/// \tparam Domain The type of the Vector representing the domain. +/// \tparam Range The type of the Vector representing the range. +/// \tparam ParallelInfo The type of the parallel information object +/// used, e.g. Dune::OwnerOverlapCommunication +/// \tparam SeqPreconditioner The underlying sequential preconditioner to use. +template +struct SmootherTraits > +{ + typedef DefaultSmootherArgs Arguments; + +}; + +} // end namespace Amg + +} // end namespace Dune + +namespace Opm{ + +/// \brief Block parallel preconditioner. +/// +/// This is essentially a wrapper that takes a sequential +/// preconditioner. In each step the sequential preconditioner +/// is applied to the whole subdomain and then all owner data +/// points are updated on all other processes from the processor +/// that knows the complete matrix row for this data point (in dune-istl +/// speak that is the one that owns the data). +/// +/// Note that this is different from the usual approach in dune-istl where +/// the application of the sequential preconditioner only takes place on +/// the (owner) partition of the process disregarding any overlap/ghost region. +/// +/// For more information see https://www.cs.colorado.edu/~cai/papers/rash.pdf +/// +/// \tparam Domain The type of the Vector representing the domain. +/// \tparam Range The type of the Vector representing the range. +/// \tparam ParallelInfo The type of the parallel information object +/// used, e.g. Dune::OwnerOverlapCommunication +/// \tparam SeqPreconditioner The underlying sequential preconditioner to use. +template > +class ParallelRestrictedOverlappingSchwarz + : public Dune::Preconditioner { + friend class Dune::Amg + ::ConstructionTraits >; +public: + //! \brief The domain type of the preconditioner. + typedef Domain domain_type; + //! \brief The range type of the preconditioner. + typedef Range range_type; + //! \brief The field type of the preconditioner. + typedef typename Domain::field_type field_type; + //! \brief The type of the communication object. + typedef ParallelInfo communication_type; + + // define the category + enum { + //! \brief The category the precondtioner is part of. + category=Dune::SolverCategory::overlapping + }; + + /*! \brief Constructor. + + constructor gets all parameters to operate the prec. + \param p The sequential preconditioner. + \param c The communication object for syncing overlap and copy + data points. (E.~g. OwnerOverlapCommunication ) + */ + ParallelRestrictedOverlappingSchwarz (SeqPreconditioner& p, const communication_type& c) + : preconditioner_(p), communication_(c) + { } + + /*! + \brief Prepare the preconditioner. + + \copydoc Preconditioner::pre(X&,Y&) + */ + virtual void pre (Domain& x, Range& b) + { + communication_.copyOwnerToAll(x,x); // make dirichlet values consistent + preconditioner_.pre(x,b); + } + + /*! + \brief Apply the preconditioner + + \copydoc Preconditioner::apply(X&,const Y&) + */ + virtual void apply (Domain& v, const Range& d) + { + apply(v, d); + } + + template + void apply (Domain& v, const Range& d) + { + // hack us a mutable d to prevent copying. + Range& md = const_cast(d); + communication_.copyOwnerToAll(md,md); + preconditioner_.template apply(v,d); + communication_.copyOwnerToAll(v,v); + // Make sure that d is the same as at the beginning of apply. + communication_.project(md); + } + + /*! + \brief Clean up. + + \copydoc Preconditioner::post(X&) + */ + virtual void post (Range& x) + { + preconditioner_.post(x); + } + +private: + //! \brief a sequential preconditioner + SeqPreconditioner& preconditioner_; + + //! \brief the communication object + const communication_type& communication_; +}; + + +} // end namespace OPM +#endif diff --git a/opm/autodiff/RedistributeDataHandles.hpp b/opm/autodiff/RedistributeDataHandles.hpp index 5ec65988f..397f53942 100644 --- a/opm/autodiff/RedistributeDataHandles.hpp +++ b/opm/autodiff/RedistributeDataHandles.hpp @@ -344,7 +344,7 @@ void distributeGridAndData( Dune::CpGrid& grid, global_grid.switchToGlobalView(); // distribute the grid and switch to the distributed view - grid.loadBalance(eclipseState); + grid.loadBalance(eclipseState, geology.transmissibility().data()); grid.switchToDistributedView(); std::vector compressedToCartesianIdx; Opm::createGlobalCellArray(grid, compressedToCartesianIdx);