From 70cff4b342f8d22dea856cb530c2a4a9adf86f95 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Mon, 23 Mar 2020 22:36:22 +0100 Subject: [PATCH] Refactor to a single function that calculates true impes weights. Two identical versions seemed like a lot of maintenance for no gain. --- opm/simulators/linalg/ISTLSolverEbos.hpp | 36 ++-------------- .../linalg/ISTLSolverEbosFlexible.hpp | 37 ++--------------- .../linalg/getQuasiImpesWeights.hpp | 41 +++++++++++++++++++ 3 files changed, 47 insertions(+), 67 deletions(-) diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 621f3e826..d13a99a28 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -902,40 +902,10 @@ protected: Vector getStorageWeights() const { Vector weights(rhs_->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); - 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; - } + Opm::Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(), + elemCtx, simulator_.model(), + ThreadManager::threadId()); return weights; } diff --git a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp index 82f6c6c68..26ab2b563 100644 --- a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp +++ b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp @@ -254,41 +254,10 @@ 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; - } + Opm::Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(), + elemCtx, simulator_.model(), + ThreadManager::threadId()); return weights; } diff --git a/opm/simulators/linalg/getQuasiImpesWeights.hpp b/opm/simulators/linalg/getQuasiImpesWeights.hpp index e782b909b..c95d0aca9 100644 --- a/opm/simulators/linalg/getQuasiImpesWeights.hpp +++ b/opm/simulators/linalg/getQuasiImpesWeights.hpp @@ -83,6 +83,47 @@ namespace Amg return weights; } + template + void getTrueImpesWeights(int pressureVarIndex, Vector& weights, const GridView& gridView, + ElementContext& elemCtx, const Model& model, std::size_t threadId) + { + using VectorBlockType = typename Vector::block_type; + using Matrix = typename std::decay_t; + using MatrixBlockType = typename Matrix::MatrixBlock; + constexpr int numEq = VectorBlockType::size(); + using Evaluation = typename std::decay_t + ::block_type; + VectorBlockType rhs(0.0); + rhs[pressureVarIndex] = 1.0; + int index = 0; + auto elemIt = gridView.template begin(); + const auto& elemEndIt = gridView.template end(); + for (; elemIt != elemEndIt; ++elemIt) { + elemCtx.updatePrimaryStencil(*elemIt); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + Dune::FieldVector storage; + model.localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,/*spaceIdx=*/0, /*timeIdx=*/0); + auto extrusionFactor = elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor(); + auto scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor; + auto 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; + } + } + } + VectorBlockType bweights; + MatrixBlockType block_transpose = Details::transposeDenseMatrix(block); + block_transpose.solve(bweights, rhs); + bweights /= 1000.0; // given normal densities this scales weights to about 1. + weights[index] = bweights; + ++index; + } + } } // namespace Amg } // namespace Opm