Make shouldCreateSolver() and getWeightsCalculator() separate functions.

OPM · Aug 3, 2020 · 3a52a1d · 3a52a1d
1 parent 2337395
commit 3a52a1d
Showing 1 changed file with 59 additions and 52 deletions.
diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp
@@ -330,9 +330,46 @@ namespace Opm
         }
 
         void prepareFlexibleSolver()
+        {
+
+            std::function<Vector()> weightsCalculator = getWeightsCalculator();
+
+            if (shouldCreateSolver()) {
+                if (isParallel()) {
+#if HAVE_MPI
+                    if (useWellConn_) {
+                        using ParOperatorType = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
+                        linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
+                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator);
+                    } else {
+                        using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
+                        wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
+                        linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
+                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator);
+                    }
+#endif
+                } else {
+                    using SeqLinearOperator = Dune::MatrixAdapter<Matrix, Vector, Vector>;
+                    linearOperatorForFlexibleSolver_ = std::make_unique<SeqLinearOperator>(getMatrix());
+                    flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator);
+                }
+            }
+            else
+            {
+                flexibleSolver_->preconditioner().update();
+            }
+        }
+
+
+        /// Return true if we should (re)create the whole solver,
+        /// instead of just calling update() on the preconditioner.
+        bool shouldCreateSolver() const
         {
             // Decide if we should recreate the solver or just do
             // a minimal preconditioner update.
+            if (!flexibleSolver_) {
+                return true;
+            }
             const int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
             bool recreate_solver = false;
             if (this->parameters_.cpr_reuse_setup_ == 0) {
@@ -353,66 +390,36 @@ namespace Opm
                 assert(recreate_solver == false);
                 // Never recreate solver.
             }
+            return recreate_solver;
+        }
 
+
+        /// Return an appropriate weight function if a cpr preconditioner is asked for.
+        std::function<Vector()> getWeightsCalculator() const
+        {
             std::function<Vector()> weightsCalculator;
 
             auto preconditionerType = prm_.get("preconditioner.type", "cpr");
-            if( preconditionerType  == "cpr" ||
-                preconditionerType == "cprt"
-                )
-            {
-                bool transpose = false;
-                if(preconditionerType == "cprt"){
-                    transpose = true;
-                }
-
-                auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes");
-                auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1);
-                if(weightsType == "quasiimpes") {
+            if (preconditionerType == "cpr" || preconditionerType == "cprt") {
+                const bool transpose = preconditionerType == "cprt";
+                const auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes");
+                const auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1);
+                if (weightsType == "quasiimpes") {
                     // weighs will be created as default in the solver
-                    weightsCalculator =
-                        [this, transpose, pressureIndex](){
-                            return Opm::Amg::getQuasiImpesWeights<Matrix,
-                                                                  Vector>(getMatrix(),
-                                                                          pressureIndex,
-                                                                          transpose);
-                        };
-
-                }else if(weightsType == "trueimpes"  ){
-                    weightsCalculator =
-                        [this, pressureIndex](){
-                            return this->getTrueImpesWeights(pressureIndex);
-                        };
-                }else{
-                    OPM_THROW(std::invalid_argument, "Weights type " << weightsType << "not implemented for cpr."
-                              << " Please use quasiimpes or trueimpes.");
-                }
-            }
-
-            if (recreate_solver || !flexibleSolver_) {
-                if (isParallel()) {
-#if HAVE_MPI
-                    if (useWellConn_) {
-                        using ParOperatorType = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Comm>;
-                        linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
-                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator);
-                    } else {
-                        using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
-                        wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
-                        linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
-                        flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator);
-                    }
-#endif
+                    weightsCalculator = [this, transpose, pressureIndex]() {
+                        return Opm::Amg::getQuasiImpesWeights<Matrix, Vector>(this->getMatrix(), pressureIndex, transpose);
+                    };
+                } else if (weightsType == "trueimpes") {
+                    weightsCalculator = [this, pressureIndex]() {
+                        return this->getTrueImpesWeights(pressureIndex);
+                    };
                 } else {
-                    using SeqLinearOperator = Dune::MatrixAdapter<Matrix, Vector, Vector>;
-                    linearOperatorForFlexibleSolver_ = std::make_unique<SeqLinearOperator>(getMatrix());
-                    flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator);
+                    OPM_THROW(std::invalid_argument,
+                              "Weights type " << weightsType << "not implemented for cpr."
+                              << " Please use quasiimpes or trueimpes.");
                 }
             }
-            else
-            {
-                flexibleSolver_->preconditioner().update();
-            }
+            return weightsCalculator;
         }