-- add true impes

-- moved code to hpp to avoid unnesseary logic

opm/simulators/linalg/ISTLSolverEbos.cpp

@@ -81,43 +81,6 @@ void makeOverlapRowsInvalid(Matrix& matrix,
     }
 }
 
-/// Return an appropriate weight function if a cpr preconditioner is asked for.
-template<class Vector, class Matrix>
-std::function<Vector()> getWeightsCalculator(const PropertyTree& prm,
-                                             const Matrix& matrix,
-                                             std::size_t pressureIndex,
-                                             std::function<Vector()> trueFunc)
-{
-    std::function<Vector()> weightsCalculator;
-
-    using namespace std::string_literals;
-
-    auto preconditionerType = prm.get("preconditioner.type"s, "cpr"s);
-    if (preconditionerType == "cpr" || preconditionerType == "cprt"
-        || preconditionerType == "cprw" || preconditionerType == "cprwt") {
-        const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
-        const auto weightsType = prm.get("preconditioner.weight_type"s, "quasiimpes"s);
-        if (weightsType == "quasiimpes") {
-            // weights will be created as default in the solver
-            // assignment p = pressureIndex prevent compiler warning about
-            // capturing variable with non-automatic storage duration
-            weightsCalculator = [matrix, transpose, pressureIndex]() {
-                return Amg::getQuasiImpesWeights<Matrix, Vector>(matrix,
-                                                                 pressureIndex,
-                                                                 transpose);
-            };
-        } else if (weightsType == "trueimpes") {
-            weightsCalculator = trueFunc;
-        } else {
-            OPM_THROW(std::invalid_argument,
-                      "Weights type " + weightsType +
-                      "not implemented for cpr."
-                      " Please use quasiimpes or trueimpes.");
-        }
-    }
-    return weightsCalculator;
-}
-
 template<class Matrix, class Vector, class Comm>
 void FlexibleSolverInfo<Matrix,Vector,Comm>::create(const Matrix& matrix,
                                                     bool parallel,
@@ -152,8 +115,7 @@ void FlexibleSolverInfo<Matrix,Vector,Comm>::create(const Matrix& matrix,
         }
     }
 
-    std::function<Vector()> weightsCalculator =
+    std::function<Vector()> weightsCalculator = trueFunc;
-        getWeightsCalculator<Vector>(prm, matrix, pressureIndex, trueFunc);
 
     if (parallel) {
 #if HAVE_MPI

opm/simulators/linalg/ISTLSolverEbos.hpp

@@ -459,16 +459,11 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
         {
             OPM_TIMEBLOCK(flexibleSolverPrepare);
             if (shouldCreateSolver()) {
-                std::function<Vector()> trueFunc =
-                    [this]
-                    {
-                        return this->getTrueImpesWeights(pressureIndex);
-                    };
-
                 if (!useWellConn_) {
                     auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
                     flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
                 }
+                std::function<Vector()> trueFunc = this->getWeightsCalculator(prm_[activeSolverNum_], getMatrix(), pressureIndex);
                 OPM_TIMEBLOCK(flexibleSolverCreate);
                 flexibleSolver_[activeSolverNum_].create(getMatrix(),
                                                          isParallel(),
@@ -539,18 +534,63 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
         // Weights to make approximate pressure equations.
         // Calculated from the storage terms (only) of the
         // conservation equations, ignoring all other terms.
-        Vector getTrueImpesWeights(int pressureVarIndex) const
+        std::function<Vector()> getWeightsCalculator(const PropertyTree& prm,
+                                                     const Matrix& matrix,
+                                                     std::size_t pressureIndex) const
         {
-            OPM_TIMEBLOCK(getTrueImpesWeights);
+            std::function<Vector()> weightsCalculator;
-            Vector weights(rhs_->size());
+
-            ElementContext elemCtx(simulator_);
+            using namespace std::string_literals;
-            Amg::getTrueImpesWeights(pressureVarIndex, weights,
+
-                                     simulator_.vanguard().gridView(),
+            auto preconditionerType = prm.get("preconditioner.type"s, "cpr"s);
-                                     elemCtx, simulator_.model(),
+            if (preconditionerType == "cpr" || preconditionerType == "cprt"
-                                     ThreadManager::threadId());
+                || preconditionerType == "cprw" || preconditionerType == "cprwt") {
-            return weights;
+                const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
+                const auto weightsType = prm.get("preconditioner.weight_type"s, "quasiimpes"s);
+                if (weightsType == "quasiimpes") {
+                    // weights will be created as default in the solver
+                    // assignment p = pressureIndex prevent compiler warning about
+                    // capturing variable with non-automatic storage duration
+                    weightsCalculator = [matrix, transpose, pressureIndex]() {
+                        return Amg::getQuasiImpesWeights<Matrix, Vector>(matrix,
+                                                                         pressureIndex,
+                                                                         transpose);
+                    };
+                } else if ( weightsType == "trueimpes" ) {
+                    weightsCalculator =
+                        [this, pressureIndex]
+                        {
+                            Vector weights(rhs_->size());
+                            ElementContext elemCtx(simulator_);
+                            Amg::getTrueImpesWeights(pressureIndex, weights,
+                                                             simulator_.vanguard().gridView(),
+                                                             elemCtx, simulator_.model(),
+                                                             ThreadManager::threadId());
+                            return weights;
+                        };
+                } else if  (weightsType == "trueimpesanalytic" ) {
+                    weightsCalculator =
+                        [this, pressureIndex]
+                        {
+                            Vector weights(rhs_->size());
+                            ElementContext elemCtx(simulator_);
+                            Amg::getTrueImpesWeightsAnalytic(pressureIndex, weights,
+                                                             simulator_.vanguard().gridView(),
+                                                             elemCtx, simulator_.model(),
+                                                             ThreadManager::threadId());
+                            return weights;
+                        };
+                } else {
+                    OPM_THROW(std::invalid_argument,
+                              "Weights type " + weightsType +
+                              "not implemented for cpr."
+                              " Please use quasiimpes or trueimpes.");
+                }
+            }
+            return weightsCalculator;
         }
 
+
         Matrix& getMatrix()
         {
             return *matrix_;

opm/simulators/linalg/getQuasiImpesWeights.hpp

@@ -23,7 +23,7 @@
 #include <dune/common/fvector.hh>
 
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
-
+#include <opm/material/common/MathToolbox.hpp>
 #include <algorithm>
 #include <cmath>
 
@@ -132,6 +132,112 @@ namespace Amg
         }
         OPM_END_PARALLEL_TRY_CATCH("getTrueImpesWeights() failed: ", elemCtx.simulator().vanguard().grid().comm());
     }
+
+    template <class Vector, class GridView, class ElementContext, class Model>
+    void getTrueImpesWeightsAnalytic(int pressureVarIndex,
+                                     Vector& weights,
+                                     const GridView& gridView,
+                                     ElementContext& elemCtx,
+                                     const Model& model,
+                                     std::size_t threadId)
+    {
+        // The sequential residual is a linear combination of the
+        // mass balance residuals, with coefficients equal to (for
+        // water, oil, gas):
+        //    1/bw,
+        //    (1/bo - rs/bg)/(1-rs*rv)
+        //    (1/bg - rv/bo)/(1-rs*rv)
+        // These coefficients must be applied for both the residual and
+        // Jacobian.
+        using FluidSystem = typename Model::FluidSystem;
+        using LhsEval = double;
+        using Indices = typename Model::Indices;
+
+        using PrimaryVariables = typename Model::PrimaryVariables;
+        using VectorBlockType = typename Vector::block_type;
+        constexpr int numEq = VectorBlockType::size();
+        using Evaluation =
+            typename std::decay_t<decltype(model.localLinearizer(threadId).localResidual().residual(0))>::block_type;
+        using Toolbox = MathToolbox<Evaluation>;
+        VectorBlockType rhs(0.0);
+        const auto& solution = model.solution(/*timeIdx*/ 0);
+        OPM_BEGIN_PARALLEL_TRY_CATCH();
+         for (const auto& elem : elements(gridView)) {
+            elemCtx.updatePrimaryStencil(elem);
+            elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+            const auto& index = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
+            const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
+            const auto& fs = intQuants.fluidState();
+            VectorBlockType bweights;
+
+            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(
+                    FluidSystem::solventComponentIndex(FluidSystem::waterPhaseIdx));
+                bweights[FluidSystem::waterPhaseIdx]
+                    = Toolbox::template decay<LhsEval>(1 / fs.invB(FluidSystem::waterPhaseIdx));
+            }
+
+            double denominator = 1.0;
+            double rs = Toolbox::template decay<double>(fs.Rs());
+            double rv = Toolbox::template decay<double>(fs.Rv());
+            const auto& priVars = solution[index];
+            if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
+                rs = 0.0;
+            }
+            if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
+                rv = 0.0;
+            }
+            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
+                && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                denominator = Toolbox::template decay<LhsEval>(1 - rs * rv);
+            }
+
+            if (not(denominator > 0)) {
+                std::cout << "Probably negative compressibility" << std::endl;
+            }
+            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(
+                    FluidSystem::solventComponentIndex(FluidSystem::oilPhaseIdx));
+                bweights[activeCompIdx] = Toolbox::template decay<LhsEval>(
+                    (1 / fs.invB(FluidSystem::oilPhaseIdx) - rs / fs.invB(FluidSystem::gasPhaseIdx))
+                    / denominator);
+            }
+            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(
+                    FluidSystem::solventComponentIndex(FluidSystem::gasPhaseIdx));
+                bweights[activeCompIdx] = Toolbox::template decay<LhsEval>(
+                    (1 / fs.invB(FluidSystem::gasPhaseIdx) - rv / fs.invB(FluidSystem::oilPhaseIdx))
+                    / denominator);
+            }
+            if(false){
+                // usure about the best for undersaturated
+
+                if (
+                    (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv)
+                    ||
+                    (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs)
+                    ) {
+                    // Probably not need bweiths may be initialized to zero anyway
+                    unsigned activeCompIdx = -10;
+                    if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
+                        // only water and gas pressent
+                        activeCompIdx = Indices::canonicalToActiveComponentIndex(
+                        FluidSystem::solventComponentIndex(FluidSystem::oilPhaseIdx));
+                        bweights[activeCompIdx] = 0.0;
+                                                                                         }
+                    if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
+                        // only water and oil pressent
+                        activeCompIdx = Indices::canonicalToActiveComponentIndex(
+                            FluidSystem::solventComponentIndex(FluidSystem::gasPhaseIdx));
+                        bweights[activeCompIdx] = 0.0;
+                }
+
+            }
+            }
+            weights[index] = bweights;
+        }
+        OPM_END_PARALLEL_TRY_CATCH("getTrueImpesAnalyticWeights() failed: ", elemCtx.simulator().vanguard().grid().comm());
+    }
 } // namespace Amg
 
 } // namespace Opm