From f05a9fdb2526800334bd9a2eea0ee57402a27a5e Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 21 Mar 2019 22:15:22 +0100 Subject: [PATCH 01/25] Version of cpr amg which can reuse setup and also change smoothers of fine and coarse system by changing tags --- CMakeLists.txt | 10 + CMakeLists_files.cmake | 4 + flow/flow_blackoil_dunecpr.cpp | 125 ++++ flow/flow_tag.hpp | 212 ++++++ opm/autodiff/BlackoilAmgCpr.hpp | 597 ++++++++++++++++ opm/autodiff/ISTLSolverEbosCpr.hpp | 350 ++++++++++ opm/autodiff/amgcpr.hh | 1027 ++++++++++++++++++++++++++++ opm/autodiff/twolevelmethodcpr.hh | 567 +++++++++++++++ 8 files changed, 2892 insertions(+) create mode 100644 flow/flow_blackoil_dunecpr.cpp create mode 100644 flow/flow_tag.hpp create mode 100644 opm/autodiff/BlackoilAmgCpr.hpp create mode 100644 opm/autodiff/ISTLSolverEbosCpr.hpp create mode 100644 opm/autodiff/amgcpr.hh create mode 100644 opm/autodiff/twolevelmethodcpr.hh 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..2c8468078 --- /dev/null +++ b/flow/flow_blackoil_dunecpr.cpp @@ -0,0 +1,125 @@ +/* + 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); +NEW_PROP_TAG(CprSmootherFine); +NEW_PROP_TAG(CprSmootherCoarse); +//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_PROP(EclFlowProblemSimple, CprSmootherFine) + { + private: + typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; + typedef typename SparseMatrixAdapter::IstlMatrix Matrix; + typedef Dune::Amg::SequentialInformation POrComm; + + public: + typedef Opm::ParallelOverlappingILU0 type; + //typedef Dune::SeqILU0 type; +}; + +SET_PROP(EclFlowProblemSimple, CprSmootherCoarse) + { + private: + typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; + typedef typename SparseMatrixAdapter::IstlMatrix Matrix; + typedef Dune::Amg::SequentialInformation POrComm; + + public: + typedef Opm::ParallelOverlappingILU0 type; + //typedef Dune::SeqILU0 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..c86c30fa5 --- /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); + 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/BlackoilAmgCpr.hpp b/opm/autodiff/BlackoilAmgCpr.hpp new file mode 100644 index 000000000..55436c00e --- /dev/null +++ b/opm/autodiff/BlackoilAmgCpr.hpp @@ -0,0 +1,597 @@ +/* + 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 Dune +{ + namespace Amg + { + template + class UnSymmetricCriterion; + } +} + +namespace Dune +{ + + template + class MatrixBlock; + +} + +namespace Opm +{ + namespace Detail + { + template + std::unique_ptr scaleMatrixDRSPtr(const Operator& op, + const Communication& comm, + std::size_t pressureEqnIndex, + const Vector& weights, + const Opm::CPRParameter& param) + { + using Matrix = typename Operator::matrix_type; + using Block = typename Matrix::block_type; + using BlockVector = typename Vector::block_type; + std::unique_ptr matrix(new Matrix(op.getmat())); + if (param.cpr_use_drs_) { + const auto endi = matrix->end(); + for (auto i = matrix->begin(); i != endi; ++i) { + const BlockVector& bw = weights[i.index()]; + const auto endj = (*i).end(); + for (auto j = (*i).begin(); j != endj; ++j) { + Block& block = *j; + BlockVector& bvec = block[pressureEqnIndex]; + // should introduce limits which also change the weights + block.mtv(bw, bvec); + } + } + } + return matrix;//, createOperator(op, *matrix, comm)); + } + } + + template + class OneComponentAggregationLevelTransferPolicyCpr; + + template + class OneComponentAggregationLevelTransferPolicyCpr + : public Dune::Amg::LevelTransferPolicyCpr::value> + { + typedef Dune::Amg::AggregatesMap AggregatesMap; + public: + using CoarseOperator = typename Detail::ScalarType::value; + typedef Dune::Amg::LevelTransferPolicy FatherType; + typedef Communication ParallelInformation; + + public: + OneComponentAggregationLevelTransferPolicyCpr(const Criterion& crit, const Communication& comm) + : criterion_(crit), communication_(&const_cast(comm)) + {} + + void createCoarseLevelSystem(const Operator& fineOperator) + { + prolongDamp_ = 1; + + using CoarseMatrix = typename CoarseOperator::matrix_type; + const auto& fineLevelMatrix = fineOperator.getmat(); + coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise)); + auto createIter = coarseLevelMatrix_->createbegin(); + + for ( const auto& row: fineLevelMatrix ) + { + for ( auto col = row.begin(), cend = row.end(); col != cend; ++col) + { + createIter.insert(col.index()); + } + ++createIter; + } + + auto coarseRow = coarseLevelMatrix_->begin(); + for ( const auto& row: fineLevelMatrix ) + { + auto coarseCol = coarseRow->begin(); + + for ( auto col = row.begin(), cend = row.end(); col != cend; ++col, ++coarseCol ) + { + assert( col.index() == coarseCol.index() ); + *coarseCol = (*col)[COMPONENT_INDEX][VARIABLE_INDEX]; + } + ++coarseRow; + } + coarseLevelCommunication_.reset(communication_, [](Communication*){}); + + + this->lhs_.resize(this->coarseLevelMatrix_->M()); + this->rhs_.resize(this->coarseLevelMatrix_->N()); + using OperatorArgs = typename Dune::Amg::ConstructionTraits::Arguments; + OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_); + this->operator_.reset(Dune::Amg::ConstructionTraits::construct(oargs)); + } + + // compleately unsafe!!!!!! + void calculateCoarseEntries(const Operator& fineOperator)//const M& fineMatrix) + { + 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]; + } + } + } + + + //template + // void calculateCoarseEntriesOld(const Operator& fineOperator)//const M& fineMatrix) + // { + // const auto& fineMatrix = fineOperator.getmat(); + // *coarseLevelMatrix_ = 0; + // for(auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); + // row != rowEnd; ++row) + // { + // const auto& i = (*aggregatesMap_)[row.index()]; + // if(i != AggregatesMap::ISOLATED) + // { + // for(auto entry = row->begin(), entryEnd = row->end(); + // entry != entryEnd; ++entry) + // { + // const auto& j = (*aggregatesMap_)[entry.index()]; + // if ( j != AggregatesMap::ISOLATED ) + // { + // (*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][COMPONENT_INDEX]; + // } + // } + // } + // } + // } + + void moveToCoarseLevel(const typename FatherType::FineRangeType& fine) + { + // Set coarse vector to zero + this->rhs_=0; + + auto end = fine.end(), begin=fine.begin(); + + for(auto block=begin; block != end; ++block) + { + this->rhs_[block-begin] = (*block)[COMPONENT_INDEX]; + } + + + this->lhs_=0; + } + + void moveToFineLevel(typename FatherType::FineDomainType& fine) + { + + auto end=fine.end(), begin=fine.begin(); + + for(auto block=begin; block != end; ++block) + { + (*block)[COMPONENT_INDEX] = this->lhs_[block-begin]; + } + + } + + OneComponentAggregationLevelTransferPolicyCpr* clone() const + { + return new OneComponentAggregationLevelTransferPolicyCpr(*this); + } + + const Communication& getCoarseLevelCommunication() const + { + return *coarseLevelCommunication_; + } + private: + typename Operator::matrix_type::field_type prolongDamp_; + //std::shared_ptr aggregatesMap_; + Criterion criterion_; + Communication* communication_; + std::shared_ptr coarseLevelCommunication_; + std::shared_ptr coarseLevelMatrix_; + }; + + namespace Detail + { + /** + * @brief A policy class for solving the coarse level system using one step of AMG. + * @tparam O The type of the linear operator used. + * @tparam S The type of the smoother used in AMG. + * @tparam C The type of the crition used for the aggregation within AMG. + * @tparam C1 The type of the information about the communication. Either + * Dune::OwnerOverlapCopyCommunication or Dune::SequentialInformation. + */ + template + class OneStepAMGCoarseSolverPolicyNoSolve + { + public: + typedef P LevelTransferPolicy; + /** @brief The type of the linear operator used. */ + typedef O Operator; + /** @brief The type of the range and domain of the operator. */ + typedef typename O::range_type X; + /** @brief The type of the crition used for the aggregation within AMG.*/ + typedef C Criterion; + /** @brief The type of the communication used for AMG.*/ + typedef typename P::ParallelInformation Communication; + /** @brief The type of the smoother used in AMG. */ + typedef S Smoother; + /** @brief The type of the arguments used for constructing the smoother. */ + typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; + /** @brief The type of the AMG construct on the coarse level.*/ + typedef Dune::Amg::AMGCPR AMGType; + /** + * @brief Constructs the coarse solver policy. + * @param args The arguments used for constructing the smoother. + * @param c The crition used for the aggregation within AMG. + */ + OneStepAMGCoarseSolverPolicyNoSolve(const CPRParameter* param, const SmootherArgs& args, const Criterion& c) + : param_(param), smootherArgs_(args), criterion_(c) + {} + /** @brief Copy constructor. */ + OneStepAMGCoarseSolverPolicyNoSolve(const OneStepAMGCoarseSolverPolicyNoSolve& other) + : param_(other.param_), coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_), + criterion_(other.criterion_) + {} + private: + /** + * @brief A wrapper that makes an inverse operator out of AMG. + * + * The operator will use one step of AMG to approximately solve + * the coarse level system. + */ + struct AMGInverseOperator : public Dune::InverseOperator + { + AMGInverseOperator(const CPRParameter* param, + typename AMGType::Operator& op, + const Criterion& crit, + const typename AMGType::SmootherArgs& args, + const Communication& comm) + : param_(param), amg_(),crit_(crit), op_(op),args_(args), comm_(comm) + { + amg_.reset(new AMGType(op, crit,args, comm)); + } + + void updateAmgPreconditioner(typename AMGType::Operator& op){ + //op_ = op; + //amg_->recalculateHierarchy(); + amg_->updateSolver(crit_, op, comm_); + //amg_.reset(new AMGType(op, crit_,args_, comm_)); + //amg_->recalculateGalerkin(); + } +#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; + } +#endif + + + void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) + { + DUNE_UNUSED_PARAMETER(reduction); + DUNE_UNUSED_PARAMETER(res); +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + auto sp = Dune::createScalarProduct(comm_, op_.category()); +#else + using Chooser = Dune::ScalarProductChooser; + auto sp = Chooser::construct(comm_); +#endif + Dune::Preconditioner* prec = amg_.get(); + // Linear solver parameters + const double tolerance = param_->cpr_solver_tol_; + const int maxit = param_->cpr_max_ell_iter_; + const int verbosity = ( param_->cpr_solver_verbose_ && + comm_.communicator().rank()==0 ) ? 1 : 0; + if ( param_->cpr_ell_solvetype_ == 0 ) + { + // Category of preconditioner will be checked at compile time. Therefore we need + // to cast to the derived class +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::BiCGSTABSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); +#else + Dune::BiCGSTABSolver solver(const_cast(op_), *sp, + reinterpret_cast(*prec), + tolerance, maxit, verbosity); +#endif + solver.apply(x,b,res); + + } + else if (param_->cpr_ell_solvetype_ == 1) + { + // Category of preconditioner will be checked at compile time. Therefore we need + // to cast to the derived class +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::CGSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); +#else + Dune::CGSolver solver(const_cast(op_), *sp, + reinterpret_cast(*prec), + tolerance, maxit, verbosity); +#endif + solver.apply(x,b,res); + } + else + { + // X v(x); + // prec->pre(x,b); + // op_->applyscaleadd(-1,x,b); + // v = 0; + // prec->apply(v,b); + // x += v; + // op_->applyscaleadd(-1,x,b); + // prec->post(x); +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::LoopSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); +#else + Dune::LoopSolver solver(const_cast(op_), *sp, + reinterpret_cast(*prec), + tolerance, maxit, verbosity); +#endif + solver.apply(x,b,res); + } +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + +#else + delete sp; +#endif + } + + void apply(X& x, X& b, Dune::InverseOperatorResult& res) + { + return apply(x,b,1e-8,res); + } + + ~AMGInverseOperator() + {} + AMGInverseOperator(const AMGInverseOperator& other) + : x_(other.x_), amg_(other.amg_) + { + } + private: + const CPRParameter* param_; + X x_; + std::unique_ptr amg_; + //std::unique_ptr op_; + typename AMGType::Operator& op_; + Criterion crit_; + typename AMGType::SmootherArgs args_; + const Communication& comm_; + }; + + public: + /** @brief The type of solver constructed for the coarse level. */ + typedef AMGInverseOperator CoarseLevelSolver; + + /** + * @brief Constructs a coarse level solver. + * + * @param transferPolicy The policy describing the transfer between levels. + * @return A pointer to the constructed coarse level solver. + */ + template + void setCoarseOperator(LTP& transferPolicy){ + coarseOperator_= transferPolicy.getCoarseLevelOperator(); + } + template + CoarseLevelSolver* createCoarseLevelSolver(LTP& transferPolicy) + { + coarseOperator_=transferPolicy.getCoarseLevelOperator(); + const LevelTransferPolicy& transfer = + reinterpret_cast(transferPolicy); + AMGInverseOperator* inv = new AMGInverseOperator(param_, + *coarseOperator_, + criterion_, + smootherArgs_, + transfer.getCoarseLevelCommunication()); + + return inv; //std::shared_ptr >(inv); + + } + //void recalculateGalerkin(){ + // coarseOperator_.recalculateHierarchy(); + //} + private: + /** @brief The coarse level operator. */ + std::shared_ptr coarseOperator_; + /** @brief The parameters for the CPR preconditioner. */ + const CPRParameter* param_; + /** @brief The arguments used to construct the smoother. */ + SmootherArgs smootherArgs_; + /** @brief The coarsening criterion. */ + Criterion criterion_; + }; + + + } // end namespace Detail + + + /** + * \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 first + * coarse level system, either by simply extracting the coupling between the components at COMPONENT_INDEX + * in the matrix blocks or by extracting them and 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 the 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 CoarseSmoother = typename Detail::ScalarType::value; + using FineCriterion = + typename Detail::OneComponentCriterionType::value; + using CoarseCriterion = typename Detail::ScalarType::value; + using LevelTransferPolicy = + OneComponentAggregationLevelTransferPolicyCpr; + using CoarseSolverPolicy = + Detail::OneStepAMGCoarseSolverPolicyNoSolve; + 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::scaleMatrixDRSPtr(fineOperator, comm, + COMPONENT_INDEX, weights_, param)), + scaledMatrixOperator_(Detail::createOperatorPtr(fineOperator, *scaledMatrix_, comm)), + smoother_(Detail::constructSmoother(*scaledMatrixOperator_, + smargs, comm)), + levelTransferPolicy_(criterion, comm), + coarseSolverPolicy_(¶m, smargs, criterion), + twoLevelMethod_(*scaledMatrixOperator_, + smoother_, + levelTransferPolicy_, + coarseSolverPolicy_, 0, 1) + { + } + void updatePreconditioner(const typename TwoLevelMethod::FineDomainType& weights, + const Operator& fineOperator, + const SmootherArgs& smargs, + const Communication& comm){ + weights_ = weights; + *scaledMatrix_ = *Detail::scaleMatrixDRSPtr(fineOperator, comm, + COMPONENT_INDEX, weights_, param_); + //*scaledMatrixOperator_ = *Detail::createOperatorPtr(fineOperator,*scaledMatrix_,comm); + smoother_ .reset(Detail::constructSmoother(*scaledMatrixOperator_, + smargs, comm)); + twoLevelMethod_.updatePreconditioner(*scaledMatrixOperator_, + smoother_, + coarseSolverPolicy_); + } + + void pre(typename TwoLevelMethod::FineDomainType& x, + typename TwoLevelMethod::FineRangeType& b) + { + twoLevelMethod_.pre(x,b); + } + + void post(typename TwoLevelMethod::FineDomainType& x) + { + twoLevelMethod_.post(x); + } + + void apply(typename TwoLevelMethod::FineDomainType& v, + const typename TwoLevelMethod::FineRangeType& d) + { + auto scaledD = d; + Detail::scaleVectorDRS(scaledD, COMPONENT_INDEX, param_, weights_); + twoLevelMethod_.apply(v, scaledD); + } + private: + const CPRParameter& param_; + //const typename TwoLevelMethod::FineDomainType& weights_; + typename TwoLevelMethod::FineDomainType weights_;//make copy + std::unique_ptr scaledMatrix_; + std::unique_ptr scaledMatrixOperator_; + //Operator scaledMatrixOperator_; + //std::tuple, Operator> + std::shared_ptr smoother_; + LevelTransferPolicy levelTransferPolicy_; + CoarseSolverPolicy coarseSolverPolicy_; + TwoLevelMethod twoLevelMethod_; + //BlockVector weights_; + }; + +} // end namespace Opm +#endif diff --git a/opm/autodiff/ISTLSolverEbosCpr.hpp b/opm/autodiff/ISTLSolverEbosCpr.hpp new file mode 100644 index 000000000..c3d24de56 --- /dev/null +++ b/opm/autodiff/ISTLSolverEbosCpr.hpp @@ -0,0 +1,350 @@ +/* + 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 +BEGIN_PROPERTIES +NEW_PROP_TAG(CprSmootherFine); +NEW_PROP_TAG(CprSmootherCoarse); +END_PROPERTIES + +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 + { + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; + typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, EclWellModel) WellModel; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename SparseMatrixAdapter::IstlMatrix Matrix; + typedef typename GET_PROP_TYPE(TypeTag, CprSmootherFine) CprSmootherFine; + typedef typename GET_PROP_TYPE(TypeTag, CprSmootherCoarse) CprSmootherCoarse; + //typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType; + //typedef typename Vector::block_type BlockVector; + //typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + //typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager; + //typedef typename GridView::template Codim<0>::Entity Element; + //typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + enum { pressureEqnIndex = BlackOilDefaultIndexTraits::waterCompIdx }; + enum { pressureVarIndex = Indices::pressureSwitchIdx }; + + static const int numEq = Indices::numEq; + typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false> OperatorSerial; + typedef ISTLSolverEbos SuperClass; + typedef Dune::Amg::SequentialInformation POrComm; + //typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType; + typedef Dune::MatrixAdapter MatrixAdapter; + +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + +#else + static constexpr int category = Dune::SolverCategory::sequential; + typedef Dune::ScalarProductChooser ScalarProductChooser; +#endif +//Operator MatrixOperator = Dune::MatrixAdapter + //typedef Opm::ParallelOverlappingILU0 Smoother; + typedef CprSmootherFine Smoother; + //ParallelInformation = Dune::Amg::SequentialInformation + //typedef Dune::Amg::AMG DuneAmg; + using CouplingMetric = Opm::Amg::Element; + using CritBase = Dune::Amg::SymmetricCriterion; + using Criterion = Dune::Amg::CoarsenCriterion; + typedef BlackoilAmgCpr BLACKOILAMG; + + + public: + typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType; + + 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. + ISTLSolverEbosCpr(const Simulator& simulator) + : SuperClass(simulator) + { + } + + // nothing to clean here + void eraseMatrix() { + this->matrix_for_preconditioner_.reset(); + } + + void prepare(const SparseMatrixAdapter& M, Vector& b){ + int newton_iteration = this->simulator_.model().newtonMethod().numIterations(); + // double dt = this->simulator_.timeStepSize(); + 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() ) + { + // parallel implemantation si as before + // typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true ,TypeTag> Operator; + + // auto ebosJacIgnoreOverlap = Matrix(*(this->matrix_)); + // //remove ghost rows in local matrix + // this->makeOverlapRowsInvalid(ebosJacIgnoreOverlap); + + // //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables + // //to be certain that correct matrix is used for preconditioning. + // Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel, + // this->parallelInformation_ ); + // assert( opA.comm() ); + // //SuperClass::solve( opA, x, *(this->rhs_), *(opA.comm()) ); + // typedef Dune::OwnerOverlapCopyCommunication& comm = *(opA.comm()); + // const size_t size = opA.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); + // // Construct operator, scalar product and vectors needed. + // Dune::InverseOperatorResult result; + // SuperClass::constructPreconditionerAndSolve(opA, x, *(this->rhs_), comm, result); + // SuperClass::checkConvergence(result); + + } + else +#endif + { + const WellModel& wellModel = this->simulator_.problem().wellModel(); + //typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false ,TypeTag> OperatorSerial; + opASerial_.reset(new OperatorSerial(*(this->matrix_), *(this->matrix_), wellModel)); + + //Dune::Amg::SequentialInformation info; + typedef Dune::Amg::SequentialInformation POrComm; + POrComm parallelInformation_arg; + typedef OperatorSerial LinearOperator; + + //SuperClass::constructPreconditionerAndSolve(opA, x, *(this->rhs_), info, result); + + +#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 + Vector& istlb = *(this->rhs_); + parallelInformation_arg.copyOwnerToAll(istlb, istlb); + + + + if( ! std::is_same< LinearOperator, MatrixAdapter > :: value ) + { + // create new operator in case linear operator and matrix operator differ + opA_.reset( new MatrixAdapter( opASerial_->getmat()));//, parallelInformation_arg ) ); + } + + const double relax = this->parameters_.ilu_relaxation_; + const MILU_VARIANT ilu_milu = this->parameters_.ilu_milu_; + using Matrix = typename MatrixAdapter::matrix_type; + //using CouplingMetric = Dune::Amg::Diagonal; + //using CritBase = Dune::Amg::SymmetricCriterion; + //using Criterion = Dune::Amg::CoarsenCriterion; + //using AMG = typename ISTLUtility + // ::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureIndex >::AMG; + + //std::unique_ptr< AMG > amg; + // Construct preconditioner. + //Criterion crit(15, 2000); + //SuperClass::constructAMGPrecond< Criterion >( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu ); + // ISTLUtility::template createAMGPreconditionerPointer( *opA_, + // relax, + // parallelInformation_arg, + // amg_, + // this->parameters_, + // this->weights_ ); + //using AMG = BlackoilAmg; + POrComm& comm = parallelInformation_arg; + 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 + typedef typename BLACKOILAMG::Smoother Smoother; + typedef typename BLACKOILAMG::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); + //ISTLUtility::setILUParameters(smootherArgs, params); + //smootherArgs.setN(params.cpr_ilu_n_); smootherArgs.setMilu(params.cpr_ilu_milu_); + + MatrixAdapter& opARef = *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 BLACKOILAMG( params, this->weights_, opARef, criterion, smootherArgs, comm ) ); + }else{ + if(this->parameters_.cpr_solver_verbose_){ + std::cout << " Only update amg solver " << std::endl; + } + amg_->updatePreconditioner(this->weights_,opARef, smootherArgs, comm); + } + // Solve. + //SuperClass::solve(linearOperator, x, istlb, *sp, *amg, result); + //references seems to do something els than refering + + int verbosity_linsolve = ( this->isIORank_ ) ? this->parameters_.linear_solver_verbosity_ : 0; + LinearOperator& opASerialRef = *opASerial_; + linsolve_.reset(new Dune::BiCGSTABSolver(opASerialRef, *sp_, *amg_, + this->parameters_.linear_solver_reduction_, + this->parameters_.linear_solver_maxiter_, + verbosity_linsolve)); + + + } + } + + bool solve(Vector& x) { + //SuperClass::solve(x); + if( this->isParallel() ){ + // for now only call the superclass + bool converged = SuperClass::solve(x); + return converged; + }else{ + // 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_; + } + + + /// Solve the system of linear equations Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] residual residual object containing A and b. + /// \return the solution x + + /// \copydoc NewtonIterationBlackoilInterface::iterations + int iterations () const { return this->iterations_; } + + /// \copydoc NewtonIterationBlackoilInterface::parallelInformation + const boost::any& parallelInformation() const { return this->parallelInformation_; } + + protected: + + //using Matrix = typename MatrixAdapter::matrix_type; + //using CouplingMetric = Dune::Amg::Diagonal; + //using CritBase = Dune::Amg::SymmetricCriterion; + //using Criterion = Dune::Amg::CoarsenCriterion; + //using AMG = typename ISTLUtility + // ::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureIndex >::AMG; + //Operator MatrixOperator = Dune::MatrixAdapter + //Smoother ParallelOverLappingILU0 + //ParallelInformation = Dune::Amg::SequentialInformation +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + typedef std::shared_ptr< Dune::ScalarProduct > SPPointer; +#else + typedef std::unique_ptr SPPointer; +#endif + std::unique_ptr< MatrixAdapter > opA_; + std::unique_ptr< OperatorSerial > opASerial_; + std::unique_ptr< BLACKOILAMG > amg_; + SPPointer sp_; + std::shared_ptr< Dune::BiCGSTABSolver > linsolve_; + //std::shared_ptr< Dune::LinearOperator > op_; + //std::shared_ptr< Dune::Preconditioner > prec_; + //std::shared_ptr< Dune::ScalarProduct > sp_; + //Vector solution_; + + }; // end ISTLSolver + +} // namespace Opm +#endif diff --git a/opm/autodiff/amgcpr.hh b/opm/autodiff/amgcpr.hh new file mode 100644 index 000000000..59645ca43 --- /dev/null +++ b/opm/autodiff/amgcpr.hh @@ -0,0 +1,1027 @@ +// -*- 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 + +#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); + + //! Category of the preconditioner (see SolverCategory::Category) + virtual SolverCategory::Category category() const + { + return category_; + } + + /** \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()); + } + + template + void updateSolver(C& criterion, Operator& matrix, const PI& pinfo) + { + Timer watch; + //this->createHierarchies(criterion, const_cast(matrix), pinfo); + smoothers_.reset(new Hierarchy); + solver_.reset(); + coarseSmoother_.reset(); + scalarProduct_.reset(); + buildHierarchy_= true; + coarsesolverconverged = true; + //this->createHierarchies(criterion, matrix, pinfo); + //matrices_.reset(new OperatorHierarchy(matrix, pinfo)); + //matrices_->template build >(criterion); + 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 "<maxlevels()<<" levels " + <levels()==matrices_->maxlevels()) + { + // We have the carsest level. Create the coarse Solver + SmootherArgs sargs(smootherArgs_); + sargs.iterations = 1; + + typename ConstructionTraits::Arguments cargs; + cargs.setArgs(sargs); + + cargs.setMatrix(matrices_->matrices().coarsest()->getmat()); + cargs.setComm(*matrices_->parallelInformation().coarsest()); + + coarseSmoother_.reset(ConstructionTraits::construct(cargs)); + scalarProduct_ = createScalarProduct(cargs.getComm(),category()); + + 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 && false) + { + solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false)); + if(verbosity_>0 ) + std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl; + } + else + { + // solver_.reset(new BiCGSTABSolver(const_cast(*matrices_->matrices().coarsest()), + // *scalarProduct_, + // *coarseSmoother_, 1E-2, 1000, 0)); + solver_.reset(new LoopSolver(const_cast(*matrices_->matrices().coarsest()), + *scalarProduct_, + *coarseSmoother_, 1E-2, 1, 0)); + } + } + + // if(verbosity_>0 && + // matrices_->parallelInformation().finest()->communicator().rank()==0){ + // std::cout<<"Recalculating galerkin and coarse somothers "<maxlevels()<<" levels " + // <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)); + scalarProduct_ = createScalarProduct(cargs.getComm(),category()); + + 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()), + *scalarProduct_, + *coarseSmoother_, 1E-2, 1000, 0)); + else + solver_.reset(); + }else + { + solver_.reset(new BiCGSTABSolver(const_cast(*matrices_->matrices().coarsest()), + *scalarProduct_, + *coarseSmoother_, 1E-2, 1000, 0)); + + } + } + } + + // if(verbosity_>0 && + // matrices_->parallelInformation().finest()->communicator().rank()==0){ + // std::cout<<"Recalculating galerkin and coarse somothers "<maxlevels()<<" levels " + // < + void createHierarchies(C& criterion, Operator& matrix, + const PI& pinfo); + /** + * @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_; + /** @brief The solver category. */ + SolverCategory::Category category_; + /** @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_), + category_(amg.category_), + 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_(), +// #warning should category be retrieved from matrices? + category_(SolverCategory::category(*smoothers_->coarsest())), + 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_(), + category_(SolverCategory::category(pinfo)), + verbosity_(criterion.debugLevel()) + { + if(SolverCategory::category(matrix) != SolverCategory::category(pinfo)) + DUNE_THROW(InvalidSolverCategory, "Matrix and Communication must have the same SolverCategory!"); + // TODO: reestablish compile time checks. + //static_assert(static_cast(PI::category)==static_cast(S::category), + // "Matrix and Solver must match in terms of category!"); + 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::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 "<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)); + // scalarProduct_ = createScalarProduct(cargs.getComm(),category()); + + // 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()), + // *scalarProduct_, + // *coarseSmoother_, 1E-2, 1000, 0)); + // else + // solver_.reset(); + // }else + // { + // solver_.reset(new BiCGSTABSolver(const_cast(*matrices_->matrices().coarsest()), + // *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); + // // }); + // } + // } + // } + + // if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) + // std::cout<<"Building hierarchy of "<maxlevels()<<" levels " + // <<"(inclusive coarse solver) took "< + 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..fe608b4ed --- /dev/null +++ b/opm/autodiff/twolevelmethodcpr.hh @@ -0,0 +1,567 @@ +// -*- 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 + +#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; + + //template + virtual void calculateCoarseEntries(const FineOperatorType& fineOperator) = 0; + //virtual void recalculateGalerkin(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="<createCoarseLevelSystem(*operator_); + policy_->calculateCoarseEntries(*operator_); + //policy_->calculateCoarseEntries(7); + coarsePolicy.setCoarseOperator(*policy_); + //delete coarseSolver_; + //coarseSolver_ = coarsePolicy.createCoarseLevelSolver(*policy_); + coarseSolver_->updateAmgPreconditioner(*(policy_->getCoarseLevelOperator())); + }else{ + // we should probably not be heere + policy_->createCoarseLevelSystem(*operator_); + coarseSolver_ = coarsePolicy.createCoarseLevelSolver(*policy_); + } + } + + + void pre(FineDomainType& x, FineRangeType& b) + { + smoother_->pre(x,b); + } + + void post(FineDomainType& x) + { + DUNE_UNUSED_PARAMETER(x); + } + + void apply(FineDomainType& v, const FineRangeType& d) + { + FineDomainType u(v); + FineRangeType rhs(d); + LevelContext context; + SequentialInformation info; + context.pinfo=&info; + context.lhs=&u; + context.update=&v; + context.smoother=smoother_; + context.rhs=&rhs; + context.matrix=operator_; + // Presmoothing + presmooth(context, preSteps_); + //Coarse grid correction + policy_->moveToCoarseLevel(*context.rhs); + InverseOperatorResult res; + coarseSolver_->apply(policy_->getCoarseLevelLhs(), policy_->getCoarseLevelRhs(), res); + *context.lhs=0; + policy_->moveToFineLevel(*context.lhs); + *context.update += *context.lhs; + // Postsmoothing + postsmooth(context, postSteps_); + + } +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) +// //! Category of the preconditioner (see SolverCategory::Category) + virtual SolverCategory::Category category() const + { + return SolverCategory::sequential; + } +#endif +private: + /** + * @brief Struct containing the level information. + */ + struct LevelContext + { + /** @brief The type of the smoother used. */ + typedef S SmootherType; + /** @brief A pointer to the smoother. */ + std::shared_ptr smoother; + /** @brief The left hand side passed to the and returned by the smoother. */ + FineDomainType* lhs; + /* + * @brief The right hand side holding the current residual. + * + * This is passed to the smoother as the right hand side. + */ + FineRangeType* rhs; + /** + * @brief The total update calculated by the preconditioner. + * + * I.e. all update from smoothing and coarse grid correction summed up. + */ + FineDomainType* update; + /** @parallel information */ + SequentialInformation* pinfo; + /** + * @brief The matrix that we are solving. + * + * Needed to update the residual. + */ + const FineOperatorType* matrix; + }; + FineOperatorType* operator_; + /** @brief The coarse level solver. */ + CoarseLevelSolver* coarseSolver_; + /** @brief The fine level smoother. */ + std::shared_ptr smoother_; + /** @brief Policy for prolongation, restriction, and coarse level system creation. */ + LevelTransferPolicyCpr* policy_; + /** @brief The number of presmoothing steps to apply. */ + std::size_t preSteps_; + /** @brief The number of postsmoothing steps to apply. */ + std::size_t postSteps_; +}; +}// end namespace Amg +}// end namespace Dune + +/** @} */ +#endif From ed24b3fcadfd3c7df4229749143590b668efad61 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 21 Mar 2019 22:32:52 +0100 Subject: [PATCH 02/25] changed code to use direct solver on coarses scale for amg. This should runtime stettings for this.. --- opm/autodiff/amgcpr.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/amgcpr.hh b/opm/autodiff/amgcpr.hh index 59645ca43..1eb1c72ac 100644 --- a/opm/autodiff/amgcpr.hh +++ b/opm/autodiff/amgcpr.hh @@ -213,7 +213,7 @@ namespace Dune 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 && false) + if( SolverSelector::isDirectSolver) { solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false)); if(verbosity_>0 ) From 4f1603407d024088d82db182454b50995f821850 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 4 Apr 2019 12:20:30 +0200 Subject: [PATCH 03/25] Added timing of prepare call for linear solvers --- opm/autodiff/BlackoilModelEbos.hpp | 9 +++++++-- opm/core/simulator/SimulatorReport.cpp | 10 ++++++++++ opm/core/simulator/SimulatorReport.hpp | 1 + 3 files changed, 18 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 92f562672..570bc652b 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -297,13 +297,15 @@ namespace Opm { // equations wellModel().linearize(ebosSimulator().model().linearizer().jacobian(), ebosSimulator().model().linearizer().residual()); - + 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 +475,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 +911,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/core/simulator/SimulatorReport.cpp b/opm/core/simulator/SimulatorReport.cpp index 83854ac4a..83878f7a3 100644 --- a/opm/core/simulator/SimulatorReport.cpp +++ b/opm/core/simulator/SimulatorReport.cpp @@ -33,6 +33,7 @@ namespace Opm total_time(0.0), solver_time(0.0), assemble_time(0.0), + linear_solve_setup_time(0.0), linear_solve_time(0.0), update_time(0.0), output_write_time(0.0), @@ -49,6 +50,7 @@ namespace Opm { pressure_time += sr.pressure_time; transport_time += sr.transport_time; + linear_solve_setup_time += sr.linear_solve_setup_time; linear_solve_time += sr.linear_solve_time; solver_time += sr.solver_time; assemble_time += sr.assemble_time; @@ -119,6 +121,14 @@ namespace Opm } os << std::endl; + t = linear_solve_setup_time + (failureReport ? failureReport->linear_solve_setup_time : 0.0); + os << " Linear solve setup time (seconds): " << t; + if (failureReport) { + os << " (Failed: " << failureReport->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; From 28dfef1ecf3ef8cf16ce48be1be91701a775d4ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 4 Apr 2019 15:52:39 +0200 Subject: [PATCH 04/25] Compile fix: do not use undefined member. --- opm/autodiff/ISTLSolverEbosCpr.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/ISTLSolverEbosCpr.hpp b/opm/autodiff/ISTLSolverEbosCpr.hpp index c3d24de56..0939788aa 100644 --- a/opm/autodiff/ISTLSolverEbosCpr.hpp +++ b/opm/autodiff/ISTLSolverEbosCpr.hpp @@ -274,7 +274,11 @@ namespace Opm //SuperClass::solve(linearOperator, x, istlb, *sp, *amg, result); //references seems to do something els than refering - int verbosity_linsolve = ( this->isIORank_ ) ? this->parameters_.linear_solver_verbosity_ : 0; + int verbosity_linsolve = 0; + if (comm.communicator().rank() == 0) { + verbosity_linsolve = this->parameters_.linear_solver_verbosity_; + } + LinearOperator& opASerialRef = *opASerial_; linsolve_.reset(new Dune::BiCGSTABSolver(opASerialRef, *sp_, *amg_, this->parameters_.linear_solver_reduction_, From b39815edb565276658dace556efe5d1f34b7bdd9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 4 Apr 2019 16:05:41 +0200 Subject: [PATCH 05/25] Minor formatting and removing unused code. --- opm/autodiff/amgcpr.hh | 141 ++++++----------------------------------- 1 file changed, 18 insertions(+), 123 deletions(-) diff --git a/opm/autodiff/amgcpr.hh b/opm/autodiff/amgcpr.hh index 1eb1c72ac..75157db72 100644 --- a/opm/autodiff/amgcpr.hh +++ b/opm/autodiff/amgcpr.hh @@ -167,34 +167,36 @@ namespace Dune matrices_->recalculateGalerkin(NegateSet()); } + /** + * @brief Update the coarse solver and the hierarchies. + */ template - void updateSolver(C& criterion, Operator& matrix, const PI& pinfo) + void updateSolver(C& criterion, Operator& /* matrix */, const PI& pinfo) { Timer watch; - //this->createHierarchies(criterion, const_cast(matrix), pinfo); smoothers_.reset(new Hierarchy); solver_.reset(); coarseSmoother_.reset(); scalarProduct_.reset(); buildHierarchy_= true; coarsesolverconverged = true; - //this->createHierarchies(criterion, matrix, pinfo); - //matrices_.reset(new OperatorHierarchy(matrix, pinfo)); - //matrices_->template build >(criterion); 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 "<maxlevels()<<" levels " - <0 && matrices_->parallelInformation().finest()->communicator().rank()==0) { + std::cout << "Recalculating galerkin and coarse somothers "<< matrices_->maxlevels() << " levels " + << watch.elapsed() << " seconds." << std::endl; } } - void setupCoarseSolver(){ + + void setupCoarseSolver() + { setupCoarseSolverSerial(); } - void setupCoarseSolverSerial(){ + + void setupCoarseSolverSerial() + { if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) { // We have the carsest level. Create the coarse Solver @@ -229,16 +231,10 @@ namespace Dune *coarseSmoother_, 1E-2, 1, 0)); } } - - // if(verbosity_>0 && - // matrices_->parallelInformation().finest()->communicator().rank()==0){ - // std::cout<<"Recalculating galerkin and coarse somothers "<maxlevels()<<" levels " - // <levels()==matrices_->maxlevels() && ( ! matrices_->redistributeInformation().back().isSetup() || matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) ) @@ -312,13 +308,8 @@ namespace Dune } } } - - // if(verbosity_>0 && - // matrices_->parallelInformation().finest()->communicator().rank()==0){ - // std::cout<<"Recalculating galerkin and coarse somothers "<maxlevels()<<" levels " - // <coarsenSmoother(*smoothers_, smootherArgs_); setupCoarseSolver(); if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) - std::cout<<"Building hierarchy of "<maxlevels()<<" levels " - <<"(inclusive coarse solver) took "<maxlevels()<<" levels " + <<"(inclusive coarse solver) took "<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)); - // scalarProduct_ = createScalarProduct(cargs.getComm(),category()); - - // 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()), - // *scalarProduct_, - // *coarseSmoother_, 1E-2, 1000, 0)); - // else - // solver_.reset(); - // }else - // { - // solver_.reset(new BiCGSTABSolver(const_cast(*matrices_->matrices().coarsest()), - // *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); - // // }); - // } - // } - // } - - // if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) - // std::cout<<"Building hierarchy of "<maxlevels()<<" levels " - // <<"(inclusive coarse solver) took "< From 13c6463f89a4466d6e91aec94f899b5266b309c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 4 Apr 2019 16:45:23 +0200 Subject: [PATCH 06/25] Refactor to be closer to the original amg.hh. Should minimize the diff when comparing. --- opm/autodiff/amgcpr.hh | 257 +++++++++++++++++++---------------------- 1 file changed, 119 insertions(+), 138 deletions(-) diff --git a/opm/autodiff/amgcpr.hh b/opm/autodiff/amgcpr.hh index 75157db72..696756f02 100644 --- a/opm/autodiff/amgcpr.hh +++ b/opm/autodiff/amgcpr.hh @@ -171,144 +171,7 @@ namespace Dune * @brief Update the coarse solver and the hierarchies. */ template - void 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; - } - } - - void setupCoarseSolver() - { - setupCoarseSolverSerial(); - } - - void setupCoarseSolverSerial() - { - if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) - { - // We have the carsest level. Create the coarse Solver - SmootherArgs sargs(smootherArgs_); - sargs.iterations = 1; - - typename ConstructionTraits::Arguments cargs; - cargs.setArgs(sargs); - - cargs.setMatrix(matrices_->matrices().coarsest()->getmat()); - cargs.setComm(*matrices_->parallelInformation().coarsest()); - - coarseSmoother_.reset(ConstructionTraits::construct(cargs)); - scalarProduct_ = createScalarProduct(cargs.getComm(),category()); - - 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) - { - solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false)); - if(verbosity_>0 ) - std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl; - } - else - { - // solver_.reset(new BiCGSTABSolver(const_cast(*matrices_->matrices().coarsest()), - // *scalarProduct_, - // *coarseSmoother_, 1E-2, 1000, 0)); - solver_.reset(new LoopSolver(const_cast(*matrices_->matrices().coarsest()), - *scalarProduct_, - *coarseSmoother_, 1E-2, 1, 0)); - } - } - } - - void setupCoarseSolverPar() - { - 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)); - scalarProduct_ = createScalarProduct(cargs.getComm(),category()); - - 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()), - *scalarProduct_, - *coarseSmoother_, 1E-2, 1000, 0)); - else - solver_.reset(); - }else - { - solver_.reset(new BiCGSTABSolver(const_cast(*matrices_->matrices().coarsest()), - *scalarProduct_, - *coarseSmoother_, 1E-2, 1000, 0)); - - } - } - } - } + void updateSolver(C& criterion, Operator& /* matrix */, const PI& pinfo); /** * @brief Check whether the coarse solver used is a direct solver. @@ -326,6 +189,9 @@ namespace Dune template void createHierarchies(C& criterion, Operator& matrix, const PI& pinfo); + + void setupCoarseSolver(); + /** * @brief A struct that holds the context of the current level. * @@ -523,6 +389,27 @@ namespace Dune 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, @@ -541,6 +428,100 @@ namespace Dune <<"(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)); + scalarProduct_ = createScalarProduct(cargs.getComm(),category()); + + 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()), + *scalarProduct_, + *coarseSmoother_, 1E-2, 1000, 0)); + else + solver_.reset(); + }else + { + solver_.reset(new BiCGSTABSolver(const_cast(*matrices_->matrices().coarsest()), + *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) From 29837be2e306bc9be49b8a807245bc40b58c8d59 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 Apr 2019 10:59:28 +0200 Subject: [PATCH 07/25] Changes to minimize differences to original istl versions. Also added note pointing to original istl files. --- opm/autodiff/amgcpr.hh | 3 ++ opm/autodiff/twolevelmethodcpr.hh | 48 ++++++++++++++----------------- 2 files changed, 24 insertions(+), 27 deletions(-) diff --git a/opm/autodiff/amgcpr.hh b/opm/autodiff/amgcpr.hh index 696756f02..0d871e782 100644 --- a/opm/autodiff/amgcpr.hh +++ b/opm/autodiff/amgcpr.hh @@ -3,6 +3,9 @@ #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 diff --git a/opm/autodiff/twolevelmethodcpr.hh b/opm/autodiff/twolevelmethodcpr.hh index fe608b4ed..f6415d24a 100644 --- a/opm/autodiff/twolevelmethodcpr.hh +++ b/opm/autodiff/twolevelmethodcpr.hh @@ -3,6 +3,9 @@ #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 @@ -119,9 +122,10 @@ public: */ virtual void createCoarseLevelSystem(const FineOperatorType& fineOperator)=0; - //template + /** + * @brief ???. + */ virtual void calculateCoarseEntries(const FineOperatorType& fineOperator) = 0; - //virtual void recalculateGalerkin(FineOperatorType& fineOperator)=0; /** @brief Clone the current object. */ virtual LevelTransferPolicyCpr* clone() const =0; @@ -202,7 +206,7 @@ public: this->rhs_.resize(this->matrix_->N()); this->operator_.reset(new O(*matrix_)); } - + void moveToCoarseLevel(const typename FatherType::FineRangeType& fineRhs) { Transfer @@ -264,7 +268,7 @@ public: : coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_), criterion_(other.criterion_) {} - private: +private: /** * @brief A wrapper that makes an inverse operator out of AMG. * @@ -341,10 +345,7 @@ public: return inv; //std::shared_ptr >(inv); } - //void recalculateGalerkin(){//coarseOperator){ - //coarseOperator_.resetOperator - //coarseOperator_.recalculateHierarchy(); - //} + private: /** @brief The coarse level operator. */ std::shared_ptr coarseOperator_; @@ -423,10 +424,10 @@ public: * @param preSteps The number of smoothing steps to apply after the coarse * level correction. */ - TwoLevelMethodCpr(FineOperatorType& op, + TwoLevelMethodCpr(const FineOperatorType& op, std::shared_ptr smoother, const LevelTransferPolicyCpr& policy, + CoarseOperatorType>& policy, CoarseLevelSolverPolicy& coarsePolicy, std::size_t preSteps=1, std::size_t postSteps=1) : operator_(&op), smoother_(smoother), @@ -434,7 +435,7 @@ public: { policy_ = policy.clone(); policy_->createCoarseLevelSystem(*operator_); - coarseSolver_= coarsePolicy.createCoarseLevelSolver(*policy_); + coarseSolver_=coarsePolicy.createCoarseLevelSolver(*policy_); } TwoLevelMethodCpr(const TwoLevelMethodCpr& other) @@ -453,36 +454,29 @@ public: std::shared_ptr smoother, CoarseLevelSolverPolicy& coarsePolicy){ //assume new matrix is not reallocated the new precondition should anyway be made - //operator_ = &op;// hope fine scale operator is the same - smoother_ = smoother; - if(not(coarseSolver_ == 0)){ - //delete coarseSolver_; - //std::cout << " Only rebuild hirarchy " << std::endl; - //policy_->createCoarseLevelSystem(*operator_); + smoother_ = smoother; + if (coarseSolver_) { policy_->calculateCoarseEntries(*operator_); - //policy_->calculateCoarseEntries(7); coarsePolicy.setCoarseOperator(*policy_); - //delete coarseSolver_; - //coarseSolver_ = coarsePolicy.createCoarseLevelSolver(*policy_); coarseSolver_->updateAmgPreconditioner(*(policy_->getCoarseLevelOperator())); - }else{ + } else { // we should probably not be heere policy_->createCoarseLevelSystem(*operator_); - coarseSolver_ = coarsePolicy.createCoarseLevelSolver(*policy_); - } + coarseSolver_ = coarsePolicy.createCoarseLevelSolver(*policy_); + } } - + void pre(FineDomainType& x, FineRangeType& b) { smoother_->pre(x,b); } - + void post(FineDomainType& x) { DUNE_UNUSED_PARAMETER(x); } - + void apply(FineDomainType& v, const FineRangeType& d) { FineDomainType u(v); @@ -548,7 +542,7 @@ private: */ const FineOperatorType* matrix; }; - FineOperatorType* operator_; + const FineOperatorType* operator_; /** @brief The coarse level solver. */ CoarseLevelSolver* coarseSolver_; /** @brief The fine level smoother. */ From 5f4f8b71a87fdafb8c566f8ee7cfced95c2168b5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 Apr 2019 11:14:13 +0200 Subject: [PATCH 08/25] Silence override and member order warnings. --- opm/autodiff/BlackoilAmg.hpp | 10 +++++----- opm/autodiff/BlackoilAmgCpr.hpp | 12 ++++++------ opm/autodiff/ISTLSolverEbos.hpp | 6 +++--- opm/autodiff/ParallelOverlappingILU0.hpp | 6 +++--- 4 files changed, 17 insertions(+), 17 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index 5f60c8528..638ff0d49 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -408,7 +408,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 +518,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); } @@ -1022,18 +1022,18 @@ public: {} 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_); diff --git a/opm/autodiff/BlackoilAmgCpr.hpp b/opm/autodiff/BlackoilAmgCpr.hpp index 55436c00e..e35c53589 100644 --- a/opm/autodiff/BlackoilAmgCpr.hpp +++ b/opm/autodiff/BlackoilAmgCpr.hpp @@ -287,7 +287,7 @@ namespace Opm const Criterion& crit, const typename AMGType::SmootherArgs& args, const Communication& comm) - : param_(param), amg_(),crit_(crit), op_(op),args_(args), comm_(comm) + : param_(param), amg_(), op_(op), crit_(crit), args_(args), comm_(comm) { amg_.reset(new AMGType(op, crit,args, comm)); } @@ -308,7 +308,7 @@ namespace Opm #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); @@ -380,7 +380,7 @@ namespace Opm #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); } @@ -561,18 +561,18 @@ namespace Opm } 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_); diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index 1b56b09f2..0013e1578 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -130,7 +130,7 @@ public: #endif } - 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 +144,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,7 +157,7 @@ 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() { 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); } From c92d84432ccbb6a4296f3f4e04b7942330a98dda Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 Apr 2019 11:21:36 +0200 Subject: [PATCH 09/25] Replace tabs with spaces. --- opm/autodiff/BlackoilModelEbos.hpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 570bc652b..64aaa51b8 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -297,15 +297,17 @@ namespace Opm { // equations wellModel().linearize(ebosSimulator().model().linearizer().jacobian(), ebosSimulator().model().linearizer().residual()); - linear_solve_setup_time_ = 0.0; + + // 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_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_setup_time += linear_solve_setup_time_; report.linear_solve_time += perfTimer.stop(); report.total_linear_iterations += linearIterationsLastSolve(); @@ -475,10 +477,10 @@ namespace Opm { x = 0.0; auto& ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver(); - Dune::Timer perfTimer; + Dune::Timer perfTimer; perfTimer.start(); ebosSolver.prepare(ebosJac, ebosResid); - linear_solve_setup_time_ = perfTimer.stop(); + 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 From 5d68a308b57234d7783723eed3bf5f2a74142b18 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 Apr 2019 15:08:14 +0200 Subject: [PATCH 10/25] Fix memory bug by refactoring. Moving out of a matrix in the same function argument list as the matrix is passed to another function is wrong. Simplified by avoiding the unwieldy tuple. --- opm/autodiff/BlackoilAmg.hpp | 43 ++++++++++++++---------------------- 1 file changed, 16 insertions(+), 27 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index 638ff0d49..fc6a649de 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -77,18 +77,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 +84,23 @@ 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 +111,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 +132,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) @@ -1011,13 +999,13 @@ 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) {} @@ -1042,7 +1030,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_; From 0e1db0c6bdff4661ab0e607b029bd318d538985f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 Apr 2019 15:10:03 +0200 Subject: [PATCH 11/25] Reindent, remove unused function (after refactoring). --- opm/autodiff/BlackoilAmgCpr.hpp | 964 +++++++++++++++----------------- 1 file changed, 456 insertions(+), 508 deletions(-) diff --git a/opm/autodiff/BlackoilAmgCpr.hpp b/opm/autodiff/BlackoilAmgCpr.hpp index e35c53589..25adc1389 100644 --- a/opm/autodiff/BlackoilAmgCpr.hpp +++ b/opm/autodiff/BlackoilAmgCpr.hpp @@ -35,563 +35,511 @@ #include #include #include -namespace Dune -{ - namespace Amg - { - template - class UnSymmetricCriterion; - } -} - -namespace Dune -{ - - template - class MatrixBlock; - -} namespace Opm { - namespace Detail - { - template - std::unique_ptr scaleMatrixDRSPtr(const Operator& op, - const Communication& comm, - std::size_t pressureEqnIndex, - const Vector& weights, - const Opm::CPRParameter& param) - { - using Matrix = typename Operator::matrix_type; - using Block = typename Matrix::block_type; - using BlockVector = typename Vector::block_type; - std::unique_ptr matrix(new Matrix(op.getmat())); - if (param.cpr_use_drs_) { - const auto endi = matrix->end(); - for (auto i = matrix->begin(); i != endi; ++i) { - const BlockVector& bw = weights[i.index()]; - const auto endj = (*i).end(); - for (auto j = (*i).begin(); j != endj; ++j) { - Block& block = *j; - BlockVector& bvec = block[pressureEqnIndex]; - // should introduce limits which also change the weights - block.mtv(bw, bvec); - } - } - } - return matrix;//, createOperator(op, *matrix, comm)); - } - } - - template - class OneComponentAggregationLevelTransferPolicyCpr; - - template + + template class OneComponentAggregationLevelTransferPolicyCpr - : public Dune::Amg::LevelTransferPolicyCpr::value> + : public Dune::Amg::LevelTransferPolicyCpr::value> { - typedef Dune::Amg::AggregatesMap AggregatesMap; + typedef Dune::Amg::AggregatesMap AggregatesMap; public: - using CoarseOperator = typename Detail::ScalarType::value; - typedef Dune::Amg::LevelTransferPolicy FatherType; - typedef Communication ParallelInformation; + using CoarseOperator = typename Detail::ScalarType::value; + typedef Dune::Amg::LevelTransferPolicy FatherType; + typedef Communication ParallelInformation; public: - OneComponentAggregationLevelTransferPolicyCpr(const Criterion& crit, const Communication& comm) - : criterion_(crit), communication_(&const_cast(comm)) - {} + OneComponentAggregationLevelTransferPolicyCpr(const Criterion& crit, const Communication& comm) + : criterion_(crit), communication_(&const_cast(comm)) + {} - void createCoarseLevelSystem(const Operator& fineOperator) - { - prolongDamp_ = 1; + void createCoarseLevelSystem(const Operator& fineOperator) + { + prolongDamp_ = 1; - using CoarseMatrix = typename CoarseOperator::matrix_type; - const auto& fineLevelMatrix = fineOperator.getmat(); - coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise)); - auto createIter = coarseLevelMatrix_->createbegin(); + using CoarseMatrix = typename CoarseOperator::matrix_type; + const auto& fineLevelMatrix = fineOperator.getmat(); + coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise)); + auto createIter = coarseLevelMatrix_->createbegin(); - for ( const auto& row: fineLevelMatrix ) - { - for ( auto col = row.begin(), cend = row.end(); col != cend; ++col) - { - createIter.insert(col.index()); - } - ++createIter; - } + for ( const auto& row: fineLevelMatrix ) + { + for ( auto col = row.begin(), cend = row.end(); col != cend; ++col) + { + createIter.insert(col.index()); + } + ++createIter; + } - auto coarseRow = coarseLevelMatrix_->begin(); - for ( const auto& row: fineLevelMatrix ) - { - auto coarseCol = coarseRow->begin(); + auto coarseRow = coarseLevelMatrix_->begin(); + for ( const auto& row: fineLevelMatrix ) + { + auto coarseCol = coarseRow->begin(); - for ( auto col = row.begin(), cend = row.end(); col != cend; ++col, ++coarseCol ) - { - assert( col.index() == coarseCol.index() ); - *coarseCol = (*col)[COMPONENT_INDEX][VARIABLE_INDEX]; - } - ++coarseRow; - } - coarseLevelCommunication_.reset(communication_, [](Communication*){}); + for ( auto col = row.begin(), cend = row.end(); col != cend; ++col, ++coarseCol ) + { + assert( col.index() == coarseCol.index() ); + *coarseCol = (*col)[COMPONENT_INDEX][VARIABLE_INDEX]; + } + ++coarseRow; + } + coarseLevelCommunication_.reset(communication_, [](Communication*){}); - this->lhs_.resize(this->coarseLevelMatrix_->M()); - this->rhs_.resize(this->coarseLevelMatrix_->N()); - using OperatorArgs = typename Dune::Amg::ConstructionTraits::Arguments; - OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_); - this->operator_.reset(Dune::Amg::ConstructionTraits::construct(oargs)); - } + this->lhs_.resize(this->coarseLevelMatrix_->M()); + this->rhs_.resize(this->coarseLevelMatrix_->N()); + using OperatorArgs = typename Dune::Amg::ConstructionTraits::Arguments; + OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_); + this->operator_.reset(Dune::Amg::ConstructionTraits::construct(oargs)); + } - // compleately unsafe!!!!!! - void calculateCoarseEntries(const Operator& fineOperator)//const M& fineMatrix) - { - 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]; - } - } - } + // compleately unsafe!!!!!! + void calculateCoarseEntries(const Operator& fineOperator)//const M& fineMatrix) + { + 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]; + } + } + } - //template - // void calculateCoarseEntriesOld(const Operator& fineOperator)//const M& fineMatrix) - // { - // const auto& fineMatrix = fineOperator.getmat(); - // *coarseLevelMatrix_ = 0; - // for(auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); - // row != rowEnd; ++row) - // { - // const auto& i = (*aggregatesMap_)[row.index()]; - // if(i != AggregatesMap::ISOLATED) - // { - // for(auto entry = row->begin(), entryEnd = row->end(); - // entry != entryEnd; ++entry) - // { - // const auto& j = (*aggregatesMap_)[entry.index()]; - // if ( j != AggregatesMap::ISOLATED ) - // { - // (*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][COMPONENT_INDEX]; - // } - // } - // } - // } - // } + //template + // void calculateCoarseEntriesOld(const Operator& fineOperator)//const M& fineMatrix) + // { + // const auto& fineMatrix = fineOperator.getmat(); + // *coarseLevelMatrix_ = 0; + // for(auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); + // row != rowEnd; ++row) + // { + // const auto& i = (*aggregatesMap_)[row.index()]; + // if(i != AggregatesMap::ISOLATED) + // { + // for(auto entry = row->begin(), entryEnd = row->end(); + // entry != entryEnd; ++entry) + // { + // const auto& j = (*aggregatesMap_)[entry.index()]; + // if ( j != AggregatesMap::ISOLATED ) + // { + // (*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][COMPONENT_INDEX]; + // } + // } + // } + // } + // } - void moveToCoarseLevel(const typename FatherType::FineRangeType& fine) - { - // Set coarse vector to zero - this->rhs_=0; + void moveToCoarseLevel(const typename FatherType::FineRangeType& fine) + { + // Set coarse vector to zero + this->rhs_=0; - auto end = fine.end(), begin=fine.begin(); + auto end = fine.end(), begin=fine.begin(); - for(auto block=begin; block != end; ++block) - { - this->rhs_[block-begin] = (*block)[COMPONENT_INDEX]; - } + for(auto block=begin; block != end; ++block) + { + this->rhs_[block-begin] = (*block)[COMPONENT_INDEX]; + } - this->lhs_=0; - } + this->lhs_=0; + } - void moveToFineLevel(typename FatherType::FineDomainType& fine) - { + void moveToFineLevel(typename FatherType::FineDomainType& fine) + { - auto end=fine.end(), begin=fine.begin(); + auto end=fine.end(), begin=fine.begin(); - for(auto block=begin; block != end; ++block) - { - (*block)[COMPONENT_INDEX] = this->lhs_[block-begin]; - } + for(auto block=begin; block != end; ++block) + { + (*block)[COMPONENT_INDEX] = this->lhs_[block-begin]; + } - } + } - OneComponentAggregationLevelTransferPolicyCpr* clone() const - { - return new OneComponentAggregationLevelTransferPolicyCpr(*this); - } + OneComponentAggregationLevelTransferPolicyCpr* clone() const + { + return new OneComponentAggregationLevelTransferPolicyCpr(*this); + } - const Communication& getCoarseLevelCommunication() const - { - return *coarseLevelCommunication_; - } + const Communication& getCoarseLevelCommunication() const + { + return *coarseLevelCommunication_; + } private: - typename Operator::matrix_type::field_type prolongDamp_; - //std::shared_ptr aggregatesMap_; - Criterion criterion_; - Communication* communication_; - std::shared_ptr coarseLevelCommunication_; - std::shared_ptr coarseLevelMatrix_; + typename Operator::matrix_type::field_type prolongDamp_; + //std::shared_ptr aggregatesMap_; + Criterion criterion_; + Communication* communication_; + std::shared_ptr coarseLevelCommunication_; + std::shared_ptr coarseLevelMatrix_; }; - namespace Detail - { + namespace Detail + { + /** + * @brief A policy class for solving the coarse level system using one step of AMG. + * @tparam O The type of the linear operator used. + * @tparam S The type of the smoother used in AMG. + * @tparam C The type of the crition used for the aggregation within AMG. + * @tparam C1 The type of the information about the communication. Either + * Dune::OwnerOverlapCopyCommunication or Dune::SequentialInformation. + */ + template + class OneStepAMGCoarseSolverPolicyNoSolve + { + public: + typedef P LevelTransferPolicy; + /** @brief The type of the linear operator used. */ + typedef O Operator; + /** @brief The type of the range and domain of the operator. */ + typedef typename O::range_type X; + /** @brief The type of the crition used for the aggregation within AMG.*/ + typedef C Criterion; + /** @brief The type of the communication used for AMG.*/ + typedef typename P::ParallelInformation Communication; + /** @brief The type of the smoother used in AMG. */ + typedef S Smoother; + /** @brief The type of the arguments used for constructing the smoother. */ + typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; + /** @brief The type of the AMG construct on the coarse level.*/ + typedef Dune::Amg::AMGCPR AMGType; + /** + * @brief Constructs the coarse solver policy. + * @param args The arguments used for constructing the smoother. + * @param c The crition used for the aggregation within AMG. + */ + OneStepAMGCoarseSolverPolicyNoSolve(const CPRParameter* param, const SmootherArgs& args, const Criterion& c) + : param_(param), smootherArgs_(args), criterion_(c) + {} + /** @brief Copy constructor. */ + OneStepAMGCoarseSolverPolicyNoSolve(const OneStepAMGCoarseSolverPolicyNoSolve& other) + : param_(other.param_), coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_), + criterion_(other.criterion_) + {} + private: + /** + * @brief A wrapper that makes an inverse operator out of AMG. + * + * The operator will use one step of AMG to approximately solve + * the coarse level system. + */ + struct AMGInverseOperator : public Dune::InverseOperator + { + AMGInverseOperator(const CPRParameter* param, + typename AMGType::Operator& op, + const Criterion& crit, + const typename AMGType::SmootherArgs& args, + const Communication& comm) + : param_(param), amg_(), op_(op), crit_(crit), args_(args), comm_(comm) + { + amg_.reset(new AMGType(op, crit,args, comm)); + } + + void updateAmgPreconditioner(typename AMGType::Operator& op){ + //op_ = op; + //amg_->recalculateHierarchy(); + amg_->updateSolver(crit_, op, comm_); + //amg_.reset(new AMGType(op, crit_,args_, comm_)); + //amg_->recalculateGalerkin(); + } +#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; + } +#endif + + + void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override + { + DUNE_UNUSED_PARAMETER(reduction); + DUNE_UNUSED_PARAMETER(res); +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + auto sp = Dune::createScalarProduct(comm_, op_.category()); +#else + using Chooser = Dune::ScalarProductChooser; + auto sp = Chooser::construct(comm_); +#endif + Dune::Preconditioner* prec = amg_.get(); + // Linear solver parameters + const double tolerance = param_->cpr_solver_tol_; + const int maxit = param_->cpr_max_ell_iter_; + const int verbosity = ( param_->cpr_solver_verbose_ && + comm_.communicator().rank()==0 ) ? 1 : 0; + if ( param_->cpr_ell_solvetype_ == 0 ) + { + // Category of preconditioner will be checked at compile time. Therefore we need + // to cast to the derived class +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::BiCGSTABSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); +#else + Dune::BiCGSTABSolver solver(const_cast(op_), *sp, + reinterpret_cast(*prec), + tolerance, maxit, verbosity); +#endif + solver.apply(x,b,res); + + } + else if (param_->cpr_ell_solvetype_ == 1) + { + // Category of preconditioner will be checked at compile time. Therefore we need + // to cast to the derived class +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::CGSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); +#else + Dune::CGSolver solver(const_cast(op_), *sp, + reinterpret_cast(*prec), + tolerance, maxit, verbosity); +#endif + solver.apply(x,b,res); + } + else + { + // X v(x); + // prec->pre(x,b); + // op_->applyscaleadd(-1,x,b); + // v = 0; + // prec->apply(v,b); + // x += v; + // op_->applyscaleadd(-1,x,b); + // prec->post(x); +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + Dune::LoopSolver solver(const_cast(op_), *sp, *prec, + tolerance, maxit, verbosity); +#else + Dune::LoopSolver solver(const_cast(op_), *sp, + reinterpret_cast(*prec), + tolerance, maxit, verbosity); +#endif + solver.apply(x,b,res); + } +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) + +#else + delete sp; +#endif + } + + void apply(X& x, X& b, Dune::InverseOperatorResult& res) override + { + return apply(x,b,1e-8,res); + } + + ~AMGInverseOperator() + {} + AMGInverseOperator(const AMGInverseOperator& other) + : x_(other.x_), amg_(other.amg_) + { + } + private: + const CPRParameter* param_; + X x_; + std::unique_ptr amg_; + //std::unique_ptr op_; + typename AMGType::Operator& op_; + Criterion crit_; + typename AMGType::SmootherArgs args_; + const Communication& comm_; + }; + + public: + /** @brief The type of solver constructed for the coarse level. */ + typedef AMGInverseOperator CoarseLevelSolver; + + /** + * @brief Constructs a coarse level solver. + * + * @param transferPolicy The policy describing the transfer between levels. + * @return A pointer to the constructed coarse level solver. + */ + template + void setCoarseOperator(LTP& transferPolicy){ + coarseOperator_= transferPolicy.getCoarseLevelOperator(); + } + template + CoarseLevelSolver* createCoarseLevelSolver(LTP& transferPolicy) + { + coarseOperator_=transferPolicy.getCoarseLevelOperator(); + const LevelTransferPolicy& transfer = + reinterpret_cast(transferPolicy); + AMGInverseOperator* inv = new AMGInverseOperator(param_, + *coarseOperator_, + criterion_, + smootherArgs_, + transfer.getCoarseLevelCommunication()); + + return inv; //std::shared_ptr >(inv); + + } + //void recalculateGalerkin(){ + // coarseOperator_.recalculateHierarchy(); + //} + private: + /** @brief The coarse level operator. */ + std::shared_ptr coarseOperator_; + /** @brief The parameters for the CPR preconditioner. */ + const CPRParameter* param_; + /** @brief The arguments used to construct the smoother. */ + SmootherArgs smootherArgs_; + /** @brief The coarsening criterion. */ + Criterion criterion_; + }; + + + } // end namespace Detail + + /** - * @brief A policy class for solving the coarse level system using one step of AMG. - * @tparam O The type of the linear operator used. - * @tparam S The type of the smoother used in AMG. - * @tparam C The type of the crition used for the aggregation within AMG. - * @tparam C1 The type of the information about the communication. Either - * Dune::OwnerOverlapCopyCommunication or Dune::SequentialInformation. + * \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 OneStepAMGCoarseSolverPolicyNoSolve + template + class BlackoilAmgCpr + : public Dune::Preconditioner { public: - typedef P LevelTransferPolicy; - /** @brief The type of the linear operator used. */ - typedef O Operator; - /** @brief The type of the range and domain of the operator. */ - typedef typename O::range_type X; - /** @brief The type of the crition used for the aggregation within AMG.*/ - typedef C Criterion; - /** @brief The type of the communication used for AMG.*/ - typedef typename P::ParallelInformation Communication; - /** @brief The type of the smoother used in AMG. */ - typedef S Smoother; - /** @brief The type of the arguments used for constructing the smoother. */ - typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; - /** @brief The type of the AMG construct on the coarse level.*/ - typedef Dune::Amg::AMGCPR AMGType; - /** - * @brief Constructs the coarse solver policy. - * @param args The arguments used for constructing the smoother. - * @param c The crition used for the aggregation within AMG. - */ - OneStepAMGCoarseSolverPolicyNoSolve(const CPRParameter* param, const SmootherArgs& args, const Criterion& c) - : param_(param), smootherArgs_(args), criterion_(c) - {} - /** @brief Copy constructor. */ - OneStepAMGCoarseSolverPolicyNoSolve(const OneStepAMGCoarseSolverPolicyNoSolve& other) - : param_(other.param_), coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_), - criterion_(other.criterion_) - {} - private: - /** - * @brief A wrapper that makes an inverse operator out of AMG. - * - * The operator will use one step of AMG to approximately solve - * the coarse level system. - */ - struct AMGInverseOperator : public Dune::InverseOperator - { - AMGInverseOperator(const CPRParameter* param, - typename AMGType::Operator& op, - const Criterion& crit, - const typename AMGType::SmootherArgs& args, - const Communication& comm) - : param_(param), amg_(), op_(op), crit_(crit), args_(args), comm_(comm) - { - amg_.reset(new AMGType(op, crit,args, comm)); - } + /** \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; - void updateAmgPreconditioner(typename AMGType::Operator& op){ - //op_ = op; - //amg_->recalculateHierarchy(); - amg_->updateSolver(crit_, op, comm_); - //amg_.reset(new AMGType(op, crit_,args_, comm_)); - //amg_->recalculateGalerkin(); - } + 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 = + OneComponentAggregationLevelTransferPolicyCpr; + using CoarseSolverPolicy = + Detail::OneStepAMGCoarseSolverPolicyNoSolve; + 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; } -#endif - - - void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override - { - DUNE_UNUSED_PARAMETER(reduction); - DUNE_UNUSED_PARAMETER(res); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - auto sp = Dune::createScalarProduct(comm_, op_.category()); -#else - using Chooser = Dune::ScalarProductChooser; - auto sp = Chooser::construct(comm_); -#endif - Dune::Preconditioner* prec = amg_.get(); - // Linear solver parameters - const double tolerance = param_->cpr_solver_tol_; - const int maxit = param_->cpr_max_ell_iter_; - const int verbosity = ( param_->cpr_solver_verbose_ && - comm_.communicator().rank()==0 ) ? 1 : 0; - if ( param_->cpr_ell_solvetype_ == 0 ) - { - // Category of preconditioner will be checked at compile time. Therefore we need - // to cast to the derived class -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - Dune::BiCGSTABSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); -#else - Dune::BiCGSTABSolver solver(const_cast(op_), *sp, - reinterpret_cast(*prec), - tolerance, maxit, verbosity); -#endif - solver.apply(x,b,res); - - } - else if (param_->cpr_ell_solvetype_ == 1) - { - // Category of preconditioner will be checked at compile time. Therefore we need - // to cast to the derived class -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - Dune::CGSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); -#else - Dune::CGSolver solver(const_cast(op_), *sp, - reinterpret_cast(*prec), - tolerance, maxit, verbosity); -#endif - solver.apply(x,b,res); - } - else - { - // X v(x); - // prec->pre(x,b); - // op_->applyscaleadd(-1,x,b); - // v = 0; - // prec->apply(v,b); - // x += v; - // op_->applyscaleadd(-1,x,b); - // prec->post(x); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - Dune::LoopSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); -#else - Dune::LoopSolver solver(const_cast(op_), *sp, - reinterpret_cast(*prec), - tolerance, maxit, verbosity); -#endif - solver.apply(x,b,res); - } -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - #else - delete sp; -#endif - } - - void apply(X& x, X& b, Dune::InverseOperatorResult& res) override - { - return apply(x,b,1e-8,res); - } - - ~AMGInverseOperator() - {} - AMGInverseOperator(const AMGInverseOperator& other) - : x_(other.x_), amg_(other.amg_) - { - } - private: - const CPRParameter* param_; - X x_; - std::unique_ptr amg_; - //std::unique_ptr op_; - typename AMGType::Operator& op_; - Criterion crit_; - typename AMGType::SmootherArgs args_; - const Communication& comm_; - }; - - public: - /** @brief The type of solver constructed for the coarse level. */ - typedef AMGInverseOperator CoarseLevelSolver; - - /** - * @brief Constructs a coarse level solver. - * - * @param transferPolicy The policy describing the transfer between levels. - * @return A pointer to the constructed coarse level solver. - */ - template - void setCoarseOperator(LTP& transferPolicy){ - coarseOperator_= transferPolicy.getCoarseLevelOperator(); - } - template - CoarseLevelSolver* createCoarseLevelSolver(LTP& transferPolicy) - { - coarseOperator_=transferPolicy.getCoarseLevelOperator(); - const LevelTransferPolicy& transfer = - reinterpret_cast(transferPolicy); - AMGInverseOperator* inv = new AMGInverseOperator(param_, - *coarseOperator_, - criterion_, - smootherArgs_, - transfer.getCoarseLevelCommunication()); - - return inv; //std::shared_ptr >(inv); - - } - //void recalculateGalerkin(){ - // coarseOperator_.recalculateHierarchy(); - //} - private: - /** @brief The coarse level operator. */ - std::shared_ptr coarseOperator_; - /** @brief The parameters for the CPR preconditioner. */ - const CPRParameter* param_; - /** @brief The arguments used to construct the smoother. */ - SmootherArgs smootherArgs_; - /** @brief The coarsening criterion. */ - Criterion criterion_; - }; - - - } // end namespace Detail - - - /** - * \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 first - * coarse level system, either by simply extracting the coupling between the components at COMPONENT_INDEX - * in the matrix blocks or by extracting them and 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 the 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 CoarseSmoother = typename Detail::ScalarType::value; - using FineCriterion = - typename Detail::OneComponentCriterionType::value; - using CoarseCriterion = typename Detail::ScalarType::value; - using LevelTransferPolicy = - OneComponentAggregationLevelTransferPolicyCpr; - using CoarseSolverPolicy = - Detail::OneStepAMGCoarseSolverPolicyNoSolve; - 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 - }; + // 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::scaleMatrixDRSPtr(fineOperator, comm, - COMPONENT_INDEX, weights_, param)), - scaledMatrixOperator_(Detail::createOperatorPtr(fineOperator, *scaledMatrix_, comm)), - smoother_(Detail::constructSmoother(*scaledMatrixOperator_, - smargs, comm)), - levelTransferPolicy_(criterion, comm), - coarseSolverPolicy_(¶m, smargs, criterion), - twoLevelMethod_(*scaledMatrixOperator_, - smoother_, - levelTransferPolicy_, - coarseSolverPolicy_, 0, 1) - { - } - void updatePreconditioner(const typename TwoLevelMethod::FineDomainType& weights, - const Operator& fineOperator, - const SmootherArgs& smargs, - const Communication& comm){ - weights_ = weights; - *scaledMatrix_ = *Detail::scaleMatrixDRSPtr(fineOperator, comm, - COMPONENT_INDEX, weights_, param_); - //*scaledMatrixOperator_ = *Detail::createOperatorPtr(fineOperator,*scaledMatrix_,comm); - smoother_ .reset(Detail::constructSmoother(*scaledMatrixOperator_, - smargs, comm)); - twoLevelMethod_.updatePreconditioner(*scaledMatrixOperator_, - smoother_, - coarseSolverPolicy_); - } + /** + * \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), + coarseSolverPolicy_(¶m, smargs, criterion), + twoLevelMethod_(*scaledMatrixOperator_, + smoother_, + levelTransferPolicy_, + coarseSolverPolicy_, 0, 1) + { + } + + void updatePreconditioner(const typename TwoLevelMethod::FineDomainType& weights, + const Operator& fineOperator, + const SmootherArgs& smargs, + const Communication& comm) + { + weights_ = weights; + *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 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 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_; - typename TwoLevelMethod::FineDomainType weights_;//make copy - std::unique_ptr scaledMatrix_; - std::unique_ptr scaledMatrixOperator_; - //Operator scaledMatrixOperator_; - //std::tuple, Operator> - std::shared_ptr smoother_; - LevelTransferPolicy levelTransferPolicy_; - CoarseSolverPolicy coarseSolverPolicy_; - TwoLevelMethod twoLevelMethod_; - //BlockVector weights_; - }; + 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_; + typename TwoLevelMethod::FineDomainType weights_; //make copy + std::unique_ptr scaledMatrix_; + std::unique_ptr scaledMatrixOperator_; + std::shared_ptr smoother_; + LevelTransferPolicy levelTransferPolicy_; + CoarseSolverPolicy coarseSolverPolicy_; + TwoLevelMethod twoLevelMethod_; + }; } // end namespace Opm #endif From 032701ffc90d750584c8945bf15c9c71d3819c10 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 Apr 2019 18:07:46 +0200 Subject: [PATCH 12/25] Silence unused argument warnings. --- opm/autodiff/amgcpr.hh | 2 +- opm/autodiff/twolevelmethodcpr.hh | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/amgcpr.hh b/opm/autodiff/amgcpr.hh index 0d871e782..f5ae6aa98 100644 --- a/opm/autodiff/amgcpr.hh +++ b/opm/autodiff/amgcpr.hh @@ -394,7 +394,7 @@ namespace Dune template template - void AMGCPR::updateSolver(C& criterion, Operator& /* matrix */, const PI& pinfo) + void AMGCPR::updateSolver(C& /* criterion */, Operator& /* matrix */, const PI& /* pinfo */) { Timer watch; smoothers_.reset(new Hierarchy); diff --git a/opm/autodiff/twolevelmethodcpr.hh b/opm/autodiff/twolevelmethodcpr.hh index f6415d24a..4c23ce681 100644 --- a/opm/autodiff/twolevelmethodcpr.hh +++ b/opm/autodiff/twolevelmethodcpr.hh @@ -450,9 +450,11 @@ public: delete policy_; delete coarseSolver_; } - void updatePreconditioner(FineOperatorType& op, + + void updatePreconditioner(FineOperatorType& /* op */, std::shared_ptr smoother, - CoarseLevelSolverPolicy& coarsePolicy){ + CoarseLevelSolverPolicy& coarsePolicy) + { //assume new matrix is not reallocated the new precondition should anyway be made smoother_ = smoother; if (coarseSolver_) { @@ -466,7 +468,6 @@ public: } } - void pre(FineDomainType& x, FineRangeType& b) { smoother_->pre(x,b); From 6d18c2849da897bbad267ff84c84fbafdd81971f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 Apr 2019 18:08:07 +0200 Subject: [PATCH 13/25] Allow updating the AMG preconditioner. This adds the significant modifications from the BlackoilAmgCpr experimental class to BlackoilAmg. --- opm/autodiff/BlackoilAmg.hpp | 53 ++++++++++++++++++++++++++++-------- 1 file changed, 42 insertions(+), 11 deletions(-) diff --git a/opm/autodiff/BlackoilAmg.hpp b/opm/autodiff/BlackoilAmg.hpp index fc6a649de..2560f69ad 100644 --- a/opm/autodiff/BlackoilAmg.hpp +++ b/opm/autodiff/BlackoilAmg.hpp @@ -23,6 +23,8 @@ #include #include #include +#include +#include #include #include #include @@ -89,6 +91,7 @@ createOperator(const Dune::MatrixAdapter&, const M& matrix, const T&) { return std::make_unique< Dune::MatrixAdapter >(matrix); } + /** * \brief Creates an OverlappingSchwarzOperator as an operator, storing it in a unique_ptr. * @@ -344,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. @@ -372,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_ ) { @@ -388,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 { @@ -523,6 +531,7 @@ private: std::unique_ptr amg_; std::unique_ptr smoother_; const typename AMGType::Operator& op_; + Criterion crit_; const Communication& comm_; }; @@ -548,8 +557,13 @@ public: smootherArgs_, transfer.getCoarseLevelCommunication()); - return inv; //std::shared_ptr >(inv); + return inv; + } + template + void setCoarseOperator(LTP& transferPolicy) + { + coarseOperator_= transferPolicy.getCoarseLevelOperator(); } private: @@ -694,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: @@ -777,7 +791,7 @@ public: { coarseLevelCommunication_->freeGlobalLookup(); } - calculateCoarseEntries(fineOperator.getmat()); + calculateCoarseEntriesWithAggregatesMap(fineOperator); } else { @@ -817,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) { @@ -840,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 @@ -968,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 From 3d9c77638497d042e9fecc69f8db6efbfd76a87a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 Apr 2019 18:22:12 +0200 Subject: [PATCH 14/25] Disable most of the file, as the base version is sufficient. All features have been added to BlackoilAmg.hpp. --- opm/autodiff/BlackoilAmgCpr.hpp | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/opm/autodiff/BlackoilAmgCpr.hpp b/opm/autodiff/BlackoilAmgCpr.hpp index 25adc1389..2a9b1b644 100644 --- a/opm/autodiff/BlackoilAmgCpr.hpp +++ b/opm/autodiff/BlackoilAmgCpr.hpp @@ -38,7 +38,7 @@ namespace Opm { - + /* template class OneComponentAggregationLevelTransferPolicyCpr : public Dune::Amg::LevelTransferPolicyCpr::value> @@ -182,7 +182,8 @@ namespace Opm std::shared_ptr coarseLevelCommunication_; std::shared_ptr coarseLevelMatrix_; }; - + */ +#if 0 namespace Detail { /** @@ -399,7 +400,7 @@ namespace Opm } // end namespace Detail - +#endif // if 0 /** * \brief An algebraic twolevel or multigrid approach for solving blackoil (supports CPR with and without AMG) @@ -442,16 +443,16 @@ namespace Opm typename Detail::OneComponentCriterionType::value; using CoarseCriterion = typename Detail::ScalarType::value; using LevelTransferPolicy = - OneComponentAggregationLevelTransferPolicyCpr; + OneComponentAggregationLevelTransferPolicy; using CoarseSolverPolicy = - Detail::OneStepAMGCoarseSolverPolicyNoSolve; + Detail::OneStepAMGCoarseSolverPolicy; using TwoLevelMethod = Dune::Amg::TwoLevelMethodCpr(*scaledMatrixOperator_, smargs, comm)), - levelTransferPolicy_(criterion, comm), + levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_), coarseSolverPolicy_(¶m, smargs, criterion), twoLevelMethod_(*scaledMatrixOperator_, smoother_, From 0e5a5c52664886d6f7b313969162acae06b36178 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 Apr 2019 18:23:09 +0200 Subject: [PATCH 15/25] Minor cleanup. --- opm/autodiff/ISTLSolverEbosCpr.hpp | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/opm/autodiff/ISTLSolverEbosCpr.hpp b/opm/autodiff/ISTLSolverEbosCpr.hpp index 0939788aa..ec3878738 100644 --- a/opm/autodiff/ISTLSolverEbosCpr.hpp +++ b/opm/autodiff/ISTLSolverEbosCpr.hpp @@ -48,7 +48,7 @@ namespace Opm template class ISTLSolverEbosCpr : public ISTLSolverEbos { - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; @@ -65,15 +65,15 @@ namespace Opm //typedef typename GridView::template Codim<0>::Entity Element; //typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; - enum { pressureEqnIndex = BlackOilDefaultIndexTraits::waterCompIdx }; + enum { pressureEqnIndex = BlackOilDefaultIndexTraits::waterCompIdx }; enum { pressureVarIndex = Indices::pressureSwitchIdx }; static const int numEq = Indices::numEq; - typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false> OperatorSerial; - typedef ISTLSolverEbos SuperClass; - typedef Dune::Amg::SequentialInformation POrComm; - //typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType; - typedef Dune::MatrixAdapter MatrixAdapter; + typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false> OperatorSerial; + typedef ISTLSolverEbos SuperClass; + typedef Dune::Amg::SequentialInformation POrComm; + //typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType; + typedef Dune::MatrixAdapter MatrixAdapter; #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) @@ -89,7 +89,7 @@ namespace Opm using CouplingMetric = Opm::Amg::Element; using CritBase = Dune::Amg::SymmetricCriterion; using Criterion = Dune::Amg::CoarsenCriterion; - typedef BlackoilAmgCpr BLACKOILAMG; + typedef BlackoilAmgCpr BlackoilAmgType; public: @@ -232,8 +232,8 @@ namespace Opm //criterion.setGamma(1); // //1 V cycle 2 WW // Since DUNE 2.2 we also need to pass the smoother args instead of steps directly - typedef typename BLACKOILAMG::Smoother Smoother; - typedef typename BLACKOILAMG::Smoother Smoother; + typedef typename BlackoilAmgType::Smoother Smoother; + typedef typename BlackoilAmgType::Smoother Smoother; typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; SmootherArgs smootherArgs; smootherArgs.iterations = 1; @@ -263,7 +263,7 @@ namespace Opm } if( update_preconditioner or (amg_== 0) ){ - amg_.reset( new BLACKOILAMG( params, this->weights_, opARef, criterion, smootherArgs, comm ) ); + amg_.reset( new BlackoilAmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) ); }else{ if(this->parameters_.cpr_solver_verbose_){ std::cout << " Only update amg solver " << std::endl; @@ -340,8 +340,8 @@ namespace Opm #endif std::unique_ptr< MatrixAdapter > opA_; std::unique_ptr< OperatorSerial > opASerial_; - std::unique_ptr< BLACKOILAMG > amg_; - SPPointer sp_; + std::unique_ptr< BlackoilAmgType > amg_; + SPPointer sp_; std::shared_ptr< Dune::BiCGSTABSolver > linsolve_; //std::shared_ptr< Dune::LinearOperator > op_; //std::shared_ptr< Dune::Preconditioner > prec_; From 416cc93f964232428627211745e77d769a611923 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 8 Apr 2019 09:56:13 +0200 Subject: [PATCH 16/25] Remove code that is no longer needed. --- opm/autodiff/BlackoilAmgCpr.hpp | 363 -------------------------------- 1 file changed, 363 deletions(-) diff --git a/opm/autodiff/BlackoilAmgCpr.hpp b/opm/autodiff/BlackoilAmgCpr.hpp index 2a9b1b644..d3f61d393 100644 --- a/opm/autodiff/BlackoilAmgCpr.hpp +++ b/opm/autodiff/BlackoilAmgCpr.hpp @@ -38,370 +38,7 @@ namespace Opm { - /* - template - class OneComponentAggregationLevelTransferPolicyCpr - : public Dune::Amg::LevelTransferPolicyCpr::value> - { - typedef Dune::Amg::AggregatesMap AggregatesMap; - public: - using CoarseOperator = typename Detail::ScalarType::value; - typedef Dune::Amg::LevelTransferPolicy FatherType; - typedef Communication ParallelInformation; - - public: - OneComponentAggregationLevelTransferPolicyCpr(const Criterion& crit, const Communication& comm) - : criterion_(crit), communication_(&const_cast(comm)) - {} - void createCoarseLevelSystem(const Operator& fineOperator) - { - prolongDamp_ = 1; - - using CoarseMatrix = typename CoarseOperator::matrix_type; - const auto& fineLevelMatrix = fineOperator.getmat(); - coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise)); - auto createIter = coarseLevelMatrix_->createbegin(); - - for ( const auto& row: fineLevelMatrix ) - { - for ( auto col = row.begin(), cend = row.end(); col != cend; ++col) - { - createIter.insert(col.index()); - } - ++createIter; - } - - auto coarseRow = coarseLevelMatrix_->begin(); - for ( const auto& row: fineLevelMatrix ) - { - auto coarseCol = coarseRow->begin(); - - for ( auto col = row.begin(), cend = row.end(); col != cend; ++col, ++coarseCol ) - { - assert( col.index() == coarseCol.index() ); - *coarseCol = (*col)[COMPONENT_INDEX][VARIABLE_INDEX]; - } - ++coarseRow; - } - coarseLevelCommunication_.reset(communication_, [](Communication*){}); - - - this->lhs_.resize(this->coarseLevelMatrix_->M()); - this->rhs_.resize(this->coarseLevelMatrix_->N()); - using OperatorArgs = typename Dune::Amg::ConstructionTraits::Arguments; - OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_); - this->operator_.reset(Dune::Amg::ConstructionTraits::construct(oargs)); - } - - // compleately unsafe!!!!!! - void calculateCoarseEntries(const Operator& fineOperator)//const M& fineMatrix) - { - 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]; - } - } - } - - - //template - // void calculateCoarseEntriesOld(const Operator& fineOperator)//const M& fineMatrix) - // { - // const auto& fineMatrix = fineOperator.getmat(); - // *coarseLevelMatrix_ = 0; - // for(auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); - // row != rowEnd; ++row) - // { - // const auto& i = (*aggregatesMap_)[row.index()]; - // if(i != AggregatesMap::ISOLATED) - // { - // for(auto entry = row->begin(), entryEnd = row->end(); - // entry != entryEnd; ++entry) - // { - // const auto& j = (*aggregatesMap_)[entry.index()]; - // if ( j != AggregatesMap::ISOLATED ) - // { - // (*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][COMPONENT_INDEX]; - // } - // } - // } - // } - // } - - void moveToCoarseLevel(const typename FatherType::FineRangeType& fine) - { - // Set coarse vector to zero - this->rhs_=0; - - auto end = fine.end(), begin=fine.begin(); - - for(auto block=begin; block != end; ++block) - { - this->rhs_[block-begin] = (*block)[COMPONENT_INDEX]; - } - - - this->lhs_=0; - } - - void moveToFineLevel(typename FatherType::FineDomainType& fine) - { - - auto end=fine.end(), begin=fine.begin(); - - for(auto block=begin; block != end; ++block) - { - (*block)[COMPONENT_INDEX] = this->lhs_[block-begin]; - } - - } - - OneComponentAggregationLevelTransferPolicyCpr* clone() const - { - return new OneComponentAggregationLevelTransferPolicyCpr(*this); - } - - const Communication& getCoarseLevelCommunication() const - { - return *coarseLevelCommunication_; - } - private: - typename Operator::matrix_type::field_type prolongDamp_; - //std::shared_ptr aggregatesMap_; - Criterion criterion_; - Communication* communication_; - std::shared_ptr coarseLevelCommunication_; - std::shared_ptr coarseLevelMatrix_; - }; - */ -#if 0 - namespace Detail - { - /** - * @brief A policy class for solving the coarse level system using one step of AMG. - * @tparam O The type of the linear operator used. - * @tparam S The type of the smoother used in AMG. - * @tparam C The type of the crition used for the aggregation within AMG. - * @tparam C1 The type of the information about the communication. Either - * Dune::OwnerOverlapCopyCommunication or Dune::SequentialInformation. - */ - template - class OneStepAMGCoarseSolverPolicyNoSolve - { - public: - typedef P LevelTransferPolicy; - /** @brief The type of the linear operator used. */ - typedef O Operator; - /** @brief The type of the range and domain of the operator. */ - typedef typename O::range_type X; - /** @brief The type of the crition used for the aggregation within AMG.*/ - typedef C Criterion; - /** @brief The type of the communication used for AMG.*/ - typedef typename P::ParallelInformation Communication; - /** @brief The type of the smoother used in AMG. */ - typedef S Smoother; - /** @brief The type of the arguments used for constructing the smoother. */ - typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; - /** @brief The type of the AMG construct on the coarse level.*/ - typedef Dune::Amg::AMGCPR AMGType; - /** - * @brief Constructs the coarse solver policy. - * @param args The arguments used for constructing the smoother. - * @param c The crition used for the aggregation within AMG. - */ - OneStepAMGCoarseSolverPolicyNoSolve(const CPRParameter* param, const SmootherArgs& args, const Criterion& c) - : param_(param), smootherArgs_(args), criterion_(c) - {} - /** @brief Copy constructor. */ - OneStepAMGCoarseSolverPolicyNoSolve(const OneStepAMGCoarseSolverPolicyNoSolve& other) - : param_(other.param_), coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_), - criterion_(other.criterion_) - {} - private: - /** - * @brief A wrapper that makes an inverse operator out of AMG. - * - * The operator will use one step of AMG to approximately solve - * the coarse level system. - */ - struct AMGInverseOperator : public Dune::InverseOperator - { - AMGInverseOperator(const CPRParameter* param, - typename AMGType::Operator& op, - const Criterion& crit, - const typename AMGType::SmootherArgs& args, - const Communication& comm) - : param_(param), amg_(), op_(op), crit_(crit), args_(args), comm_(comm) - { - amg_.reset(new AMGType(op, crit,args, comm)); - } - - void updateAmgPreconditioner(typename AMGType::Operator& op){ - //op_ = op; - //amg_->recalculateHierarchy(); - amg_->updateSolver(crit_, op, comm_); - //amg_.reset(new AMGType(op, crit_,args_, comm_)); - //amg_->recalculateGalerkin(); - } -#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; - } -#endif - - - void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override - { - DUNE_UNUSED_PARAMETER(reduction); - DUNE_UNUSED_PARAMETER(res); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - auto sp = Dune::createScalarProduct(comm_, op_.category()); -#else - using Chooser = Dune::ScalarProductChooser; - auto sp = Chooser::construct(comm_); -#endif - Dune::Preconditioner* prec = amg_.get(); - // Linear solver parameters - const double tolerance = param_->cpr_solver_tol_; - const int maxit = param_->cpr_max_ell_iter_; - const int verbosity = ( param_->cpr_solver_verbose_ && - comm_.communicator().rank()==0 ) ? 1 : 0; - if ( param_->cpr_ell_solvetype_ == 0 ) - { - // Category of preconditioner will be checked at compile time. Therefore we need - // to cast to the derived class -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - Dune::BiCGSTABSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); -#else - Dune::BiCGSTABSolver solver(const_cast(op_), *sp, - reinterpret_cast(*prec), - tolerance, maxit, verbosity); -#endif - solver.apply(x,b,res); - - } - else if (param_->cpr_ell_solvetype_ == 1) - { - // Category of preconditioner will be checked at compile time. Therefore we need - // to cast to the derived class -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - Dune::CGSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); -#else - Dune::CGSolver solver(const_cast(op_), *sp, - reinterpret_cast(*prec), - tolerance, maxit, verbosity); -#endif - solver.apply(x,b,res); - } - else - { - // X v(x); - // prec->pre(x,b); - // op_->applyscaleadd(-1,x,b); - // v = 0; - // prec->apply(v,b); - // x += v; - // op_->applyscaleadd(-1,x,b); - // prec->post(x); -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - Dune::LoopSolver solver(const_cast(op_), *sp, *prec, - tolerance, maxit, verbosity); -#else - Dune::LoopSolver solver(const_cast(op_), *sp, - reinterpret_cast(*prec), - tolerance, maxit, verbosity); -#endif - solver.apply(x,b,res); - } -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - -#else - delete sp; -#endif - } - - void apply(X& x, X& b, Dune::InverseOperatorResult& res) override - { - return apply(x,b,1e-8,res); - } - - ~AMGInverseOperator() - {} - AMGInverseOperator(const AMGInverseOperator& other) - : x_(other.x_), amg_(other.amg_) - { - } - private: - const CPRParameter* param_; - X x_; - std::unique_ptr amg_; - //std::unique_ptr op_; - typename AMGType::Operator& op_; - Criterion crit_; - typename AMGType::SmootherArgs args_; - const Communication& comm_; - }; - - public: - /** @brief The type of solver constructed for the coarse level. */ - typedef AMGInverseOperator CoarseLevelSolver; - - /** - * @brief Constructs a coarse level solver. - * - * @param transferPolicy The policy describing the transfer between levels. - * @return A pointer to the constructed coarse level solver. - */ - template - void setCoarseOperator(LTP& transferPolicy){ - coarseOperator_= transferPolicy.getCoarseLevelOperator(); - } - template - CoarseLevelSolver* createCoarseLevelSolver(LTP& transferPolicy) - { - coarseOperator_=transferPolicy.getCoarseLevelOperator(); - const LevelTransferPolicy& transfer = - reinterpret_cast(transferPolicy); - AMGInverseOperator* inv = new AMGInverseOperator(param_, - *coarseOperator_, - criterion_, - smootherArgs_, - transfer.getCoarseLevelCommunication()); - - return inv; //std::shared_ptr >(inv); - - } - //void recalculateGalerkin(){ - // coarseOperator_.recalculateHierarchy(); - //} - private: - /** @brief The coarse level operator. */ - std::shared_ptr coarseOperator_; - /** @brief The parameters for the CPR preconditioner. */ - const CPRParameter* param_; - /** @brief The arguments used to construct the smoother. */ - SmootherArgs smootherArgs_; - /** @brief The coarsening criterion. */ - Criterion criterion_; - }; - - - } // end namespace Detail -#endif // if 0 - /** * \brief An algebraic twolevel or multigrid approach for solving blackoil (supports CPR with and without AMG) * From 2cf990e4044c361f81b8115f98fb59d0caadac0d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 8 Apr 2019 09:57:22 +0200 Subject: [PATCH 17/25] Do not store copy of weights inside BlackoilAmgCpr class. --- opm/autodiff/BlackoilAmgCpr.hpp | 6 ++---- opm/autodiff/ISTLSolverEbosCpr.hpp | 2 +- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/BlackoilAmgCpr.hpp b/opm/autodiff/BlackoilAmgCpr.hpp index d3f61d393..ae6ece976 100644 --- a/opm/autodiff/BlackoilAmgCpr.hpp +++ b/opm/autodiff/BlackoilAmgCpr.hpp @@ -136,12 +136,10 @@ namespace Opm { } - void updatePreconditioner(const typename TwoLevelMethod::FineDomainType& weights, - const Operator& fineOperator, + void updatePreconditioner(const Operator& fineOperator, const SmootherArgs& smargs, const Communication& comm) { - weights_ = weights; *scaledMatrix_ = *Detail::scaleMatrixDRS(fineOperator, COMPONENT_INDEX, weights_, param_); smoother_.reset(Detail::constructSmoother(*scaledMatrixOperator_, smargs, comm)); twoLevelMethod_.updatePreconditioner(*scaledMatrixOperator_, @@ -170,7 +168,7 @@ namespace Opm private: const CPRParameter& param_; - typename TwoLevelMethod::FineDomainType weights_; //make copy + const typename TwoLevelMethod::FineDomainType& weights_; std::unique_ptr scaledMatrix_; std::unique_ptr scaledMatrixOperator_; std::shared_ptr smoother_; diff --git a/opm/autodiff/ISTLSolverEbosCpr.hpp b/opm/autodiff/ISTLSolverEbosCpr.hpp index ec3878738..9b38d785a 100644 --- a/opm/autodiff/ISTLSolverEbosCpr.hpp +++ b/opm/autodiff/ISTLSolverEbosCpr.hpp @@ -268,7 +268,7 @@ namespace Opm if(this->parameters_.cpr_solver_verbose_){ std::cout << " Only update amg solver " << std::endl; } - amg_->updatePreconditioner(this->weights_,opARef, smootherArgs, comm); + amg_->updatePreconditioner(opARef, smootherArgs, comm); } // Solve. //SuperClass::solve(linearOperator, x, istlb, *sp, *amg, result); From be86b90bac82e7c12fa0e350a83f0434b4a86322 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 8 Apr 2019 13:44:48 +0200 Subject: [PATCH 18/25] Make typed protected rather than private, enabling inheritance. --- opm/autodiff/ISTLSolverEbos.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index 0013e1578..4f29fbf24 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -178,6 +178,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; From 7664d27849877386cafc1de8a03ba92d8b946ee0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 8 Apr 2019 13:45:08 +0200 Subject: [PATCH 19/25] Clean up formatting and get most types from superclass. --- opm/autodiff/ISTLSolverEbosCpr.hpp | 403 ++++++++++++----------------- 1 file changed, 164 insertions(+), 239 deletions(-) diff --git a/opm/autodiff/ISTLSolverEbosCpr.hpp b/opm/autodiff/ISTLSolverEbosCpr.hpp index 9b38d785a..efc9f98bf 100644 --- a/opm/autodiff/ISTLSolverEbosCpr.hpp +++ b/opm/autodiff/ISTLSolverEbosCpr.hpp @@ -48,53 +48,37 @@ namespace Opm template class ISTLSolverEbosCpr : public ISTLSolverEbos { - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; - typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; - typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - typedef typename GET_PROP_TYPE(TypeTag, EclWellModel) WellModel; - typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; - typedef typename SparseMatrixAdapter::IstlMatrix Matrix; - typedef typename GET_PROP_TYPE(TypeTag, CprSmootherFine) CprSmootherFine; - typedef typename GET_PROP_TYPE(TypeTag, CprSmootherCoarse) CprSmootherCoarse; - //typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType; - //typedef typename Vector::block_type BlockVector; - //typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; - //typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager; - //typedef typename GridView::template Codim<0>::Entity Element; - //typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + 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 }; - enum { pressureEqnIndex = BlackOilDefaultIndexTraits::waterCompIdx }; - enum { pressureVarIndex = Indices::pressureSwitchIdx }; + // New properties in this subclass. + using CprSmootherFine = typename GET_PROP_TYPE(TypeTag, CprSmootherFine); + using CprSmootherCoarse = typename GET_PROP_TYPE(TypeTag, CprSmootherCoarse); + + using OperatorSerial = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false>; + using POrComm = Dune::Amg::SequentialInformation; + using MatrixAdapter = Dune::MatrixAdapter; + +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - static const int numEq = Indices::numEq; - typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false> OperatorSerial; - typedef ISTLSolverEbos SuperClass; - typedef Dune::Amg::SequentialInformation POrComm; - //typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType; - typedef Dune::MatrixAdapter MatrixAdapter; - -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - #else static constexpr int category = Dune::SolverCategory::sequential; typedef Dune::ScalarProductChooser ScalarProductChooser; #endif -//Operator MatrixOperator = Dune::MatrixAdapter - //typedef Opm::ParallelOverlappingILU0 Smoother; - typedef CprSmootherFine Smoother; - //ParallelInformation = Dune::Amg::SequentialInformation - //typedef Dune::Amg::AMG DuneAmg; - using CouplingMetric = Opm::Amg::Element; - using CritBase = Dune::Amg::SymmetricCriterion; - using Criterion = Dune::Amg::CoarsenCriterion; - typedef BlackoilAmgCpr BlackoilAmgType; - + using CouplingMetric = Opm::Amg::Element; + using CritBase = Dune::Amg::SymmetricCriterion; + using Criterion = Dune::Amg::CoarsenCriterion; + using BlackoilAmgType = BlackoilAmgCpr; public: - typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType; - static void registerParameters() { FlowLinearSolverParameters::registerParameters(); @@ -103,236 +87,177 @@ namespace Opm /// Construct a system solver. /// \param[in] parallelInformation In the case of a parallel run /// with dune-istl the information about the parallelization. - ISTLSolverEbosCpr(const Simulator& simulator) - : SuperClass(simulator) - { - } - - // nothing to clean here - void eraseMatrix() { - this->matrix_for_preconditioner_.reset(); + explicit ISTLSolverEbosCpr(const Simulator& simulator) + : SuperClass(simulator) + { } - void prepare(const SparseMatrixAdapter& M, Vector& b){ - int newton_iteration = this->simulator_.model().newtonMethod().numIterations(); - // double dt = this->simulator_.timeStepSize(); - 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(); + void prepare(const SparseMatrixAdapter& M, Vector& b) + { + int newton_iteration = this->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() ) - { - // parallel implemantation si as before - // typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true ,TypeTag> Operator; + if( this->isParallel() ) { + // parallel implemantation si as before + // typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true ,TypeTag> Operator; - // auto ebosJacIgnoreOverlap = Matrix(*(this->matrix_)); - // //remove ghost rows in local matrix - // this->makeOverlapRowsInvalid(ebosJacIgnoreOverlap); + // auto ebosJacIgnoreOverlap = Matrix(*(this->matrix_)); + // //remove ghost rows in local matrix + // this->makeOverlapRowsInvalid(ebosJacIgnoreOverlap); - // //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables - // //to be certain that correct matrix is used for preconditioning. - // Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel, - // this->parallelInformation_ ); - // assert( opA.comm() ); - // //SuperClass::solve( opA, x, *(this->rhs_), *(opA.comm()) ); - // typedef Dune::OwnerOverlapCopyCommunication& comm = *(opA.comm()); - // const size_t size = opA.getmat().N(); + // //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables + // //to be certain that correct matrix is used for preconditioning. + // Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel, + // this->parallelInformation_ ); + // assert( opA.comm() ); + // //SuperClass::solve( opA, x, *(this->rhs_), *(opA.comm()) ); + // typedef Dune::OwnerOverlapCopyCommunication& comm = *(opA.comm()); + // const size_t size = opA.getmat().N(); - // const ParallelISTLInformation& info = - // boost::any_cast( this->parallelInformation_); + // 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); - // // Construct operator, scalar product and vectors needed. - // Dune::InverseOperatorResult result; - // SuperClass::constructPreconditionerAndSolve(opA, x, *(this->rhs_), comm, result); - // SuperClass::checkConvergence(result); + // // 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); + // // Construct operator, scalar product and vectors needed. + // Dune::InverseOperatorResult result; + // SuperClass::constructPreconditionerAndSolve(opA, x, *(this->rhs_), comm, result); + // SuperClass::checkConvergence(result); - } - else + } else #endif - { - const WellModel& wellModel = this->simulator_.problem().wellModel(); - //typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false ,TypeTag> OperatorSerial; - opASerial_.reset(new OperatorSerial(*(this->matrix_), *(this->matrix_), wellModel)); - - //Dune::Amg::SequentialInformation info; - typedef Dune::Amg::SequentialInformation POrComm; - POrComm parallelInformation_arg; - typedef OperatorSerial LinearOperator; - - //SuperClass::constructPreconditionerAndSolve(opA, x, *(this->rhs_), info, result); - - + { + const WellModel& wellModel = this->simulator_.problem().wellModel(); + opASerial_.reset(new OperatorSerial(*(this->matrix_), *(this->matrix_), wellModel)); + typedef Dune::Amg::SequentialInformation POrComm; + 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); + 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); + 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 - Vector& istlb = *(this->rhs_); - parallelInformation_arg.copyOwnerToAll(istlb, istlb); - - + Vector& istlb = *(this->rhs_); + parallelInformation_arg.copyOwnerToAll(istlb, istlb); - if( ! std::is_same< LinearOperator, MatrixAdapter > :: value ) - { - // create new operator in case linear operator and matrix operator differ - opA_.reset( new MatrixAdapter( opASerial_->getmat()));//, parallelInformation_arg ) ); + if( ! std::is_same< LinearOperator, MatrixAdapter > :: value ) { + // create new operator in case linear operator and matrix operator differ + opA_.reset( new MatrixAdapter( opASerial_->getmat()));//, parallelInformation_arg ) ); } - const double relax = this->parameters_.ilu_relaxation_; - const MILU_VARIANT ilu_milu = this->parameters_.ilu_milu_; - using Matrix = typename MatrixAdapter::matrix_type; - //using CouplingMetric = Dune::Amg::Diagonal; - //using CritBase = Dune::Amg::SymmetricCriterion; - //using Criterion = Dune::Amg::CoarsenCriterion; - //using AMG = typename ISTLUtility - // ::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureIndex >::AMG; - - //std::unique_ptr< AMG > amg; - // Construct preconditioner. - //Criterion crit(15, 2000); - //SuperClass::constructAMGPrecond< Criterion >( linearOperator, parallelInformation_arg, amg, opA, relax, ilu_milu ); - // ISTLUtility::template createAMGPreconditionerPointer( *opA_, - // relax, - // parallelInformation_arg, - // amg_, - // this->parameters_, - // this->weights_ ); - //using AMG = BlackoilAmg; - POrComm& comm = parallelInformation_arg; - const int verbosity = ( this->parameters_.cpr_solver_verbose_ && - comm.communicator().rank()==0 ) ? 1 : 0; + const double relax = this->parameters_.ilu_relaxation_; + const MILU_VARIANT ilu_milu = this->parameters_.ilu_milu_; + using Matrix = typename MatrixAdapter::matrix_type; + POrComm& comm = parallelInformation_arg; + 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 + // 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 - typedef typename BlackoilAmgType::Smoother Smoother; - typedef typename BlackoilAmgType::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); - //ISTLUtility::setILUParameters(smootherArgs, params); - //smootherArgs.setN(params.cpr_ilu_n_); smootherArgs.setMilu(params.cpr_ilu_milu_); + // Since DUNE 2.2 we also need to pass the smoother args instead of steps directly + typedef typename BlackoilAmgType::Smoother Smoother; + typedef typename BlackoilAmgType::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); - MatrixAdapter& opARef = *opA_; - int newton_iteration = this->simulator_.model().newtonMethod().numIterations(); - double dt = this->simulator_.timeStepSize(); - bool update_preconditioner = false; + MatrixAdapter& opARef = *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 (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 BlackoilAmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) ); - }else{ - if(this->parameters_.cpr_solver_verbose_){ - std::cout << " Only update amg solver " << std::endl; - } - amg_->updatePreconditioner(opARef, smootherArgs, comm); - } - // Solve. - //SuperClass::solve(linearOperator, x, istlb, *sp, *amg, result); - //references seems to do something els than refering + if ( update_preconditioner or (amg_== 0) ) { + amg_.reset( new BlackoilAmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) ); + } else { + if (this->parameters_.cpr_solver_verbose_) { + std::cout << " Only update amg solver " << std::endl; + } + amg_->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_; - } + int verbosity_linsolve = 0; + if (comm.communicator().rank() == 0) { + verbosity_linsolve = this->parameters_.linear_solver_verbosity_; + } - LinearOperator& opASerialRef = *opASerial_; - linsolve_.reset(new Dune::BiCGSTABSolver(opASerialRef, *sp_, *amg_, - this->parameters_.linear_solver_reduction_, - this->parameters_.linear_solver_maxiter_, - verbosity_linsolve)); - - - } + LinearOperator& opASerialRef = *opASerial_; + linsolve_.reset(new Dune::BiCGSTABSolver(opASerialRef, *sp_, *amg_, + this->parameters_.linear_solver_reduction_, + this->parameters_.linear_solver_maxiter_, + verbosity_linsolve)); + } } - bool solve(Vector& x) { - //SuperClass::solve(x); - if( this->isParallel() ){ - // for now only call the superclass - bool converged = SuperClass::solve(x); - return converged; - }else{ - // 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); - } - } + bool solve(Vector& x) + { + if (this->isParallel()) { + // for now only call the superclass + bool converged = SuperClass::solve(x); + return converged; + } else { + // 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_; } - /// Solve the system of linear equations Ax = b, with A being the - /// combined derivative matrix of the residual and b - /// being the residual itself. - /// \param[in] residual residual object containing A and b. - /// \return the solution x - - /// \copydoc NewtonIterationBlackoilInterface::iterations - int iterations () const { return this->iterations_; } - - /// \copydoc NewtonIterationBlackoilInterface::parallelInformation - const boost::any& parallelInformation() const { return this->parallelInformation_; } - protected: - - //using Matrix = typename MatrixAdapter::matrix_type; - //using CouplingMetric = Dune::Amg::Diagonal; - //using CritBase = Dune::Amg::SymmetricCriterion; - //using Criterion = Dune::Amg::CoarsenCriterion; - //using AMG = typename ISTLUtility - // ::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureIndex >::AMG; - //Operator MatrixOperator = Dune::MatrixAdapter - //Smoother ParallelOverLappingILU0 - //ParallelInformation = Dune::Amg::SequentialInformation + #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) typedef std::shared_ptr< Dune::ScalarProduct > SPPointer; #else From 07e495c9da0c3aa4d247be7dd93a8692d28aae7f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 9 Apr 2019 10:37:09 +0200 Subject: [PATCH 20/25] More minor cleanup. --- opm/autodiff/ISTLSolverEbosCpr.hpp | 63 +++++++++++++++--------------- 1 file changed, 32 insertions(+), 31 deletions(-) diff --git a/opm/autodiff/ISTLSolverEbosCpr.hpp b/opm/autodiff/ISTLSolverEbosCpr.hpp index efc9f98bf..6f52b06c3 100644 --- a/opm/autodiff/ISTLSolverEbosCpr.hpp +++ b/opm/autodiff/ISTLSolverEbosCpr.hpp @@ -95,7 +95,7 @@ namespace Opm void prepare(const SparseMatrixAdapter& M, Vector& b) { int newton_iteration = this->simulator_.model().newtonMethod().numIterations(); - if( newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_) ) { + if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) { SuperClass::matrix_.reset(new Matrix(M.istlMatrix())); } else { *SuperClass::matrix_ = M.istlMatrix(); @@ -106,33 +106,34 @@ namespace Opm #if HAVE_MPI if( this->isParallel() ) { - // parallel implemantation si as before - // typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true ,TypeTag> Operator; - - // auto ebosJacIgnoreOverlap = Matrix(*(this->matrix_)); - // //remove ghost rows in local matrix - // this->makeOverlapRowsInvalid(ebosJacIgnoreOverlap); +#if 0 + // Copied from superclass. + typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true> Operator; + auto ebosJacIgnoreOverlap = Matrix(*(this->matrix_)); + //remove ghost rows in local matrix + this->makeOverlapRowsInvalid(ebosJacIgnoreOverlap); - // //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables - // //to be certain that correct matrix is used for preconditioning. - // Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel, - // this->parallelInformation_ ); - // assert( opA.comm() ); - // //SuperClass::solve( opA, x, *(this->rhs_), *(opA.comm()) ); - // typedef Dune::OwnerOverlapCopyCommunication& comm = *(opA.comm()); - // const size_t size = opA.getmat().N(); + //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables + //to be certain that correct matrix is used for preconditioning. + Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel, + this->parallelInformation_ ); + assert( opA.comm() ); + //SuperClass::solve( opA, x, *(this->rhs_), *(opA.comm()) ); + Dune::OwnerOverlapCopyCommunication* comm = opA.comm(); + const size_t size = opA.getmat().N(); - // const ParallelISTLInformation& info = - // boost::any_cast( this->parallelInformation_); + 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); - // // Construct operator, scalar product and vectors needed. - // Dune::InverseOperatorResult result; - // SuperClass::constructPreconditionerAndSolve(opA, x, *(this->rhs_), comm, result); - // SuperClass::checkConvergence(result); + // 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); + // Construct operator, scalar product and vectors needed. + Dune::InverseOperatorResult result; + SuperClass::constructPreconditionerAndSolve(opA, x, *(this->rhs_), *comm, result); + SuperClass::checkConvergence(result); +#endif // if 0 } else #endif @@ -156,7 +157,7 @@ namespace Opm #endif Vector& istlb = *(this->rhs_); parallelInformation_arg.copyOwnerToAll(istlb, istlb); - + if( ! std::is_same< LinearOperator, MatrixAdapter > :: value ) { // create new operator in case linear operator and matrix operator differ opA_.reset( new MatrixAdapter( opASerial_->getmat()));//, parallelInformation_arg ) ); @@ -168,10 +169,10 @@ namespace Opm POrComm& comm = parallelInformation_arg; 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; + // 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); @@ -181,7 +182,7 @@ namespace Opm //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 typedef typename BlackoilAmgType::Smoother Smoother; typedef typename BlackoilAmgType::Smoother Smoother; @@ -254,7 +255,7 @@ namespace Opm } return this->converged_; } - + protected: From 6538b59d9ed694f1899c632352d6f281ee3a4d2b Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 9 Apr 2019 14:56:00 +0200 Subject: [PATCH 21/25] Added parallel CPR that reuses the matrix hierarchy. Currently this just compiles but still segfaults in parallel. --- flow/flow_blackoil_dunecpr.cpp | 29 --- opm/autodiff/ISTLSolverEbos.hpp | 3 +- opm/autodiff/ISTLSolverEbosCpr.hpp | 334 ++++++++++++++++------------- 3 files changed, 186 insertions(+), 180 deletions(-) diff --git a/flow/flow_blackoil_dunecpr.cpp b/flow/flow_blackoil_dunecpr.cpp index 2c8468078..efb95b9dc 100644 --- a/flow/flow_blackoil_dunecpr.cpp +++ b/flow/flow_blackoil_dunecpr.cpp @@ -27,8 +27,6 @@ BEGIN_PROPERTIES NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem)); NEW_PROP_TAG(FluidState); -NEW_PROP_TAG(CprSmootherFine); -NEW_PROP_TAG(CprSmootherCoarse); //SET_TYPE_PROP(EclBaseProblem, Problem, Ewoms::EclProblem); SET_PROP(EclFlowProblemSimple, FluidState) { @@ -47,33 +45,6 @@ SET_PROP(EclFlowProblemSimple, FluidState) //typedef Opm::BlackOilFluidSystemSimple type; typedef Opm::BlackOilFluidState type; }; -SET_PROP(EclFlowProblemSimple, CprSmootherFine) - { - private: - typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; - typedef typename SparseMatrixAdapter::IstlMatrix Matrix; - typedef Dune::Amg::SequentialInformation POrComm; - - public: - typedef Opm::ParallelOverlappingILU0 type; - //typedef Dune::SeqILU0 type; -}; - -SET_PROP(EclFlowProblemSimple, CprSmootherCoarse) - { - private: - typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; - typedef typename SparseMatrixAdapter::IstlMatrix Matrix; - typedef Dune::Amg::SequentialInformation POrComm; - - public: - typedef Opm::ParallelOverlappingILU0 type; - //typedef Dune::SeqILU0 type; -}; SET_BOOL_PROP(EclFlowProblemSimple,MatrixAddWellContributions,true); SET_INT_PROP(EclFlowProblemSimple,LinearSolverVerbosity,1); diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index 4f29fbf24..e329d1f34 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. @@ -167,7 +168,7 @@ public: protected: const matrix_type& A_ ; const matrix_type& A_for_precond_ ; - const WellModel& wellMod_; + const WellModel& wellMod_; std::unique_ptr< communication_type > comm_; }; diff --git a/opm/autodiff/ISTLSolverEbosCpr.hpp b/opm/autodiff/ISTLSolverEbosCpr.hpp index 6f52b06c3..b6e8e766e 100644 --- a/opm/autodiff/ISTLSolverEbosCpr.hpp +++ b/opm/autodiff/ISTLSolverEbosCpr.hpp @@ -24,10 +24,6 @@ #include #include #include -BEGIN_PROPERTIES -NEW_PROP_TAG(CprSmootherFine); -NEW_PROP_TAG(CprSmootherCoarse); -END_PROPERTIES namespace Opm { @@ -60,23 +56,25 @@ namespace Opm enum { pressureVarIndex = SuperClass::pressureVarIndex }; // New properties in this subclass. - using CprSmootherFine = typename GET_PROP_TYPE(TypeTag, CprSmootherFine); - using CprSmootherCoarse = typename GET_PROP_TYPE(TypeTag, CprSmootherCoarse); + 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 POrComm = Dune::Amg::SequentialInformation; - using MatrixAdapter = Dune::MatrixAdapter; - -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - -#else - static constexpr int category = Dune::SolverCategory::sequential; - typedef Dune::ScalarProductChooser ScalarProductChooser; -#endif - using CouplingMetric = Opm::Amg::Element; - using CritBase = Dune::Amg::SymmetricCriterion; - using Criterion = Dune::Amg::CoarsenCriterion; - using BlackoilAmgType = BlackoilAmgCpr; + using OperatorParallel = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true>; public: static void registerParameters() @@ -88,12 +86,17 @@ namespace Opm /// \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) + : 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())); @@ -103,44 +106,68 @@ namespace Opm SuperClass::rhs_ = &b; SuperClass::scaleSystem(); const WellModel& wellModel = this->simulator_.problem().wellModel(); - -#if HAVE_MPI - if( this->isParallel() ) { -#if 0 - // Copied from superclass. - typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true> Operator; - auto ebosJacIgnoreOverlap = Matrix(*(this->matrix_)); - //remove ghost rows in local matrix - this->makeOverlapRowsInvalid(ebosJacIgnoreOverlap); - - //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables - //to be certain that correct matrix is used for preconditioning. - Operator opA(ebosJacIgnoreOverlap, ebosJacIgnoreOverlap, wellModel, - this->parallelInformation_ ); - assert( opA.comm() ); - //SuperClass::solve( opA, x, *(this->rhs_), *(opA.comm()) ); - Dune::OwnerOverlapCopyCommunication* comm = opA.comm(); - const size_t size = opA.getmat().N(); - const ParallelISTLInformation& info = +#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. + opAParallel_.reset(new OperatorParallel(*(this->matrix_), *(this->matrix_), wellModel, + this->parallelInformation_ )); + } + + typedef OperatorParallel LinearOperator; + using POrComm = Dune::OwnerOverlapCopyCommunication; + POrComm* comm = opAParallel_->comm(); + + if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) { + 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); - // Construct operator, scalar product and vectors needed. - Dune::InverseOperatorResult result; - SuperClass::constructPreconditionerAndSolve(opA, x, *(this->rhs_), *comm, result); - SuperClass::checkConvergence(result); -#endif // if 0 + // 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); + } + +#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(info)); + 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< LinearOperator, 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 +#endif { - const WellModel& wellModel = this->simulator_.problem().wellModel(); - opASerial_.reset(new OperatorSerial(*(this->matrix_), *(this->matrix_), wellModel)); - typedef Dune::Amg::SequentialInformation POrComm; + + 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; @@ -154,104 +181,109 @@ namespace Opm typedef std::unique_ptr SPPointer; SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg)); sp_ = std::move(sp); -#endif - Vector& istlb = *(this->rhs_); - parallelInformation_arg.copyOwnerToAll(istlb, istlb); +#endif - if( ! std::is_same< LinearOperator, MatrixAdapter > :: value ) { + // 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 ) ); } - const double relax = this->parameters_.ilu_relaxation_; - const MILU_VARIANT ilu_milu = this->parameters_.ilu_milu_; - using Matrix = typename MatrixAdapter::matrix_type; - POrComm& comm = parallelInformation_arg; - 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 - typedef typename BlackoilAmgType::Smoother Smoother; - typedef typename BlackoilAmgType::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); - - MatrixAdapter& opARef = *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 BlackoilAmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) ); - } else { - if (this->parameters_.cpr_solver_verbose_) { - std::cout << " Only update amg solver " << std::endl; - } - amg_->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_; - } - - LinearOperator& opASerialRef = *opASerial_; - linsolve_.reset(new Dune::BiCGSTABSolver(opASerialRef, *sp_, *amg_, - this->parameters_.linear_solver_reduction_, - this->parameters_.linear_solver_maxiter_, - verbosity_linsolve)); + 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 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, *sp_, *amg_, + this->parameters_.linear_solver_reduction_, + this->parameters_.linear_solver_maxiter_, + verbosity_linsolve)); + } + bool solve(Vector& x) { - if (this->isParallel()) { - // for now only call the superclass - bool converged = SuperClass::solve(x); - return converged; - } else { - // 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); - } + // 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_; } @@ -259,21 +291,23 @@ namespace Opm protected: -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) typedef std::shared_ptr< Dune::ScalarProduct > SPPointer; -#else +#else typedef std::unique_ptr SPPointer; #endif - std::unique_ptr< MatrixAdapter > opA_; + ///! \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_; - std::unique_ptr< BlackoilAmgType > amg_; + ///! \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_; + SPPointer sp_; std::shared_ptr< Dune::BiCGSTABSolver > linsolve_; - //std::shared_ptr< Dune::LinearOperator > op_; - //std::shared_ptr< Dune::Preconditioner > prec_; - //std::shared_ptr< Dune::ScalarProduct > sp_; - //Vector solution_; - + const void* oldMat; }; // end ISTLSolver } // namespace Opm From 8057c2ba832c7f0b26d32c5a75b7226b7d9c61ec Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 9 Apr 2019 21:46:05 +0200 Subject: [PATCH 22/25] Use a shared pointer to the OwnerOverlapCommunication in WellModelMatrixAdapter. Previously the OwnerOverlapCommunication object was owned by the adapter as a unique_ptr was used. The change to shared_ptr makes it possible to to extract it and reuse even if OwnerOverlapCommunication is thrown away. --- opm/autodiff/ISTLSolverEbos.hpp | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/ISTLSolverEbos.hpp b/opm/autodiff/ISTLSolverEbos.hpp index e329d1f34..68b80ee08 100644 --- a/opm/autodiff/ISTLSolverEbos.hpp +++ b/opm/autodiff/ISTLSolverEbos.hpp @@ -130,6 +130,14 @@ 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 override { @@ -160,16 +168,16 @@ public: 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_; + std::shared_ptr< communication_type > comm_; }; /// This class solves the fully implicit black-oil system by From 29e893dc8e7f152a6c15b7290374b92e70d004f7 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 9 Apr 2019 21:50:56 +0200 Subject: [PATCH 23/25] Make sure the the OwnerOverlopCommunication object stays alive. As the AMG hierarchy is kept and uses the initial OwnerOverlopCommunication object we need to make sure that this object lives as long as the hierarchy. We assure this by storing a shared_ptr to it in ISTLSolverEbosCpr and reuse it to construct the WellModelMatrixAdapter. --- opm/autodiff/ISTLSolverEbosCpr.hpp | 48 +++++++++++++++++------------- 1 file changed, 27 insertions(+), 21 deletions(-) diff --git a/opm/autodiff/ISTLSolverEbosCpr.hpp b/opm/autodiff/ISTLSolverEbosCpr.hpp index b6e8e766e..f4680a959 100644 --- a/opm/autodiff/ISTLSolverEbosCpr.hpp +++ b/opm/autodiff/ISTLSolverEbosCpr.hpp @@ -116,30 +116,34 @@ namespace Opm 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. - opAParallel_.reset(new OperatorParallel(*(this->matrix_), *(this->matrix_), wellModel, - this->parallelInformation_ )); + 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_ )); + } } - typedef OperatorParallel LinearOperator; using POrComm = Dune::OwnerOverlapCopyCommunication; - POrComm* comm = opAParallel_->comm(); - - if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) { - 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); - } #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::overlapping; - auto sp = Dune::createScalarProduct(*comm, category); + auto sp = Dune::createScalarProduct(*comm_, category); sp_ = std::move(sp); #else constexpr int category = Dune::SolverCategory::overlapping; @@ -151,13 +155,13 @@ namespace Opm using AMGOperator = Dune::OverlappingSchwarzOperator; // If clause is always execute as as Linearoperator is WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false|true>; - if( ! std::is_same< LinearOperator, AMGOperator > :: value && + 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 )); + opA_.reset( new AMGOperator( opAParallel_->getmat(), *comm_ )); } - prepareSolver(*opAParallel_, *comm); + prepareSolver(*opAParallel_, *comm_); } else #endif @@ -308,6 +312,8 @@ namespace Opm SPPointer sp_; std::shared_ptr< Dune::BiCGSTABSolver > linsolve_; const void* oldMat; + using POrComm = Dune::OwnerOverlapCopyCommunication; + std::shared_ptr comm_; }; // end ISTLSolver } // namespace Opm From 81b7af1d8af678c6f8b009772f837234129f9fc9 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 10 Apr 2019 14:40:11 +0200 Subject: [PATCH 24/25] Support DUNE 2.5 with CPR reusing hierarchy --- opm/autodiff/ISTLSolverEbosCpr.hpp | 19 ++++++------- opm/autodiff/amgcpr.hh | 44 ++++++++++++++++++++++++++---- 2 files changed, 47 insertions(+), 16 deletions(-) diff --git a/opm/autodiff/ISTLSolverEbosCpr.hpp b/opm/autodiff/ISTLSolverEbosCpr.hpp index f4680a959..3d07f3300 100644 --- a/opm/autodiff/ISTLSolverEbosCpr.hpp +++ b/opm/autodiff/ISTLSolverEbosCpr.hpp @@ -149,7 +149,7 @@ namespace Opm constexpr int category = Dune::SolverCategory::overlapping; typedef Dune::ScalarProductChooser ScalarProductChooser; typedef std::unique_ptr SPPointer; - SPPointer sp(ScalarProductChooser::construct(info)); + SPPointer sp(ScalarProductChooser::construct(*comm_)); sp_ = std::move(sp); #endif @@ -227,6 +227,9 @@ namespace Opm // 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; @@ -273,7 +276,7 @@ namespace Opm verbosity_linsolve = this->parameters_.linear_solver_verbosity_; } - linsolve_.reset(new Dune::BiCGSTABSolver(wellOpA, *sp_, *amg_, + linsolve_.reset(new Dune::BiCGSTABSolver(wellOpA, reinterpret_cast(*sp_), reinterpret_cast(*amg_), this->parameters_.linear_solver_reduction_, this->parameters_.linear_solver_maxiter_, verbosity_linsolve)); @@ -295,11 +298,6 @@ namespace Opm protected: -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6) - typedef std::shared_ptr< Dune::ScalarProduct > SPPointer; -#else - typedef std::unique_ptr SPPointer; -#endif ///! \brief The dune-istl operator (either serial or parallel std::unique_ptr< Dune::LinearOperator > opA_; ///! \brief Serial well matrix adapter @@ -309,11 +307,12 @@ namespace Opm ///! \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_; + const void* oldMat; + using POrComm = Dune::OwnerOverlapCopyCommunication; + std::shared_ptr comm_; }; // end ISTLSolver } // namespace Opm diff --git a/opm/autodiff/amgcpr.hh b/opm/autodiff/amgcpr.hh index f5ae6aa98..838b220b1 100644 --- a/opm/autodiff/amgcpr.hh +++ b/opm/autodiff/amgcpr.hh @@ -137,11 +137,19 @@ namespace Dune /** \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); @@ -299,8 +307,10 @@ namespace Dune 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_; }; @@ -315,7 +325,9 @@ namespace Dune 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_) @@ -337,8 +349,10 @@ namespace Dune 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()); @@ -360,14 +374,15 @@ namespace Dune 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!"); - // TODO: reestablish compile time checks. - //static_assert(static_cast(PI::category)==static_cast(S::category), - // "Matrix and Solver must match in terms of category!"); +#endif createHierarchies(criterion, const_cast(matrix), pinfo); } @@ -458,7 +473,16 @@ namespace Dune } 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; @@ -498,15 +522,23 @@ namespace Dune // 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()), - *scalarProduct_, + // 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()), - *scalarProduct_, - *coarseSmoother_, 1E-2, 1000, 0)); + // 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; From 9f49e3fc26c5bef357c54d151a6c247fc9535cc7 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 10 Apr 2019 15:01:52 +0200 Subject: [PATCH 25/25] Also pass the parseContext and errorGuard to Opm::checkDeck() This is now mandated since the checkDeck interface was changed recently. Otherwise compilation errors out. --- flow/flow_tag.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flow/flow_tag.hpp b/flow/flow_tag.hpp index c86c30fa5..d879fe82f 100644 --- a/flow/flow_tag.hpp +++ b/flow/flow_tag.hpp @@ -188,7 +188,7 @@ int mainFlow(int argc, char** argv) std::shared_ptr deck = std::make_shared< Opm::Deck >( parser.parseFile(deckFilename , parseContext, errorGuard) ); if ( outputCout ) { - Opm::checkDeck(*deck, parser); + Opm::checkDeck(*deck, parser, parseContext, errorGuard); Opm::MissingFeatures::checkKeywords(*deck); } //Opm::Runspec runspec( *deck );