From bfa859c099644278a0e9fd3e4c7a99daf40ba3bf Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 21 Aug 2019 10:50:57 +0200 Subject: [PATCH] changes to make true impes work with flexible solver Cherry-picked commit 7da5ce4fec33. Which got cleaned up (deleted trailing whitespace, tabs) and made compilable with recent changes. --- .../linalg/ISTLSolverEbosFlexible.hpp | 80 ++++++++++++++++++- .../linalg/OwningBlockPreconditioner.hpp | 4 +- .../linalg/OwningTwoLevelPreconditioner.hpp | 30 +++++-- .../linalg/ParallelOverlappingILU0.hpp | 2 +- .../linalg/PreconditionerWithUpdate.hpp | 6 +- .../linalg/PressureSolverPolicy.hpp | 3 +- opm/simulators/linalg/amgcpr.hh | 8 +- 7 files changed, 117 insertions(+), 16 deletions(-) diff --git a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp index 27fde5cd7..93b6abe40 100644 --- a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp +++ b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp @@ -65,6 +65,19 @@ class ISTLSolverEbosFlexible #endif using SolverType = Dune::FlexibleSolver; + // for quasiImpesWeights + 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 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; + public: static void registerParameters() @@ -135,6 +148,30 @@ public: assert(recreate_solver == false); // Never recreate solver. } + if( prm_.get("preconditioner.type") == "cpr" || + prm_.get("preconditioner.type") == "cprt" + ) + { + bool transpose = false; + if(prm_.get("preconditioner.type") == "cprt"){ + transpose = true; + } + + + if(prm_.get("preconditioner.weight_type") == "quasiimpes") { + VectorType weights = Opm::Amg::getQuasiImpesWeights( + mat.istlMatrix(), + prm_.get("preconditioner.pressure_var_index"), transpose); + prm_.put("preconditioner.weights",weights); + + }else if(prm_.get("preconditioner.weight_type") == "trueimpes" ){ + VectorType weights = + this->getTrueImpesWeights(b, prm_.get("preconditioner.pressure_var_index")); + prm_.put("preconditioner.weights",weights); + }else{ + throw std::runtime_error("no such weights implemented for cpr"); + } + } if (recreate_solver || !solver_) { if (isParallel()) { @@ -146,7 +183,7 @@ public: } rhs_ = b; } else { - solver_->preconditioner().update(); + solver_->preconditioner().update(prm_); rhs_ = b; } } @@ -204,6 +241,47 @@ protected: } } + VectorType getTrueImpesWeights(const VectorType& b,const int pressureVarIndex) + { + VectorType weights(b.size()); + BlockVector rhs(0.0); + rhs[pressureVarIndex] = 1.0; + int index = 0; + ElementContext elemCtx(simulator_); + const auto& vanguard = simulator_.vanguard(); + auto elemIt = vanguard.gridView().template begin(); + const auto& elemEndIt = vanguard.gridView().template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const Element& elem = *elemIt; + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + const int numEq = MatrixType::block_type::rows; + Dune::FieldVector storage; + unsigned threadId = ThreadManager::threadId(); + simulator_.model().localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,/*spaceIdx=*/0, /*timeIdx=*/0); + Scalar extrusionFactor = elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor(); + Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor; + Scalar storage_scale = scvVolume / elemCtx.simulator().timeStepSize(); + MatrixBlockType block; + double pressure_scale = 50e5; + for (int ii = 0; ii < numEq; ++ii) { + for (int jj = 0; jj < numEq; ++jj) { + block[ii][jj] = storage[ii].derivative(jj)/storage_scale; + if (jj == pressureVarIndex) { + block[ii][jj] *= pressure_scale; + } + } + } + BlockVector bweights; + MatrixBlockType block_transpose = Opm::transposeDenseMatrix(block); + block_transpose.solve(bweights, rhs); + bweights /= 1000.0; // given normal densities this scales weights to about 1. + weights[index] = bweights; + ++index; + } + return weights; + } + const Simulator& simulator_; diff --git a/opm/simulators/linalg/OwningBlockPreconditioner.hpp b/opm/simulators/linalg/OwningBlockPreconditioner.hpp index 7dcd2d384..5b5860478 100644 --- a/opm/simulators/linalg/OwningBlockPreconditioner.hpp +++ b/opm/simulators/linalg/OwningBlockPreconditioner.hpp @@ -63,9 +63,9 @@ public: } // The update() function does nothing for a wrapped preconditioner. - virtual void update() override + virtual void update(const boost::property_tree::ptree& p) override { - orig_precond_.update(); + orig_precond_.update(p); } private: diff --git a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp index bbba9e21f..08b54412f 100644 --- a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp +++ b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp @@ -53,6 +53,21 @@ namespace Dune template class FlexibleSolver; +template +std::ostream& operator<<(std::ostream& out, + const BlockVector< FieldVector< T, i >, A >& vector) +{ + Dune::writeMatrixMarket(vector, out); + return out; +} + + template + std::istream& operator>>(std::istream& input, + BlockVector< FieldVector< T, i >, A >& vector) +{ + Dune::readMatrixMarket(vector, input); + return input; +} /// A version of the two-level preconditioner that is: /// - Self-contained, because it owns its policy components. @@ -73,8 +88,7 @@ public: : linear_operator_(linearoperator) , finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother"))) , comm_(nullptr) - , weights_(Opm::Amg::getQuasiImpesWeights( - linearoperator.getmat(), prm.get("pressure_var_index"), transpose)) + , weights_(prm.get("weights")) , levelTransferPolicy_(dummy_comm_, weights_, prm.get("pressure_var_index")) , coarseSolverPolicy_(prm.get_child("coarsesolver")) , twolevel_method_(linearoperator, @@ -98,8 +112,7 @@ public: : linear_operator_(linearoperator) , finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother"), comm)) , comm_(&comm) - , weights_(Opm::Amg::getQuasiImpesWeights( - linearoperator.getmat(), prm.get("pressure_var_index"), transpose)) + , weights_(prm.get("weights")) , levelTransferPolicy_(*comm_, weights_, prm.get("pressure_var_index")) , coarseSolverPolicy_(prm.get_child("coarsesolver")) , twolevel_method_(linearoperator, @@ -110,6 +123,8 @@ public: transpose ? 0 : 1) , prm_(prm) { + //Opm::Amg::getQuasiImpesWeights( + // linearoperator.getmat(), prm.get("pressure_var_index"), transpose)) if (prm.get("verbosity") > 10) { std::ofstream outfile(prm.get("weights_filename")); if (!outfile) { @@ -134,10 +149,11 @@ public: twolevel_method_.post(x); } - virtual void update() override + virtual void update(const boost::property_tree::ptree& prm) override { - Opm::Amg::getQuasiImpesWeights( - linear_operator_.getmat(), prm_.get("pressure_var_index"), transpose, weights_); + //Opm::Amg::getQuasiImpesWeights( + // linear_operator_.getmat(), prm_.get("pressure_var_index"), transpose, weights_); + weights_ = prm.get("weights"); updateImpl(comm_); } diff --git a/opm/simulators/linalg/ParallelOverlappingILU0.hpp b/opm/simulators/linalg/ParallelOverlappingILU0.hpp index c5f70e7c1..0f75f2b32 100644 --- a/opm/simulators/linalg/ParallelOverlappingILU0.hpp +++ b/opm/simulators/linalg/ParallelOverlappingILU0.hpp @@ -921,7 +921,7 @@ public: DUNE_UNUSED_PARAMETER(x); } - virtual void update() override + virtual void update(const boost::property_tree::ptree& = boost::property_tree::ptree()) override { // (For older DUNE versions the communicator might be // invalid if redistribution in AMG happened on the coarset level. diff --git a/opm/simulators/linalg/PreconditionerWithUpdate.hpp b/opm/simulators/linalg/PreconditionerWithUpdate.hpp index ba9111dfd..a0d0cede6 100644 --- a/opm/simulators/linalg/PreconditionerWithUpdate.hpp +++ b/opm/simulators/linalg/PreconditionerWithUpdate.hpp @@ -22,7 +22,7 @@ #include #include - +#include namespace Dune { @@ -31,7 +31,7 @@ template class PreconditionerWithUpdate : public Preconditioner { public: - virtual void update() = 0; + virtual void update(const boost::property_tree::ptree& pt) = 0; }; template @@ -69,7 +69,7 @@ public: } // The update() function does nothing for a wrapped preconditioner. - virtual void update() override + virtual void update(const boost::property_tree::ptree& /* pt */) override { } diff --git a/opm/simulators/linalg/PressureSolverPolicy.hpp b/opm/simulators/linalg/PressureSolverPolicy.hpp index 5692031ca..21f4a2937 100644 --- a/opm/simulators/linalg/PressureSolverPolicy.hpp +++ b/opm/simulators/linalg/PressureSolverPolicy.hpp @@ -84,7 +84,8 @@ namespace Amg void updatePreconditioner() { - linsolver_->preconditioner().update(); + pt::ptree prm; + linsolver_->preconditioner().update(prm); } private: diff --git a/opm/simulators/linalg/amgcpr.hh b/opm/simulators/linalg/amgcpr.hh index 301683cdb..4176cdfd7 100644 --- a/opm/simulators/linalg/amgcpr.hh +++ b/opm/simulators/linalg/amgcpr.hh @@ -239,8 +239,9 @@ namespace Dune /** * @brief Update the coarse solver and the hierarchies. */ - virtual void update(); + virtual void update(const boost::property_tree::ptree& prm); + virtual void update(); /** * @brief Check whether the coarse solver used is a direct solver. * @return True if the coarse level solver is a direct solver. @@ -468,6 +469,11 @@ namespace Dune update(); } + template + void AMGCPR::update(const boost::property_tree::ptree& /* prm */) + { + update(); + } template void AMGCPR::update() {