diff --git a/CMakeLists.txt b/CMakeLists.txt index 18c2c96f3..3e531de34 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -157,6 +157,16 @@ opm_add_test(flow flow/flow_ebos_oilwater_polymer.cpp flow/flow_ebos_oilwater_polymer_injectivity.cpp) +opm_add_test(flow_blackoil_dunecpr + ONLY_COMPILE + DEFAULT_ENABLE_IF ${FLOW_DEFAULT_ENABLE_IF} + SOURCES flow/flow_blackoil_dunecpr.cpp + EXE_NAME flow_blackoil_dunecpr + DEPENDS "opmsimulators" + LIBRARIES "opmsimulators") + + + if (BUILD_FLOW) install(TARGETS flow DESTINATION bin) opm_add_bash_completion(flow) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 8a34d6a26..638bf4fc8 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -118,6 +118,9 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/AquiferCarterTracy.hpp opm/autodiff/AquiferFetkovich.hpp opm/autodiff/BlackoilAmg.hpp + opm/autodiff/BlackoilAmgCpr.hpp + opm/autodiff/amgcpr.hh + opm/autodiff/twolevelmethodcpr.hh opm/autodiff/BlackoilDetails.hpp opm/autodiff/BlackoilModelParametersEbos.hpp opm/autodiff/BlackoilAquiferModel.hpp @@ -129,6 +132,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/FlowMainEbos.hpp opm/autodiff/GraphColoring.hpp opm/autodiff/ISTLSolverEbos.hpp + opm/autodiff/ISTLSolverEbosCpr.hpp opm/autodiff/IterationReport.hpp opm/autodiff/MatrixBlock.hpp opm/autodiff/moduleVersion.hpp diff --git a/flow/flow_blackoil_dunecpr.cpp b/flow/flow_blackoil_dunecpr.cpp new file mode 100644 index 000000000..efb95b9dc --- /dev/null +++ b/flow/flow_blackoil_dunecpr.cpp @@ -0,0 +1,96 @@ +/* + Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015, 2017 IRIS 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 . +*/ +#include "config.h" +#include "flow/flow_tag.hpp" +//#include +#include +//#include + +BEGIN_PROPERTIES +NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem)); +NEW_PROP_TAG(FluidState); +//SET_TYPE_PROP(EclBaseProblem, Problem, Ewoms::EclProblem); +SET_PROP(EclFlowProblemSimple, FluidState) + { + private: + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; + enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + static const bool compositionSwitchEnabled = Indices::gasEnabled; + + public: +//typedef Opm::BlackOilFluidSystemSimple type; + typedef Opm::BlackOilFluidState type; +}; + +SET_BOOL_PROP(EclFlowProblemSimple,MatrixAddWellContributions,true); +SET_INT_PROP(EclFlowProblemSimple,LinearSolverVerbosity,1); +SET_SCALAR_PROP(EclFlowProblemSimple, LinearSolverReduction, 1e-4); +SET_INT_PROP(EclFlowProblemSimple, LinearSolverMaxIter, 20); +SET_BOOL_PROP(EclFlowProblemSimple, UseAmg, true);//probably not used +SET_BOOL_PROP(EclFlowProblemSimple, UseCpr, true); +SET_INT_PROP(EclFlowProblemSimple, CprMaxEllIter, 1); +SET_INT_PROP(EclFlowProblemSimple, CprEllSolvetype, 3); +SET_INT_PROP(EclFlowProblemSimple, CprReuseSetup, 3); +SET_INT_PROP(EclFlowProblemSimple, CprSolverVerbose, 3); +SET_STRING_PROP(EclFlowProblemSimple, SystemStrategy, "quasiimpes"); +END_PROPERTIES + +namespace Ewoms { + namespace Properties { + + SET_PROP(EclFlowProblemSimple, FluidSystem) + { + private: + //typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + + public: + typedef Opm::BlackOilFluidSystem type; + }; + //NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem)); + SET_TYPE_PROP(EclFlowProblemSimple, IntensiveQuantities, Ewoms::BlackOilIntensiveQuantities); + //SET_TYPE_PROP(EclFlowProblemSimple, LinearSolverBackend, Opm::ISTLSolverEbos); + //SET_TAG_PROP(EclFlowProblemSimple, LinearSolverSplice, ParallelBiCGStabLinearSolver); + //SET_TYPE_PROP(EclFlowProblemSimple, LinearSolverBackend, Ewoms::Linear::ParallelBiCGStabSolverBackend);//not work + //SET_TYPE_PROP(EclFlowProblemSimple, LinearSolverBackend, Ewoms::Linear::SuperLUBackend)//not work + //SET_TAG_PROP(EclFlowProblem, FluidState, Opm::BlackOilFluidState); + SET_TYPE_PROP(EclFlowProblemSimple, LinearSolverBackend, Opm::ISTLSolverEbosCpr); + SET_BOOL_PROP(EclFlowProblemSimple, EnableStorageCache, true); + SET_BOOL_PROP(EclFlowProblemSimple, EnableIntensiveQuantityCache, true); + + //SET_INT_PROP(EclFlowProblemSimple, NumWellAdjoint, 1); + //SET_BOOL_PROP(EclFlowProblem, EnableStorageCache, true); + //SET_BOOL_PROP(EclFlowProblem, EnableIntensiveQuantityCache, true); + } +} + +int main(int argc, char** argv) +{ + typedef TTAG(EclFlowProblemSimple) TypeTag; + return mainFlow(argc, argv); +} diff --git a/flow/flow_tag.hpp b/flow/flow_tag.hpp new file mode 100644 index 000000000..d879fe82f --- /dev/null +++ b/flow/flow_tag.hpp @@ -0,0 +1,212 @@ +/* + Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015, 2017 IRIS 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 FLOW_TAG_HPP +#define FLOW_TAG_HPP +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +//#include +//#include +#include +#include +//#include + +#if HAVE_DUNE_FEM +#include +#else +#include +#endif + + +BEGIN_PROPERTIES + +// this is a dummy type tag that is used to setup the parameters before the actual +// simulator. +NEW_TYPE_TAG(FlowEarlyBird, INHERITS_FROM(EclFlowProblem)); + +END_PROPERTIES + + + + +namespace Opm { + template + void flowEbosSetDeck(Deck &deck, EclipseState& eclState) + { + typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard; + Vanguard::setExternalDeck(&deck, &eclState); + } + +// ----------------- Main program ----------------- + template + int flowEbosMain(int argc, char** argv) + { + // we always want to use the default locale, and thus spare us the trouble + // with incorrect locale settings. + Opm::resetLocale(); + +#if HAVE_DUNE_FEM + Dune::Fem::MPIManager::initialize(argc, argv); +#else + Dune::MPIHelper::instance(argc, argv); +#endif + Opm::FlowMainEbos mainfunc; + return mainfunc.execute(argc, argv); + } + +} + + + + + +namespace detail +{ + boost::filesystem::path simulationCaseName( const std::string& casename ) { + namespace fs = boost::filesystem; + + const auto exists = []( const fs::path& f ) -> bool { + if( !fs::exists( f ) ) return false; + + if( fs::is_regular_file( f ) ) return true; + + return fs::is_symlink( f ) + && fs::is_regular_file( fs::read_symlink( f ) ); + }; + + auto simcase = fs::path( casename ); + + if( exists( simcase ) ) { + return simcase; + } + + for( const auto& ext : { std::string("data"), std::string("DATA") } ) { + if( exists( simcase.replace_extension( ext ) ) ) { + return simcase; + } + } + + throw std::invalid_argument( "Cannot find input case " + casename ); + } +} + + +// ----------------- Main program ----------------- +template +int mainFlow(int argc, char** argv) +{ + // MPI setup. +#if HAVE_DUNE_FEM + Dune::Fem::MPIManager::initialize(argc, argv); + int mpiRank = Dune::Fem::MPIManager::rank(); +#else + // the design of the plain dune MPIHelper class is quite flawed: there is no way to + // get the instance without having the argc and argv parameters available and it is + // not possible to determine the MPI rank and size without an instance. (IOW: the + // rank() and size() methods are supposed to be static.) + const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); + int mpiRank = mpiHelper.rank(); +#endif + + // we always want to use the default locale, and thus spare us the trouble + // with incorrect locale settings. + Opm::resetLocale(); + + // this is a work-around for a catch 22: we do not know what code path to use without + // parsing the deck, but we don't know the deck without having access to the + // parameters and this requires to know the type tag to be used. To solve this, we + // use a type tag just for parsing the parameters before we instantiate the actual + // simulator object. (Which parses the parameters again, but since this is done in an + // identical manner it does not matter.) + typedef TTAG(FlowEarlyBird) PreTypeTag; + typedef GET_PROP_TYPE(PreTypeTag, Problem) PreProblem; + + PreProblem::setBriefDescription("Simple Flow, an advanced reservoir simulator for ECL-decks provided by the Open Porous Media project."); + + int status = Opm::FlowMainEbos::setupParameters_(argc, argv); + if (status != 0) + // if setupParameters_ returns a value smaller than 0, there was no error, but + // the program should abort. This is the case e.g. for the --help and the + // --print-properties parameters. + return (status >= 0)?status:0; + + bool outputCout = false; + if (mpiRank == 0) + outputCout = EWOMS_GET_PARAM(PreTypeTag, bool, EnableTerminalOutput); + + std::string deckFilename = EWOMS_GET_PARAM(PreTypeTag, std::string, EclDeckFileName); + typedef typename GET_PROP_TYPE(PreTypeTag, Vanguard) PreVanguard; + try { + deckFilename = PreVanguard::canonicalDeckPath(deckFilename).string(); + } + catch (const std::exception& e) { + Ewoms::Parameters::printUsage(PreProblem::helpPreamble(argc, const_cast(argv)), + e.what()); + return 1; + } + + // Create Deck and EclipseState. + try { + Opm::Parser parser; + typedef std::pair ParseModePair; + typedef std::vector ParseModePairs; + ParseModePairs tmp; + tmp.push_back(ParseModePair(Opm::ParseContext::PARSE_RANDOM_SLASH, Opm::InputError::IGNORE)); + tmp.push_back(ParseModePair(Opm::ParseContext::PARSE_MISSING_DIMS_KEYWORD, Opm::InputError::WARN)); + tmp.push_back(ParseModePair(Opm::ParseContext::SUMMARY_UNKNOWN_WELL, Opm::InputError::WARN)); + tmp.push_back(ParseModePair(Opm::ParseContext::SUMMARY_UNKNOWN_GROUP, Opm::InputError::WARN)); + Opm::ParseContext parseContext(tmp); + Opm::ErrorGuard errorGuard; + + std::shared_ptr deck = std::make_shared< Opm::Deck >( parser.parseFile(deckFilename , parseContext, errorGuard) ); + if ( outputCout ) { + Opm::checkDeck(*deck, parser, parseContext, errorGuard); + Opm::MissingFeatures::checkKeywords(*deck); + } + //Opm::Runspec runspec( *deck ); + //const auto& phases = runspec.phases(); + + std::shared_ptr eclipseState = std::make_shared< Opm::EclipseState > ( *deck, parseContext, errorGuard ); + Opm::flowEbosSetDeck(*deck, *eclipseState); + return Opm::flowEbosMain(argc, argv); + } + catch (const std::invalid_argument& e) + { + if (outputCout) { + std::cerr << "Failed to create valid EclipseState object." << std::endl; + std::cerr << "Exception caught: " << e.what() << std::endl; + } + throw; + } + + return EXIT_SUCCESS; +} +#endif diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index 5f60c8528..2560f69ad 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -23,6 +23,8 @@ #include #include #include +#include +#include #include #include #include @@ -77,18 +79,6 @@ namespace Opm namespace Detail { -/** - * \brief Creates a MatrixAdapter as an operator - * - * The first argument is used to specify the return type using function overloading. - * \param matrix The matrix to wrap. - */ -template -Dune::MatrixAdapter createOperator(const Dune::MatrixAdapter&, const M& matrix, const T&) -{ - return Dune::MatrixAdapter(matrix); -} - /** * \brief Creates a MatrixAdapter as an operator, storing it in a unique_ptr. * @@ -96,22 +86,24 @@ Dune::MatrixAdapter createOperator(const Dune::MatrixAdapter&, con * \param matrix The matrix to wrap. */ template -std::unique_ptr< Dune::MatrixAdapter > createOperatorPtr(const Dune::MatrixAdapter&, const M& matrix, const T&) +std::unique_ptr< Dune::MatrixAdapter > +createOperator(const Dune::MatrixAdapter&, const M& matrix, const T&) { return std::make_unique< Dune::MatrixAdapter >(matrix); } + /** - * \brief Creates an OverlappingSchwarzOperator as an operator. + * \brief Creates an OverlappingSchwarzOperator as an operator, storing it in a unique_ptr. * * The first argument is used to specify the return type using function overloading. * \param matrix The matrix to wrap. * \param comm The object encapsulating the parallelization information. */ template -Dune::OverlappingSchwarzOperator createOperator(const Dune::OverlappingSchwarzOperator&, - const M& matrix, const T& comm) +std::unique_ptr< Dune::OverlappingSchwarzOperator > +createOperator(const Dune::OverlappingSchwarzOperator&, const M& matrix, const T& comm) { - return Dune::OverlappingSchwarzOperator(matrix, comm); + return std::make_unique< Dune::OverlappingSchwarzOperator >(matrix, comm); } //! \brief Applies diagonal scaling to the discretization Matrix (Scheichl, 2003) @@ -122,10 +114,9 @@ Dune::OverlappingSchwarzOperator createOperator(const Dune::Overlapping //! \param comm The communication objecte describing the data distribution. //! \param pressureEqnIndex The index of the pressure in the matrix block //! \retun A pair of the scaled matrix and the associated operator- -template -std::tuple, Operator> -scaleMatrixDRS(const Operator& op, const Communication& comm, - std::size_t pressureEqnIndex, const Vector& weights, const Opm::CPRParameter& param) +template +std::unique_ptr +scaleMatrixDRS(const Operator& op, std::size_t pressureEqnIndex, const Vector& weights, const Opm::CPRParameter& param) { using Matrix = typename Operator::matrix_type; using Block = typename Matrix::block_type; @@ -144,7 +135,7 @@ scaleMatrixDRS(const Operator& op, const Communication& comm, } } } - return std::make_tuple(std::move(matrix), createOperator(op, *matrix, comm)); + return matrix; } //! \brief Applies diagonal scaling to the discretization Matrix (Scheichl, 2003) @@ -356,7 +347,7 @@ public: /** @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; + typedef Dune::Amg::AMGCPR AMGType; /** * @brief Constructs the coarse solver policy. * @param args The arguments used for constructing the smoother. @@ -384,7 +375,7 @@ private: const Criterion& crit, const typename AMGType::SmootherArgs& args, const Communication& comm) - : param_(param), amg_(), smoother_(), op_(op), comm_(comm) + : param_(param), amg_(), smoother_(), op_(op), crit_(crit), comm_(comm) { if ( param_->cpr_use_amg_ ) { @@ -400,6 +391,11 @@ private: } } + void updateAmgPreconditioner(typename AMGType::Operator& op) + { + amg_->updateSolver(crit_, op, comm_); + } + #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) Dune::SolverCategory::Category category() const override { @@ -408,7 +404,7 @@ private: } #endif - void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) + void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override { DUNE_UNUSED_PARAMETER(reduction); DUNE_UNUSED_PARAMETER(res); @@ -518,7 +514,7 @@ private: #endif } - void apply(X& x, X& b, Dune::InverseOperatorResult& res) + void apply(X& x, X& b, Dune::InverseOperatorResult& res) override { return apply(x,b,1e-8,res); } @@ -535,6 +531,7 @@ private: std::unique_ptr amg_; std::unique_ptr smoother_; const typename AMGType::Operator& op_; + Criterion crit_; const Communication& comm_; }; @@ -560,8 +557,13 @@ public: smootherArgs_, transfer.getCoarseLevelCommunication()); - return inv; //std::shared_ptr >(inv); + return inv; + } + template + void setCoarseOperator(LTP& transferPolicy) + { + coarseOperator_= transferPolicy.getCoarseLevelOperator(); } private: @@ -706,7 +708,7 @@ void buildCoarseSparseMatrix(M& coarseMatrix, G& fineGraph, const V& visitedMap, */ template class OneComponentAggregationLevelTransferPolicy - : public Dune::Amg::LevelTransferPolicy::value> + : public Dune::Amg::LevelTransferPolicyCpr::value> { typedef Dune::Amg::AggregatesMap AggregatesMap; public: @@ -789,7 +791,7 @@ public: { coarseLevelCommunication_->freeGlobalLookup(); } - calculateCoarseEntries(fineOperator.getmat()); + calculateCoarseEntriesWithAggregatesMap(fineOperator); } else { @@ -829,10 +831,10 @@ public: this->operator_.reset(Dune::Amg::ConstructionTraits::construct(oargs)); } - template - void calculateCoarseEntries(const M& fineMatrix) + void calculateCoarseEntriesWithAggregatesMap(const Operator& fineOperator) { - *coarseLevelMatrix_ = 0; + const auto& fineMatrix = fineOperator.getmat(); + *coarseLevelMatrix_ = 0; for(auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); row != rowEnd; ++row) { @@ -852,6 +854,23 @@ public: } } + virtual void calculateCoarseEntries(const Operator& fineOperator) + { + const auto& fineMatrix = fineOperator.getmat(); + *coarseLevelMatrix_ = 0; + for(auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); + row != rowEnd; ++row) + { + const auto& i = row.index(); + for(auto entry = row->begin(), entryEnd = row->end(); + entry != entryEnd; ++entry) + { + const auto& j = entry.index(); + (*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][VARIABLE_INDEX]; + } + } + } + void moveToCoarseLevel(const typename FatherType::FineRangeType& fine) { // Set coarse vector to zero @@ -980,9 +999,9 @@ protected: CoarseCriterion, LevelTransferPolicy>; using TwoLevelMethod = - Dune::Amg::TwoLevelMethod; + Dune::Amg::TwoLevelMethodCpr; public: #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) Dune::SolverCategory::Category category() const override @@ -1011,29 +1030,29 @@ public: const SmootherArgs& smargs, const Communication& comm) : param_(param), weights_(weights), - scaledMatrixOperator_(Detail::scaleMatrixDRS(fineOperator, comm, - COMPONENT_INDEX, weights, param)), - smoother_(Detail::constructSmoother(std::get<1>(scaledMatrixOperator_), + scaledMatrix_(Detail::scaleMatrixDRS(fineOperator, COMPONENT_INDEX, weights, param)), + scaledMatrixOperator_(Detail::createOperator(fineOperator, *scaledMatrix_, comm)), + smoother_(Detail::constructSmoother(*scaledMatrixOperator_, smargs, comm)), levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_), coarseSolverPolicy_(¶m, smargs, criterion), - twoLevelMethod_(std::get<1>(scaledMatrixOperator_), smoother_, + twoLevelMethod_(*scaledMatrixOperator_, smoother_, levelTransferPolicy_, coarseSolverPolicy_, 0, 1) {} void pre(typename TwoLevelMethod::FineDomainType& x, - typename TwoLevelMethod::FineRangeType& b) + typename TwoLevelMethod::FineRangeType& b) override { twoLevelMethod_.pre(x,b); } - void post(typename TwoLevelMethod::FineDomainType& x) + void post(typename TwoLevelMethod::FineDomainType& x) override { twoLevelMethod_.post(x); } void apply(typename TwoLevelMethod::FineDomainType& v, - const typename TwoLevelMethod::FineRangeType& d) + const typename TwoLevelMethod::FineRangeType& d) override { auto scaledD = d; Detail::scaleVectorDRS(scaledD, COMPONENT_INDEX, param_, weights_); @@ -1042,7 +1061,8 @@ public: private: const CPRParameter& param_; const typename TwoLevelMethod::FineDomainType& weights_; - std::tuple, Operator> scaledMatrixOperator_; + std::unique_ptr scaledMatrix_; + std::unique_ptr scaledMatrixOperator_; std::shared_ptr smoother_; LevelTransferPolicy levelTransferPolicy_; CoarseSolverPolicy coarseSolverPolicy_; diff --git a/opm/autodiff/BlackoilAmgCpr.hpp b/opm/autodiff/BlackoilAmgCpr.hpp new file mode 100644 index 000000000..ae6ece976 --- /dev/null +++ b/opm/autodiff/BlackoilAmgCpr.hpp @@ -0,0 +1,181 @@ +/* + 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_AMGCPR_HEADER_INCLUDED +#define OPM_AMGCPR_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + /** + * \brief An algebraic twolevel or multigrid approach for solving blackoil (supports CPR with and without AMG) + * + * This preconditioner first decouples the component used for coarsening using a simple scaling + * approach (e.g. Scheichl, Masson 2013,\see scaleMatrixDRS). Then it constructs the + * coarse level system. The coupling is defined by the weights corresponding to the element located at + * (COMPONENT_INDEX, VARIABLE_INDEX) in the block matrix. Then the coarse level system is constructed + * either by extracting these elements, or by applying aggregation to them directly. This coarse level + * can be solved either by AMG or by ILU. The preconditioner is configured using CPRParameter. + * \tparam O The type of the operator (encapsulating a BCRSMatrix). + * \tparam S The type of the smoother. + * \tparam C The type of coarsening criterion to use. + * \tparam P The type of the class describing the parallelization. + * \tparam COMPONENT_INDEX The index of the component to use for coarsening (usually water). + * \tparam VARIABLE_INDEX The index of the variable to use for coarsening (usually pressure). + */ + template + class BlackoilAmgCpr + : public Dune::Preconditioner + { + public: + /** \brief The type of the operator (encapsulating a BCRSMatrix). */ + using Operator = O; + /** \brief The type of coarsening criterion to use. */ + using Criterion = C; + /** \brief The type of the class describing the parallelization. */ + using Communication = P; + /** \brief The type of the smoother. */ + using Smoother = S; + /** \brief The type of the smoother arguments for construction. */ + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + + protected: + using Matrix = typename Operator::matrix_type; + using CoarseOperator = typename Detail::ScalarType::value; + using CoarseSmoother = typename Detail::ScalarType::value; + using FineCriterion = + typename Detail::OneComponentCriterionType::value; + using CoarseCriterion = typename Detail::ScalarType::value; + using LevelTransferPolicy = + OneComponentAggregationLevelTransferPolicy; + using CoarseSolverPolicy = + Detail::OneStepAMGCoarseSolverPolicy; + using TwoLevelMethod = + Dune::Amg::TwoLevelMethodCpr; + public: +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::SolverCategory::Category category() const override + { + return std::is_same::value ? + Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping; + } +#else + // define the category + enum { + //! \brief The category the precondtioner is part of. + category = Operator::category + }; +#endif + + /** + * \brief Constructor. + * \param param The parameters used for configuring the solver. + * \param fineOperator The operator of the fine level. + * \param criterion The criterion describing the coarsening approach. + * \param smargs The arguments for constructing the smoother. + * \param comm The information about the parallelization. + */ + BlackoilAmgCpr(const CPRParameter& param, + const typename TwoLevelMethod::FineDomainType& weights, + const Operator& fineOperator, const Criterion& criterion, + const SmootherArgs& smargs, const Communication& comm) + : param_(param), + weights_(weights), + scaledMatrix_(Detail::scaleMatrixDRS(fineOperator, COMPONENT_INDEX, weights_, param)), + scaledMatrixOperator_(Detail::createOperator(fineOperator, *scaledMatrix_, comm)), + smoother_(Detail::constructSmoother(*scaledMatrixOperator_, + smargs, comm)), + levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_), + coarseSolverPolicy_(¶m, smargs, criterion), + twoLevelMethod_(*scaledMatrixOperator_, + smoother_, + levelTransferPolicy_, + coarseSolverPolicy_, 0, 1) + { + } + + void updatePreconditioner(const Operator& fineOperator, + const SmootherArgs& smargs, + const Communication& comm) + { + *scaledMatrix_ = *Detail::scaleMatrixDRS(fineOperator, COMPONENT_INDEX, weights_, param_); + smoother_.reset(Detail::constructSmoother(*scaledMatrixOperator_, smargs, comm)); + twoLevelMethod_.updatePreconditioner(*scaledMatrixOperator_, + smoother_, + coarseSolverPolicy_); + } + + void pre(typename TwoLevelMethod::FineDomainType& x, + typename TwoLevelMethod::FineRangeType& b) override + { + twoLevelMethod_.pre(x,b); + } + + void post(typename TwoLevelMethod::FineDomainType& x) override + { + twoLevelMethod_.post(x); + } + + void apply(typename TwoLevelMethod::FineDomainType& v, + const typename TwoLevelMethod::FineRangeType& d) override + { + auto scaledD = d; + Detail::scaleVectorDRS(scaledD, COMPONENT_INDEX, param_, weights_); + twoLevelMethod_.apply(v, scaledD); + } + + private: + const CPRParameter& param_; + const typename TwoLevelMethod::FineDomainType& weights_; + std::unique_ptr scaledMatrix_; + std::unique_ptr scaledMatrixOperator_; + std::shared_ptr smoother_; + LevelTransferPolicy levelTransferPolicy_; + CoarseSolverPolicy coarseSolverPolicy_; + TwoLevelMethod twoLevelMethod_; + }; + +} // end namespace Opm +#endif diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 92f562672..64aaa51b8 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -298,12 +298,16 @@ namespace Opm { wellModel().linearize(ebosSimulator().model().linearizer().jacobian(), ebosSimulator().model().linearizer().residual()); + // Solve the linear system. + linear_solve_setup_time_ = 0.0; try { solveJacobianSystem(x); + report.linear_solve_setup_time += linear_solve_setup_time_; report.linear_solve_time += perfTimer.stop(); report.total_linear_iterations += linearIterationsLastSolve(); } catch (...) { + report.linear_solve_setup_time += linear_solve_setup_time_; report.linear_solve_time += perfTimer.stop(); report.total_linear_iterations += linearIterationsLastSolve(); @@ -473,7 +477,10 @@ namespace Opm { x = 0.0; auto& ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver(); + Dune::Timer perfTimer; + perfTimer.start(); ebosSolver.prepare(ebosJac, ebosResid); + linear_solve_setup_time_ = perfTimer.stop(); ebosSolver.setResidual(ebosResid); // actually, the error needs to be calculated after setResidual in order to // account for parallelization properly. since the residual of ECFV @@ -906,7 +913,7 @@ namespace Opm { double dsMax() const { return param_.ds_max_; } double drMaxRel() const { return param_.dr_max_rel_; } double maxResidualAllowed() const { return param_.max_residual_allowed_; } - + double linear_solve_setup_time_; public: std::vector wasSwitched_; }; diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index 1b56b09f2..68b80ee08 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -76,6 +76,7 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) // Implementation for ISTL-matrix based operator //===================================================================== + /*! \brief Adapter to turn a matrix into a linear operator. @@ -129,8 +130,16 @@ public: } #endif } + WellModelMatrixAdapter (const M& A, + const M& A_for_precond, + const WellModel& wellMod, + std::shared_ptr comm ) + : A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), comm_(comm) + { + } - virtual void apply( const X& x, Y& y ) const + + virtual void apply( const X& x, Y& y ) const override { A_.mv( x, y ); @@ -144,7 +153,7 @@ public: } // y += \alpha * A * x - virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const + virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override { A_.usmv(alpha,x,y); @@ -157,18 +166,18 @@ public: #endif } - virtual const matrix_type& getmat() const { return A_for_precond_; } + virtual const matrix_type& getmat() const override { return A_for_precond_; } - communication_type* comm() + std::shared_ptr comm() { - return comm_.operator->(); + return comm_; } protected: const matrix_type& A_ ; const matrix_type& A_for_precond_ ; - const WellModel& wellMod_; - std::unique_ptr< communication_type > comm_; + const WellModel& wellMod_; + std::shared_ptr< communication_type > comm_; }; /// This class solves the fully implicit black-oil system by @@ -178,6 +187,7 @@ protected: template class ISTLSolverEbos { + protected: typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; diff --git a/opm/autodiff/ISTLSolverEbosCpr.hpp b/opm/autodiff/ISTLSolverEbosCpr.hpp new file mode 100644 index 000000000..3d07f3300 --- /dev/null +++ b/opm/autodiff/ISTLSolverEbosCpr.hpp @@ -0,0 +1,319 @@ +/* + Copyright 2016 IRIS 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_ISTLSOLVERCPR_EBOS_HEADER_INCLUDED +#define OPM_ISTLSOLVERCPR_EBOS_HEADER_INCLUDED + +#include +#include +#include +#include + +namespace Opm +{ +//===================================================================== +// Implementation for ISTL-matrix based operator +//===================================================================== + + + /// This class solves the fully implicit black-oil system by + /// solving the reduced system (after eliminating well variables) + /// as a block-structured matrix (one block for all cell variables) for a fixed + /// number of cell variables np . + /// \tparam MatrixBlockType The type of the matrix block used. + /// \tparam VectorBlockType The type of the vector block used. + /// \tparam pressureIndex The index of the pressure component in the vector + /// vector block. It is used to guide the AMG coarsening. + /// Default is zero. + template + class ISTLSolverEbosCpr : public ISTLSolverEbos + { + protected: + // Types and indices from superclass. + using SuperClass = ISTLSolverEbos; + using Matrix = typename SuperClass::Matrix; + using Vector = typename SuperClass::Vector; + using WellModel = typename SuperClass::WellModel; + using Simulator = typename SuperClass::Simulator; + using SparseMatrixAdapter = typename SuperClass::SparseMatrixAdapter; + enum { pressureEqnIndex = SuperClass::pressureEqnIndex }; + enum { pressureVarIndex = SuperClass::pressureVarIndex }; + + // New properties in this subclass. + using Preconditioner = Dune::Preconditioner; + using MatrixAdapter = Dune::MatrixAdapter; + + using CouplingMetric = Opm::Amg::Element; + using CritBase = Dune::Amg::SymmetricCriterion; + using Criterion = Dune::Amg::CoarsenCriterion; + + using ParallelMatrixAdapter = Dune::OverlappingSchwarzOperator >; + using CprSmootherFine = Opm::ParallelOverlappingILU0; + using CprSmootherCoarse = CprSmootherFine; + using BlackoilAmgType = BlackoilAmgCpr; + using ParallelCprSmootherFine = Opm::ParallelOverlappingILU0 >; + using ParallelCprSmootherCoarse = ParallelCprSmootherFine; + using ParallelBlackoilAmgType = BlackoilAmgCpr, pressureEqnIndex, pressureVarIndex>; + + using OperatorSerial = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false>; + using OperatorParallel = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true>; + + public: + static void registerParameters() + { + FlowLinearSolverParameters::registerParameters(); + } + + /// Construct a system solver. + /// \param[in] parallelInformation In the case of a parallel run + /// with dune-istl the information about the parallelization. + explicit ISTLSolverEbosCpr(const Simulator& simulator) + : SuperClass(simulator), oldMat() + { + extractParallelGridInformationToISTL(this->simulator_.vanguard().grid(), this->parallelInformation_); + detail::findOverlapRowsAndColumns(this->simulator_.vanguard().grid(), this->overlapRowAndColumns_); + } + + void prepare(const SparseMatrixAdapter& M, Vector& b) + { + if (oldMat != nullptr) + std::cout << "old was "<simulator_.model().newtonMethod().numIterations(); + if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) { + SuperClass::matrix_.reset(new Matrix(M.istlMatrix())); + } else { + *SuperClass::matrix_ = M.istlMatrix(); + } + SuperClass::rhs_ = &b; + SuperClass::scaleSystem(); + const WellModel& wellModel = this->simulator_.problem().wellModel(); + +#if HAVE_MPI + if( this->isParallel() ) { + + //remove ghost rows in local matrix without doing a copy. + this->makeOverlapRowsInvalid(*(this->matrix_)); + + if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) { + //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables + //to be certain that correct matrix is used for preconditioning. + if( ! comm_ ) + { + opAParallel_.reset(new OperatorParallel(*(this->matrix_), *(this->matrix_), wellModel, + this->parallelInformation_ )); + comm_ = opAParallel_->comm(); + assert(comm_->indexSet().size()==0); + const size_t size = opAParallel_->getmat().N(); + + const ParallelISTLInformation& info = + boost::any_cast( this->parallelInformation_); + + // As we use a dune-istl with block size np the number of components + // per parallel is only one. + info.copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), + size, 1); + } + else + { + opAParallel_.reset(new OperatorParallel(*(this->matrix_), *(this->matrix_), wellModel, + comm_ )); + } + } + + using POrComm = Dune::OwnerOverlapCopyCommunication; + +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::overlapping; + auto sp = Dune::createScalarProduct(*comm_, category); + sp_ = std::move(sp); +#else + constexpr int category = Dune::SolverCategory::overlapping; + typedef Dune::ScalarProductChooser ScalarProductChooser; + typedef std::unique_ptr SPPointer; + SPPointer sp(ScalarProductChooser::construct(*comm_)); + sp_ = std::move(sp); +#endif + + using AMGOperator = Dune::OverlappingSchwarzOperator; + // If clause is always execute as as Linearoperator is WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false|true>; + if( ! std::is_same< OperatorParallel, AMGOperator > :: value && + ( newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_) ) ) { + // create new operator in case linear operator and matrix operator differ + opA_.reset( new AMGOperator( opAParallel_->getmat(), *comm_ )); + } + + prepareSolver(*opAParallel_, *comm_); + + } else +#endif + { + + if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) { + opASerial_.reset(new OperatorSerial(*(this->matrix_), *(this->matrix_), wellModel)); + } + + using POrComm = Dune::Amg::SequentialInformation; + POrComm parallelInformation_arg; + typedef OperatorSerial LinearOperator; + +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::sequential; + auto sp = Dune::createScalarProduct(parallelInformation_arg, category); + sp_ = std::move(sp); +#else + constexpr int category = Dune::SolverCategory::sequential; + typedef Dune::ScalarProductChooser ScalarProductChooser; + typedef std::unique_ptr SPPointer; + SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg)); + sp_ = std::move(sp); +#endif + + // If clause is always execute as as Linearoperator is WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false|true>; + if( ! std::is_same< LinearOperator, MatrixAdapter > :: value && + ( newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_) ) ) { + // create new operator in case linear operator and matrix operator differ + opA_.reset( new MatrixAdapter( opASerial_->getmat()));//, parallelInformation_arg ) ); + } + + prepareSolver(*opASerial_, parallelInformation_arg); + } + } + + template + void prepareSolver(Operator& wellOpA, Comm& comm) + { + + Vector& istlb = *(this->rhs_); + comm.copyOwnerToAll(istlb, istlb); + + const double relax = this->parameters_.ilu_relaxation_; + const MILU_VARIANT ilu_milu = this->parameters_.ilu_milu_; + using Matrix = typename MatrixAdapter::matrix_type; + const int verbosity = ( this->parameters_.cpr_solver_verbose_ && + comm.communicator().rank()==0 ) ? 1 : 0; + + // TODO: revise choice of parameters + // int coarsenTarget = 4000; + int coarsenTarget = 1200; + Criterion criterion(15, coarsenTarget); + criterion.setDebugLevel( this->parameters_.cpr_solver_verbose_ ); // no debug information, 1 for printing hierarchy information + criterion.setDefaultValuesIsotropic(2); + criterion.setNoPostSmoothSteps( 1 ); + criterion.setNoPreSmoothSteps( 1 ); + //new guesses by hmbn + //criterion.setAlpha(0.01); // criterion for connection strong 1/3 is default + //criterion.setMaxLevel(2); // + //criterion.setGamma(1); // //1 V cycle 2 WW + + // Since DUNE 2.2 we also need to pass the smoother args instead of steps directly + using AmgType = typename std::conditional::value, + BlackoilAmgType, ParallelBlackoilAmgType>::type; + using SpType = typename std::conditional::value, + Dune::SeqScalarProduct, + Dune::OverlappingSchwarzScalarProduct >::type; + using OperatorType = typename std::conditional::value, + MatrixAdapter, ParallelMatrixAdapter>::type; + typedef typename AmgType::Smoother Smoother; + typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; + SmootherArgs smootherArgs; + smootherArgs.iterations = 1; + smootherArgs.relaxationFactor = relax; + const Opm::CPRParameter& params(this->parameters_); // strange conversion + ISTLUtility::setILUParameters(smootherArgs, ilu_milu); + + auto& opARef = reinterpret_cast(*opA_); + int newton_iteration = this->simulator_.model().newtonMethod().numIterations(); + double dt = this->simulator_.timeStepSize(); + bool update_preconditioner = false; + + if (this->parameters_.cpr_reuse_setup_ < 1) { + update_preconditioner = true; + } + if (this->parameters_.cpr_reuse_setup_ < 2) { + if (newton_iteration < 1) { + update_preconditioner = true; + } + } + if (this->parameters_.cpr_reuse_setup_ < 3) { + if (this->iterations() > 10) { + update_preconditioner = true; + } + } + + if ( update_preconditioner or (amg_== 0) ) { + amg_.reset( new AmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) ); + } else { + if (this->parameters_.cpr_solver_verbose_) { + std::cout << " Only update amg solver " << std::endl; + } + reinterpret_cast(amg_.get())->updatePreconditioner(opARef, smootherArgs, comm); + } + // Solve. + //SuperClass::solve(linearOperator, x, istlb, *sp, *amg, result); + //references seems to do something els than refering + + int verbosity_linsolve = 0; + if (comm.communicator().rank() == 0) { + verbosity_linsolve = this->parameters_.linear_solver_verbosity_; + } + + linsolve_.reset(new Dune::BiCGSTABSolver(wellOpA, reinterpret_cast(*sp_), reinterpret_cast(*amg_), + this->parameters_.linear_solver_reduction_, + this->parameters_.linear_solver_maxiter_, + verbosity_linsolve)); + } + + bool solve(Vector& x) + { + // Solve system. + Dune::InverseOperatorResult result; + Vector& istlb = *(this->rhs_); + linsolve_->apply(x, istlb, result); + SuperClass::checkConvergence(result); + if (this->parameters_.scale_linear_system_) { + this->scaleSolution(x); + } + return this->converged_; + } + + + protected: + + ///! \brief The dune-istl operator (either serial or parallel + std::unique_ptr< Dune::LinearOperator > opA_; + ///! \brief Serial well matrix adapter + std::unique_ptr< OperatorSerial > opASerial_; + ///! \brief Parallel well matrix adapter + std::unique_ptr< OperatorParallel > opAParallel_; + ///! \brief The preconditoner to use (either serial or parallel CPR with AMG) + std::unique_ptr< Preconditioner > amg_; + + using SPPointer = std::shared_ptr< Dune::ScalarProduct >; + SPPointer sp_; + std::shared_ptr< Dune::BiCGSTABSolver > linsolve_; + const void* oldMat; + using POrComm = Dune::OwnerOverlapCopyCommunication; + std::shared_ptr comm_; + }; // end ISTLSolver + +} // namespace Opm +#endif diff --git a/opm/autodiff/ParallelOverlappingILU0.hpp b/opm/autodiff/ParallelOverlappingILU0.hpp index 2698841ef..d65a41900 100644 --- a/opm/autodiff/ParallelOverlappingILU0.hpp +++ b/opm/autodiff/ParallelOverlappingILU0.hpp @@ -724,7 +724,7 @@ public: \copydoc Preconditioner::pre(X&,Y&) */ - virtual void pre (Domain& x, Range& b) + virtual void pre (Domain& x, Range& b) override { DUNE_UNUSED_PARAMETER(x); DUNE_UNUSED_PARAMETER(b); @@ -735,7 +735,7 @@ public: \copydoc Preconditioner::apply(X&,const Y&) */ - virtual void apply (Domain& v, const Range& d) + virtual void apply (Domain& v, const Range& d) override { Range& md = reorderD(d); Domain& mv = reorderV(v); @@ -806,7 +806,7 @@ public: \copydoc Preconditioner::post(X&) */ - virtual void post (Range& x) + virtual void post (Range& x) override { DUNE_UNUSED_PARAMETER(x); } diff --git a/opm/autodiff/amgcpr.hh b/opm/autodiff/amgcpr.hh new file mode 100644 index 000000000..838b220b1 --- /dev/null +++ b/opm/autodiff/amgcpr.hh @@ -0,0 +1,938 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef DUNE_AMG_AMG_CPR_HH +#define DUNE_AMG_AMG_CPR_HH + +// NOTE: This file is a modified version of dune/istl/paamg/amg.hh from +// dune-istl release 2.6.0. Modifications have been kept as minimal as possible. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Dune +{ + namespace Amg + { + /** + * @defgroup ISTL_PAAMG Parallel Algebraic Multigrid + * @ingroup ISTL_Prec + * @brief A Parallel Algebraic Multigrid based on Agglomeration. + */ + + /** + * @addtogroup ISTL_PAAMG + * + * @{ + */ + + /** @file + * @author Markus Blatt + * @brief The AMG preconditioner. + */ + + template + class KAMG; + + template + class KAmgTwoGrid; + + /** + * @brief Parallel algebraic multigrid based on agglomeration. + * + * \tparam M The matrix type + * \tparam X The vector type + * \tparam S The smoother type + * \tparam A An allocator for X + * + * \todo drop the smoother template parameter and replace with dynamic construction + */ + template > + class AMGCPR : public Preconditioner + { + template + friend class KAMG; + + friend class KAmgTwoGrid; + + public: + /** @brief The matrix operator type. */ + typedef M Operator; + /** + * @brief The type of the parallel information. + * Either OwnerOverlapCommunication or another type + * describing the parallel data distribution and + * providing communication methods. + */ + typedef PI ParallelInformation; + /** @brief The operator hierarchy type. */ + typedef MatrixHierarchy OperatorHierarchy; + /** @brief The parallal data distribution hierarchy type. */ + typedef typename OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy; + + /** @brief The domain type. */ + typedef X Domain; + /** @brief The range type. */ + typedef X Range; + /** @brief the type of the coarse solver. */ + typedef InverseOperator CoarseSolver; + /** + * @brief The type of the smoother. + * + * One of the preconditioners implementing the Preconditioner interface. + * Note that the smoother has to fit the ParallelInformation.*/ + typedef S Smoother; + + /** @brief The argument type for the construction of the smoother. */ + typedef typename SmootherTraits::Arguments SmootherArgs; + + /** + * @brief Construct a new amg with a specific coarse solver. + * @param matrices The already set up matix hierarchy. + * @param coarseSolver The set up solver to use on the coarse + * grid, must match the coarse matrix in the matrix hierarchy. + * @param smootherArgs The arguments needed for thesmoother to use + * for pre and post smoothing. + * @param parms The parameters for the AMG. + */ + AMGCPR(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, + const SmootherArgs& smootherArgs, const Parameters& parms); + + /** + * @brief Construct an AMG with an inexact coarse solver based on the smoother. + * + * As coarse solver a preconditioned CG method with the smoother as preconditioner + * will be used. The matrix hierarchy is built automatically. + * @param fineOperator The operator on the fine level. + * @param criterion The criterion describing the coarsening strategy. E. g. SymmetricCriterion + * or UnsymmetricCriterion, and providing the parameters. + * @param smootherArgs The arguments for constructing the smoothers. + * @param pinfo The information about the parallel distribution of the data. + */ + template + AMGCPR(const Operator& fineOperator, const C& criterion, + const SmootherArgs& smootherArgs=SmootherArgs(), + const ParallelInformation& pinfo=ParallelInformation()); + + /** + * @brief Copy constructor. + */ + AMGCPR(const AMGCPR& amg); + + ~AMGCPR(); + + /** \copydoc Preconditioner::pre */ + void pre(Domain& x, Range& b); + + /** \copydoc Preconditioner::apply */ + void apply(Domain& v, const Range& d); + +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + //! Category of the preconditioner (see SolverCategory::Category) + virtual SolverCategory::Category category() const + { + return category_; + } +#else + enum { + //! \brief The category the preconditioner is part of. + category = std::is_same::value? + Dune::SolverCategory::sequential:Dune::SolverCategory::overlapping + }; +#endif + + /** \copydoc Preconditioner::post */ + void post(Domain& x); + + /** + * @brief Get the aggregate number of each unknown on the coarsest level. + * @param cont The random access container to store the numbers in. + */ + template + void getCoarsestAggregateNumbers(std::vector& cont); + + std::size_t levels(); + + std::size_t maxlevels(); + + /** + * @brief Recalculate the matrix hierarchy. + * + * It is assumed that the coarsening for the changed fine level + * matrix would yield the same aggregates. In this case it suffices + * to recalculate all the Galerkin products for the matrices of the + * coarser levels. + */ + void recalculateHierarchy() + { + matrices_->recalculateGalerkin(NegateSet()); + } + + /** + * @brief Update the coarse solver and the hierarchies. + */ + template + void updateSolver(C& criterion, Operator& /* matrix */, const PI& pinfo); + + /** + * @brief Check whether the coarse solver used is a direct solver. + * @return True if the coarse level solver is a direct solver. + */ + bool usesDirectCoarseLevelSolver() const; + + private: + /** + * @brief Create matrix and smoother hierarchies. + * @param criterion The coarsening criterion. + * @param matrix The fine level matrix operator. + * @param pinfo The fine level parallel information. + */ + template + void createHierarchies(C& criterion, Operator& matrix, + const PI& pinfo); + + void setupCoarseSolver(); + + /** + * @brief A struct that holds the context of the current level. + * + * These are the iterators to the smoother, matrix, parallel information, + * and so on needed for the computations on the current level. + */ + struct LevelContext + { + typedef Smoother SmootherType; + /** + * @brief The iterator over the smoothers. + */ + typename Hierarchy::Iterator smoother; + /** + * @brief The iterator over the matrices. + */ + typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix; + /** + * @brief The iterator over the parallel information. + */ + typename ParallelInformationHierarchy::Iterator pinfo; + /** + * @brief The iterator over the redistribution information. + */ + typename OperatorHierarchy::RedistributeInfoList::const_iterator redist; + /** + * @brief The iterator over the aggregates maps. + */ + typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates; + /** + * @brief The iterator over the left hand side. + */ + typename Hierarchy::Iterator lhs; + /** + * @brief The iterator over the updates. + */ + typename Hierarchy::Iterator update; + /** + * @brief The iterator over the right hand sided. + */ + typename Hierarchy::Iterator rhs; + /** + * @brief The level index. + */ + std::size_t level; + }; + + + /** + * @brief Multigrid cycle on a level. + * @param levelContext the iterators of the current level. + */ + void mgc(LevelContext& levelContext); + + void additiveMgc(); + + /** + * @brief Move the iterators to the finer level + * @param levelContext the iterators of the current level. + * @param processedFineLevel Whether the process computed on + * fine level or not. + */ + void moveToFineLevel(LevelContext& levelContext,bool processedFineLevel); + + /** + * @brief Move the iterators to the coarser level. + * @param levelContext the iterators of the current level + */ + bool moveToCoarseLevel(LevelContext& levelContext); + + /** + * @brief Initialize iterators over levels with fine level. + * @param levelContext the iterators of the current level + */ + void initIteratorsWithFineLevel(LevelContext& levelContext); + + /** @brief The matrix we solve. */ + std::shared_ptr matrices_; + /** @brief The arguments to construct the smoother */ + SmootherArgs smootherArgs_; + /** @brief The hierarchy of the smoothers. */ + std::shared_ptr > smoothers_; + /** @brief The solver of the coarsest level. */ + std::shared_ptr solver_; + /** @brief The right hand side of our problem. */ + Hierarchy* rhs_; + /** @brief The left approximate solution of our problem. */ + Hierarchy* lhs_; + /** @brief The total update for the outer solver. */ + Hierarchy* update_; + /** @brief The type of the scalar product for the coarse solver. */ + using ScalarProduct = Dune::ScalarProduct; + /** @brief Scalar product on the coarse level. */ + std::shared_ptr scalarProduct_; + /** @brief Gamma, 1 for V-cycle and 2 for W-cycle. */ + std::size_t gamma_; + /** @brief The number of pre and postsmoothing steps. */ + std::size_t preSteps_; + /** @brief The number of postsmoothing steps. */ + std::size_t postSteps_; + bool buildHierarchy_; + bool additive; + bool coarsesolverconverged; + std::shared_ptr coarseSmoother_; +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + /** @brief The solver category. */ + SolverCategory::Category category_; +#endif + /** @brief The verbosity level. */ + std::size_t verbosity_; + }; + + template + inline AMGCPR::AMGCPR(const AMGCPR& amg) + : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_), + smoothers_(amg.smoothers_), solver_(amg.solver_), + rhs_(), lhs_(), update_(), + scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_), + preSteps_(amg.preSteps_), postSteps_(amg.postSteps_), + buildHierarchy_(amg.buildHierarchy_), + additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged), + coarseSmoother_(amg.coarseSmoother_), +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + category_(amg.category_), +#endif + verbosity_(amg.verbosity_) + { + if(amg.rhs_) + rhs_=new Hierarchy(*amg.rhs_); + if(amg.lhs_) + lhs_=new Hierarchy(*amg.lhs_); + if(amg.update_) + update_=new Hierarchy(*amg.update_); + } + + template + AMGCPR::AMGCPR(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, + const SmootherArgs& smootherArgs, + const Parameters& parms) + : matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs), + smoothers_(new Hierarchy), solver_(&coarseSolver), + rhs_(), lhs_(), update_(), scalarProduct_(0), + gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()), + postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false), + additive(parms.getAdditive()), coarsesolverconverged(true), + coarseSmoother_(), +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) +// #warning should category be retrieved from matrices? + category_(SolverCategory::category(*smoothers_->coarsest())), +#endif + verbosity_(parms.debugLevel()) + { + assert(matrices_->isBuilt()); + + // build the necessary smoother hierarchies + matrices_->coarsenSmoother(*smoothers_, smootherArgs_); + } + + template + template + AMGCPR::AMGCPR(const Operator& matrix, + const C& criterion, + const SmootherArgs& smootherArgs, + const PI& pinfo) + : smootherArgs_(smootherArgs), + smoothers_(new Hierarchy), solver_(), + rhs_(), lhs_(), update_(), scalarProduct_(), + gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()), + postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true), + additive(criterion.getAdditive()), coarsesolverconverged(true), + coarseSmoother_(), +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + category_(SolverCategory::category(pinfo)), +#endif + verbosity_(criterion.debugLevel()) + { +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + if(SolverCategory::category(matrix) != SolverCategory::category(pinfo)) + DUNE_THROW(InvalidSolverCategory, "Matrix and Communication must have the same SolverCategory!"); +#endif + createHierarchies(criterion, const_cast(matrix), pinfo); + } + + + template + AMGCPR::~AMGCPR() + { + if(buildHierarchy_) { + if(solver_) + solver_.reset(); + if(coarseSmoother_) + coarseSmoother_.reset(); + } + if(lhs_) + delete lhs_; + lhs_=nullptr; + if(update_) + delete update_; + update_=nullptr; + if(rhs_) + delete rhs_; + rhs_=nullptr; + } + + template + template + void AMGCPR::updateSolver(C& /* criterion */, Operator& /* matrix */, const PI& /* pinfo */) + { + Timer watch; + smoothers_.reset(new Hierarchy); + solver_.reset(); + coarseSmoother_.reset(); + scalarProduct_.reset(); + buildHierarchy_= true; + coarsesolverconverged = true; + smoothers_.reset(new Hierarchy); + matrices_->recalculateGalerkin(NegateSet()); + matrices_->coarsenSmoother(*smoothers_, smootherArgs_); + setupCoarseSolver(); + if (verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) { + std::cout << "Recalculating galerkin and coarse somothers "<< matrices_->maxlevels() << " levels " + << watch.elapsed() << " seconds." << std::endl; + } + } + + template + template + void AMGCPR::createHierarchies(C& criterion, Operator& matrix, + const PI& pinfo) + { + Timer watch; + matrices_.reset(new OperatorHierarchy(matrix, pinfo)); + + matrices_->template build >(criterion); + + // build the necessary smoother hierarchies + matrices_->coarsenSmoother(*smoothers_, smootherArgs_); + setupCoarseSolver(); + if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) + std::cout<<"Building hierarchy of "<maxlevels()<<" levels " + <<"(inclusive coarse solver) took "< + void AMGCPR::setupCoarseSolver() + { + // test whether we should solve on the coarse level. That is the case if we + // have that level and if there was a redistribution on this level then our + // communicator has to be valid (size()>0) as the smoother might try to communicate + // in the constructor. + if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels() + && ( ! matrices_->redistributeInformation().back().isSetup() || + matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) ) + { + // We have the carsest level. Create the coarse Solver + SmootherArgs sargs(smootherArgs_); + sargs.iterations = 1; + + typename ConstructionTraits::Arguments cargs; + cargs.setArgs(sargs); + if(matrices_->redistributeInformation().back().isSetup()) { + // Solve on the redistributed partitioning + cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat()); + cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed()); + }else{ + cargs.setMatrix(matrices_->matrices().coarsest()->getmat()); + cargs.setComm(*matrices_->parallelInformation().coarsest()); + } + + coarseSmoother_.reset(ConstructionTraits::construct(cargs)); + +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + scalarProduct_ = createScalarProduct(cargs.getComm(),category()); +#else + typedef Dune::ScalarProductChooser + ScalarProductChooser; + // the scalar product. + scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm())); +#endif + + + typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector; + + // Use superlu if we are purely sequential or with only one processor on the coarsest level. + if( SolverSelector::isDirectSolver && + (std::is_same::value // sequential mode + || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor + || (matrices_->parallelInformation().coarsest().isRedistributed() + && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1 + && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) ) + ) + { // redistribute and 1 proc + if(matrices_->parallelInformation().coarsest().isRedistributed()) + { + if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0) + { + // We are still participating on this level + solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false)); + } + else + solver_.reset(); + } + else + { + solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false)); + } + if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0) + std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl; + } + else + { + if(matrices_->parallelInformation().coarsest().isRedistributed()) + { + if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0) + // We are still participating on this level + + // we have to allocate these types using the rebound allocator + // in order to ensure that we fulfill the alignement requirements + solver_.reset(new BiCGSTABSolver(const_cast(matrices_->matrices().coarsest().getRedistributed()), + // Cast needed for Dune <=2.5 + reinterpret_cast::value, + Dune::SeqScalarProduct, + Dune::OverlappingSchwarzScalarProduct >::type&>(*scalarProduct_), + *coarseSmoother_, 1E-2, 1000, 0)); + else + solver_.reset(); + }else + { + solver_.reset(new BiCGSTABSolver(const_cast(*matrices_->matrices().coarsest()), + // Cast needed for Dune <=2.5 + reinterpret_cast::value, + Dune::SeqScalarProduct, + Dune::OverlappingSchwarzScalarProduct >::type&>(*scalarProduct_), + *coarseSmoother_, 1E-2, 1000, 0)); + // // we have to allocate these types using the rebound allocator + // // in order to ensure that we fulfill the alignement requirements + // using Alloc = typename A::template rebind>::other; + // Alloc alloc; + // auto p = alloc.allocate(1); + // alloc.construct(p, + // const_cast(*matrices_->matrices().coarsest()), + // *scalarProduct_, + // *coarseSmoother_, 1E-2, 1000, 0); + // solver_.reset(p,[](BiCGSTABSolver* p){ + // Alloc alloc; + // alloc.destroy(p); + // alloc.deallocate(p,1); + // }); + } + } + } + } + + template + void AMGCPR::pre(Domain& x, Range& b) + { + // Detect Matrix rows where all offdiagonal entries are + // zero and set x such that A_dd*x_d=b_d + // Thus users can be more careless when setting up their linear + // systems. + typedef typename M::matrix_type Matrix; + typedef typename Matrix::ConstRowIterator RowIter; + typedef typename Matrix::ConstColIterator ColIter; + typedef typename Matrix::block_type Block; + Block zero; + zero=typename Matrix::field_type(); + + const Matrix& mat=matrices_->matrices().finest()->getmat(); + for(RowIter row=mat.begin(); row!=mat.end(); ++row) { + bool isDirichlet = true; + bool hasDiagonal = false; + Block diagonal; + for(ColIter col=row->begin(); col!=row->end(); ++col) { + if(row.index()==col.index()) { + diagonal = *col; + hasDiagonal = false; + }else{ + if(*col!=zero) + isDirichlet = false; + } + } + if(isDirichlet && hasDiagonal) + diagonal.solve(x[row.index()], b[row.index()]); + } + + if(smoothers_->levels()>0) + smoothers_->finest()->pre(x,b); + else + // No smoother to make x consistent! Do it by hand + matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x); + Range* copy = new Range(b); + if(rhs_) + delete rhs_; + rhs_ = new Hierarchy(copy); + Domain* dcopy = new Domain(x); + if(lhs_) + delete lhs_; + lhs_ = new Hierarchy(dcopy); + dcopy = new Domain(x); + if(update_) + delete update_; + update_ = new Hierarchy(dcopy); + matrices_->coarsenVector(*rhs_); + matrices_->coarsenVector(*lhs_); + matrices_->coarsenVector(*update_); + + // Preprocess all smoothers + typedef typename Hierarchy::Iterator Iterator; + typedef typename Hierarchy::Iterator RIterator; + typedef typename Hierarchy::Iterator DIterator; + Iterator coarsest = smoothers_->coarsest(); + Iterator smoother = smoothers_->finest(); + RIterator rhs = rhs_->finest(); + DIterator lhs = lhs_->finest(); + if(smoothers_->levels()>0) { + + assert(lhs_->levels()==rhs_->levels()); + assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels()); + assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()maxlevels()); + + if(smoother!=coarsest) + for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs) + smoother->pre(*lhs,*rhs); + smoother->pre(*lhs,*rhs); + } + + + // The preconditioner might change x and b. So we have to + // copy the changes to the original vectors. + x = *lhs_->finest(); + b = *rhs_->finest(); + + } + template + std::size_t AMGCPR::levels() + { + return matrices_->levels(); + } + template + std::size_t AMGCPR::maxlevels() + { + return matrices_->maxlevels(); + } + + /** \copydoc Preconditioner::apply */ + template + void AMGCPR::apply(Domain& v, const Range& d) + { + LevelContext levelContext; + + if(additive) { + *(rhs_->finest())=d; + additiveMgc(); + v=*lhs_->finest(); + }else{ + // Init all iterators for the current level + initIteratorsWithFineLevel(levelContext); + + + *levelContext.lhs = v; + *levelContext.rhs = d; + *levelContext.update=0; + levelContext.level=0; + + mgc(levelContext); + + if(postSteps_==0||matrices_->maxlevels()==1) + levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update); + + v=*levelContext.update; + } + + } + + template + void AMGCPR::initIteratorsWithFineLevel(LevelContext& levelContext) + { + levelContext.smoother = smoothers_->finest(); + levelContext.matrix = matrices_->matrices().finest(); + levelContext.pinfo = matrices_->parallelInformation().finest(); + levelContext.redist = + matrices_->redistributeInformation().begin(); + levelContext.aggregates = matrices_->aggregatesMaps().begin(); + levelContext.lhs = lhs_->finest(); + levelContext.update = update_->finest(); + levelContext.rhs = rhs_->finest(); + } + + template + bool AMGCPR + ::moveToCoarseLevel(LevelContext& levelContext) + { + + bool processNextLevel=true; + + if(levelContext.redist->isSetup()) { + levelContext.redist->redistribute(static_cast(*levelContext.rhs), + levelContext.rhs.getRedistributed()); + processNextLevel = levelContext.rhs.getRedistributed().size()>0; + if(processNextLevel) { + //restrict defect to coarse level right hand side. + typename Hierarchy::Iterator fineRhs = levelContext.rhs++; + ++levelContext.pinfo; + Transfer + ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs, + static_cast(fineRhs.getRedistributed()), + *levelContext.pinfo); + } + }else{ + //restrict defect to coarse level right hand side. + typename Hierarchy::Iterator fineRhs = levelContext.rhs++; + ++levelContext.pinfo; + Transfer + ::restrictVector(*(*levelContext.aggregates), + *levelContext.rhs, static_cast(*fineRhs), + *levelContext.pinfo); + } + + if(processNextLevel) { + // prepare coarse system + ++levelContext.lhs; + ++levelContext.update; + ++levelContext.matrix; + ++levelContext.level; + ++levelContext.redist; + + if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()maxlevels()) { + // next level is not the globally coarsest one + ++levelContext.smoother; + ++levelContext.aggregates; + } + // prepare the update on the next level + *levelContext.update=0; + } + return processNextLevel; + } + + template + void AMGCPR + ::moveToFineLevel(LevelContext& levelContext, bool processNextLevel) + { + if(processNextLevel) { + if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()maxlevels()) { + // previous level is not the globally coarsest one + --levelContext.smoother; + --levelContext.aggregates; + } + --levelContext.redist; + --levelContext.level; + //prolongate and add the correction (update is in coarse left hand side) + --levelContext.matrix; + + //typename Hierarchy::Iterator coarseLhs = lhs--; + --levelContext.lhs; + --levelContext.pinfo; + } + if(levelContext.redist->isSetup()) { + // Need to redistribute during prolongateVector + levelContext.lhs.getRedistributed()=0; + Transfer + ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs, + levelContext.lhs.getRedistributed(), + matrices_->getProlongationDampingFactor(), + *levelContext.pinfo, *levelContext.redist); + }else{ + *levelContext.lhs=0; + Transfer + ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs, + matrices_->getProlongationDampingFactor(), + *levelContext.pinfo); + } + + + if(processNextLevel) { + --levelContext.update; + --levelContext.rhs; + } + + *levelContext.update += *levelContext.lhs; + } + + template + bool AMGCPR::usesDirectCoarseLevelSolver() const + { + return IsDirectSolver< CoarseSolver>::value; + } + + template + void AMGCPR::mgc(LevelContext& levelContext){ + if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) { + // Solve directly + InverseOperatorResult res; + res.converged=true; // If we do not compute this flag will not get updated + if(levelContext.redist->isSetup()) { + levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed()); + if(levelContext.rhs.getRedistributed().size()>0) { + // We are still participating in the computation + levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(), + levelContext.rhs.getRedistributed()); + solver_->apply(levelContext.update.getRedistributed(), + levelContext.rhs.getRedistributed(), res); + } + levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed()); + levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update); + }else{ + levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs); + solver_->apply(*levelContext.update, *levelContext.rhs, res); + } + + if (!res.converged) + coarsesolverconverged = false; + }else{ + // presmoothing + presmooth(levelContext, preSteps_); + +#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION + bool processNextLevel = moveToCoarseLevel(levelContext); + + if(processNextLevel) { + // next level + for(std::size_t i=0; imatrices().finest()) { + coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged); + if(!coarsesolverconverged){ + //DUNE_THROW(MathError, "Coarse solver did not converge"); + } + } + // postsmoothing + postsmooth(levelContext, postSteps_); + + } + } + + template + void AMGCPR::additiveMgc(){ + + // restrict residual to all levels + typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest(); + typename Hierarchy::Iterator rhs=rhs_->finest(); + typename Hierarchy::Iterator lhs = lhs_->finest(); + typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin(); + + for(typename Hierarchy::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) { + ++pinfo; + Transfer + ::restrictVector(*(*aggregates), *rhs, static_cast(*fineRhs), *pinfo); + } + + // pinfo is invalid, set to coarsest level + //pinfo = matrices_->parallelInformation().coarsest + // calculate correction for all levels + lhs = lhs_->finest(); + typename Hierarchy::Iterator smoother = smoothers_->finest(); + + for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) { + // presmoothing + *lhs=0; + smoother->apply(*lhs, *rhs); + } + + // Coarse level solve +#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION + InverseOperatorResult res; + pinfo->copyOwnerToAll(*rhs, *rhs); + solver_->apply(*lhs, *rhs, res); + + if(!res.converged) + DUNE_THROW(MathError, "Coarse solver did not converge"); +#else + *lhs=0; +#endif + // Prologate and add up corrections from all levels + --pinfo; + --aggregates; + + for(typename Hierarchy::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) { + Transfer + ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo); + } + } + + + /** \copydoc Preconditioner::post */ + template + void AMGCPR::post(Domain& x) + { + DUNE_UNUSED_PARAMETER(x); + // Postprocess all smoothers + typedef typename Hierarchy::Iterator Iterator; + typedef typename Hierarchy::Iterator DIterator; + Iterator coarsest = smoothers_->coarsest(); + Iterator smoother = smoothers_->finest(); + DIterator lhs = lhs_->finest(); + if(smoothers_->levels()>0) { + if(smoother != coarsest || matrices_->levels()maxlevels()) + smoother->post(*lhs); + if(smoother!=coarsest) + for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs) + smoother->post(*lhs); + smoother->post(*lhs); + } + //delete &(*lhs_->finest()); + delete lhs_; + lhs_=nullptr; + //delete &(*update_->finest()); + delete update_; + update_=nullptr; + //delete &(*rhs_->finest()); + delete rhs_; + rhs_=nullptr; + } + + template + template + void AMGCPR::getCoarsestAggregateNumbers(std::vector& cont) + { + matrices_->getCoarsestAggregatesOnFinest(cont); + } + + } // end namespace Amg +} // end namespace Dune + +#endif diff --git a/opm/autodiff/twolevelmethodcpr.hh b/opm/autodiff/twolevelmethodcpr.hh new file mode 100644 index 000000000..4c23ce681 --- /dev/null +++ b/opm/autodiff/twolevelmethodcpr.hh @@ -0,0 +1,562 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#ifndef DUNE_ISTL_TWOLEVELMETHODCPR_HH +#define DUNE_ISTL_TWOLEVELMETHODCPR_HH + +// NOTE: This file is a modified version of dune/istl/paamg/twolevelmethod.hh from +// dune-istl release 2.6.0. Modifications have been kept as minimal as possible. + +#include + +#include +//#include "amg.hh" +//#include"galerkin.hh" +#include +#include +#include + +#include + +/** + * @addtogroup ISTL_PAAMG + * @{ + * @file + * @author Markus Blatt + * @brief Algebraic twolevel methods. + */ +namespace Dune +{ +namespace Amg +{ + +/** + * @brief Abstract base class for transfer between levels and creation + * of the coarse level system. + * + * @tparam FO The type of the linear operator of the finel level system. Has to be + * derived from AssembledLinearOperator. + * @tparam CO The type of the linear operator of the coarse level system. Has to be + * derived from AssembledLinearOperator. + */ +template +class LevelTransferPolicyCpr +{ +public: + /** + * @brief The linear operator of the finel level system. Has to be + * derived from AssembledLinearOperator. + */ + typedef FO FineOperatorType; + /** + * @brief The type of the range of the fine level operator. + */ + typedef typename FineOperatorType::range_type FineRangeType; + /** + * @brief The type of the domain of the fine level operator. + */ + typedef typename FineOperatorType::domain_type FineDomainType; + /** + * @brief The linear operator of the finel level system. Has to be + * derived from AssembledLinearOperator. + */ + typedef CO CoarseOperatorType; + /** + * @brief The type of the range of the coarse level operator. + */ + typedef typename CoarseOperatorType::range_type CoarseRangeType; + /** + * @brief The type of the domain of the coarse level operator. + */ + typedef typename CoarseOperatorType::domain_type CoarseDomainType; + /** + * @brief Get the coarse level operator. + * @return A shared pointer to the coarse level system. + */ + std::shared_ptr& getCoarseLevelOperator() + { + return operator_; + } + /** + * @brief Get the coarse level right hand side. + * @return The coarse level right hand side. + */ + CoarseRangeType& getCoarseLevelRhs() + { + return rhs_; + } + + /** + * @brief Get the coarse level left hand side. + * @return The coarse level leftt hand side. + */ + CoarseDomainType& getCoarseLevelLhs() + { + return lhs_; + } + /** + * @brief Transfers the data to the coarse level. + * + * Restricts the residual to the right hand side of the + * coarse level system and initialies the left hand side + * of the coarse level system. These can afterwards be accessed + * usinf getCoarseLevelRhs() and getCoarseLevelLhs(). + * @param fineDefect The current residual of the fine level system. + */ + virtual void moveToCoarseLevel(const FineRangeType& fineRhs)=0; + /** + * @brief Updates the fine level linear system after the correction + * of the coarse levels system. + * + * After returning from this function the coarse level correction + * will have been added to fine level system. + * @param[inout] fineLhs The left hand side of the fine level to update + * with the coarse level correction. + */ + virtual void moveToFineLevel(FineDomainType& fineLhs)=0; + /** + * @brief Algebraically creates the coarse level system. + * + * After returning from this function the coarse level operator + * can be accessed using getCoarseLevelOperator(). + * @param fineOperator The operator of the fine level system. + */ + virtual void createCoarseLevelSystem(const FineOperatorType& fineOperator)=0; + + /** + * @brief ???. + */ + virtual void calculateCoarseEntries(const FineOperatorType& fineOperator) = 0; + + /** @brief Clone the current object. */ + virtual LevelTransferPolicyCpr* clone() const =0; + + /** @brief Destructor. */ + virtual ~LevelTransferPolicyCpr(){} + + protected: + /** @brief The coarse level rhs. */ + CoarseRangeType rhs_; + /** @brief The coarse level lhs. */ + CoarseDomainType lhs_; + /** @brief the coarse level linear operator. */ + std::shared_ptr operator_; +}; + +/** + * @brief A LeveTransferPolicy that used aggregation to construct the coarse level system. + * @tparam O The type of the fine and coarse level operator. + * @tparam C The criterion that describes the aggregation procedure. + */ +template +class AggregationLevelTransferPolicyCpr + : public LevelTransferPolicyCpr +{ + typedef Dune::Amg::AggregatesMap AggregatesMap; +public: + typedef LevelTransferPolicyCpr FatherType; + typedef C Criterion; + typedef SequentialInformation ParallelInformation; + + AggregationLevelTransferPolicyCpr(const Criterion& crit) + : criterion_(crit) + {} + + void createCoarseLevelSystem(const O& fineOperator) + { + prolongDamp_ = criterion_.getProlongationDampingFactor(); + GalerkinProduct productBuilder; + typedef typename Dune::Amg::MatrixGraph MatrixGraph; + typedef typename Dune::Amg::PropertiesGraph PropertiesGraph; + MatrixGraph mg(fineOperator.getmat()); + PropertiesGraph pg(mg,Dune::IdentityMap(),Dune::IdentityMap()); + typedef NegateSet OverlapFlags; + + aggregatesMap_.reset(new AggregatesMap(pg.maxVertex()+1)); + + int noAggregates, isoAggregates, oneAggregates, skippedAggregates; + + std::tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) = + aggregatesMap_->buildAggregates(fineOperator.getmat(), pg, criterion_, true); + std::cout<<"no aggregates="<linear_solve_setup_time << "; " + << 100*failureReport->linear_solve_setup_time/t << "%)"; + } + os << std::endl; + t = update_time + (failureReport ? failureReport->update_time : 0.0); os << " Update time (seconds): " << t; if (failureReport) { diff --git a/opm/core/simulator/SimulatorReport.hpp b/opm/core/simulator/SimulatorReport.hpp index e86f4e7ef..e1ffae61a 100644 --- a/opm/core/simulator/SimulatorReport.hpp +++ b/opm/core/simulator/SimulatorReport.hpp @@ -33,6 +33,7 @@ namespace Opm double total_time; double solver_time; double assemble_time; + double linear_solve_setup_time; double linear_solve_time; double update_time; double output_write_time;