use analytic weights for cprw by default - except when solvent is enabled

opm/simulators/flow/BlackoilModelNldd.hpp

@@ -99,6 +99,9 @@ class BlackoilModelNldd
 
     static constexpr int numEq = Indices::numEq;
 
+    enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
+    static constexpr bool canUseAnalyticWeightsForCprw = !enableSolvent;
+
     //! \brief The constructor sets up the subdomains.
     //! \param model BlackOil model to solve for
     //! \param param param Model parameters
@@ -203,7 +206,7 @@ class BlackoilModelNldd
             const auto& eclState = model_.simulator().vanguard().eclState();
             FlowLinearSolverParameters loc_param;
             loc_param.is_nldd_local_solver_ = true;
-            loc_param.init(eclState.getSimulationConfig().useCPR());
+            loc_param.init(eclState.getSimulationConfig().useCPR(), canUseAnalyticWeightsForCprw);
             // Override solver type with umfpack if small domain.
             if (domains_[index].cells.size() < 200) {
                 loc_param.linsolver_ = "umfpack";

opm/simulators/linalg/FlowLinearSolverParameters.cpp

@@ -29,7 +29,7 @@
 
 namespace Opm {
 
-void FlowLinearSolverParameters::init(bool cprRequestedInDataFile)
+void FlowLinearSolverParameters::init(bool cprRequestedInDataFile, bool cprwUseAnalyticWeights)
 {
     // TODO: these parameters have undocumented non-trivial dependencies
     relaxed_linear_solver_reduction_ = Parameters::Get<Parameters::RelaxedLinearSolverReduction>();
@@ -48,6 +48,7 @@ void FlowLinearSolverParameters::init(bool cprRequestedInDataFile)
     cpr_reuse_setup_  = Parameters::Get<Parameters::CprReuseSetup>();
     cpr_reuse_interval_  = Parameters::Get<Parameters::CprReuseInterval>();
     gpu_aware_mpi_ = Parameters::Get<Parameters::GpuAwareMpi>();
+    cprw_use_analytic_weights_ = cprwUseAnalyticWeights;
 
     if (!Parameters::IsSet<Parameters::LinearSolver>() && cprRequestedInDataFile) {
         linsolver_ = "cpr";

opm/simulators/linalg/FlowLinearSolverParameters.hpp

@@ -122,10 +122,11 @@ struct FlowLinearSolverParameters
     bool gpu_aware_mpi_;
     bool verify_gpu_aware_mpi_;
     bool cpr_weights_thread_parallel_;
+    bool cprw_use_analytic_weights_;
 
     FlowLinearSolverParameters() { reset(); }
 
-    void init(bool cprRequestedInDataFile);
+    void init(bool cprRequestedInDataFile, bool cprwUseAnalyticWeights);
 
     static void registerParameters();
 

opm/simulators/linalg/ISTLSolver.hpp

@@ -166,9 +166,11 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
         using ElementChunksType = ElementChunks<GridView, Dune::Partitions::All>;
 
         constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
-
+        enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>()};
         enum { enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>() };
-        constexpr static bool isIncompatibleWithCprw = enablePolymerMolarWeight;
+
+        static constexpr bool isIncompatibleWithCprw = enablePolymerMolarWeight;
+        static constexpr bool canUseAnalyticWeightsForCprw = !enableSolvent;
 
 #if HAVE_MPI
         using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
@@ -212,7 +214,9 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
               matrix_(nullptr)
         {
             parameters_.resize(1);
-            parameters_[0].init(simulator_.vanguard().eclState().getSimulationConfig().useCPR());
+            parameters_[0].init(simulator_.vanguard().eclState().getSimulationConfig().useCPR(),
+                canUseAnalyticWeightsForCprw);
+
             initialize();
         }
 
@@ -228,7 +232,6 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
                               "Choose a different option, for example --linear-solver=ilu0");
                 }
             }
-
             if (parameters_[0].linsolver_ == "hybrid") {
                 // Experimental hybrid configuration.
                 // When chosen, will set up two solvers, one with CPRW
@@ -238,7 +241,7 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
                 parameters_.clear();
                 {
                     FlowLinearSolverParameters para;
-                    para.init(false);
+                    para.init(false, canUseAnalyticWeightsForCprw);
                     para.linsolver_ = "cprw";
                     parameters_.push_back(para);
                     prm_.push_back(setupPropertyTree(parameters_[0],
@@ -247,7 +250,7 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
                 }
                 {
                     FlowLinearSolverParameters para;
-                    para.init(false);
+                    para.init(false, canUseAnalyticWeightsForCprw);
                     para.linsolver_ = "ilu0";
                     parameters_.push_back(para);
                     prm_.push_back(setupPropertyTree(parameters_[1],

opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp

@@ -77,7 +77,9 @@ class ISTLSolverGPUISTL : public AbstractISTLSolver<TypeTag>
     using GPUVectorInt = Opm::gpuistl::GpuVector<int>;
 
     constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
-
+
+    enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
+    static constexpr bool canUseAnalyticWeightsForCprw = !enableSolvent;
 
 #if HAVE_MPI
     using CommunicationType = Dune::OwnerOverlapCopyCommunication<int, int>;
@@ -119,7 +121,8 @@ class ISTLSolverGPUISTL : public AbstractISTLSolver<TypeTag>
 #else
         m_comm = std::make_shared<CommunicationType>(simulator.gridView().comm());
 #endif
-        m_parameters.init(simulator.vanguard().eclState().getSimulationConfig().useCPR());
+        m_parameters.init(simulator.vanguard().eclState().getSimulationConfig().useCPR(),
+            canUseAnalyticWeightsForCprw);
         m_propertyTree = setupPropertyTree(m_parameters,
                                            Parameters::IsSet<Parameters::LinearSolverMaxIter>(),
                                            Parameters::IsSet<Parameters::LinearSolverReduction>());

opm/simulators/linalg/setupPropertyTree.cpp

@@ -290,7 +290,12 @@ setupCPRW(const std::string& /*conf*/, const FlowLinearSolverParameters& p)
     prm.put("preconditioner.type", "cprw"s);
     prm.put("preconditioner.use_well_weights", "false"s);
     prm.put("preconditioner.add_wells", "true"s);
-    prm.put("preconditioner.weight_type", "trueimpes"s);
+    // Can't use trueimpesanalytic when e.g., solvent is active
+    if (p.cprw_use_analytic_weights_) {
+        prm.put("preconditioner.weight_type", "trueimpesanalytic"s);
+    } else {
+        prm.put("preconditioner.weight_type", "trueimpes"s);
+    }
     prm.put("preconditioner.pre_smooth", 0);
     prm.put("preconditioner.post_smooth", 1);
     prm.put("preconditioner.finesmoother.type", "paroverilu0"s);