From 9ccb67ead617c5d5753ad3c7c777fa115de97033 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Fri, 23 Aug 2024 13:08:08 +0200
Subject: [PATCH] Intercept CPRW if MICP or polymer injectivity is used.

These modeling options are currently incompatible with our CPRW implementation.
---
 opm/simulators/linalg/ISTLSolver.hpp | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/opm/simulators/linalg/ISTLSolver.hpp b/opm/simulators/linalg/ISTLSolver.hpp
index b04383bd3..80c9a73ee 100644
--- a/opm/simulators/linalg/ISTLSolver.hpp
+++ b/opm/simulators/linalg/ISTLSolver.hpp
@@ -35,6 +35,7 @@
 #include <opm/models/utils/propertysystem.hh>
 #include <opm/simulators/flow/BlackoilModelParameters.hpp>
 #include <opm/simulators/flow/FlowBaseVanguard.hpp>
+#include <opm/simulators/flow/FlowProblemProperties.hpp>
 #include <opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp>
 #include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
 #include <opm/simulators/linalg/matrixblock.hh>
@@ -159,6 +160,10 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
         using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
         constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
 
+        enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
+        enum { enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>() };
+        constexpr static bool isIncompatibleWithCprw = enableMICP || enablePolymerMolarWeight;
+
 #if HAVE_MPI
         using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
 #else
@@ -211,6 +216,21 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
         {
             OPM_TIMEBLOCK(IstlSolver);
 
+            if (isIncompatibleWithCprw) {
+                // Some model variants are incompatible with the CPRW linear solver.
+                if (parameters_[0].linsolver_ == "cprw" || parameters_[0].linsolver_ == "hybrid") {
+                    std::string incompatible_model = "Unknown";
+                    if (enableMICP) {
+                        incompatible_model = "MICP";
+                    } else if (enablePolymerMolarWeight) {
+                        incompatible_model = "Polymer injectivity";
+                    }
+                    OPM_THROW(std::runtime_error,
+                              incompatible_model + " model is incompatible with the CPRW linear solver.\n"
+                              "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