diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake
index 38c9711f0..663f0570e 100755
--- a/compareECLFiles.cmake
+++ b/compareECLFiles.cmake
@@ -363,21 +363,19 @@ add_test_compareECLFiles(CASENAME spe1_metric_vfp1
                          REL_TOL ${rel_tol}
                          DIR vfpprod_spe1)
 
-if(BUILD_FLOW_VARIANTS)
-  add_test_compareECLFiles(CASENAME spe1_water
-                           FILENAME SPE1CASE1_WATER
-                           SIMULATOR flow_onephase
-                           ABS_TOL ${abs_tol}
-                           REL_TOL ${rel_tol}
-                           DIR spe1)
+add_test_compareECLFiles(CASENAME spe1_water
+                         FILENAME SPE1CASE1_WATER
+                         SIMULATOR flow
+                         ABS_TOL ${abs_tol}
+                         REL_TOL ${rel_tol}
+                         DIR spe1)
 
-  add_test_compareECLFiles(CASENAME spe1_thermal_onephase
-                           FILENAME SPE1CASE2_THERMAL_ONEPHASE
-                           SIMULATOR flow_onephase_energy
-                           ABS_TOL ${abs_tol}
-                           REL_TOL ${rel_tol}
-                           DIR spe1)
-endif()
+add_test_compareECLFiles(CASENAME spe1_thermal_onephase
+                         FILENAME SPE1CASE2_THERMAL_ONEPHASE
+                         SIMULATOR flow
+                         ABS_TOL ${abs_tol}
+                         REL_TOL ${rel_tol}
+                         DIR spe1)
 
 add_test_compareECLFiles(CASENAME spe1_spider
                            FILENAME SPIDER_CAKESLICE
@@ -1431,10 +1429,9 @@ if(MPI_FOUND)
                                        DIR spe1
                                        TEST_ARGS --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-8)
 
-if(BUILD_FLOW_VARIANTS)
   add_test_compare_parallel_simulation(CASENAME spe1_thermal
                                        FILENAME SPE1CASE2_THERMAL_ONEPHASE
-                                       SIMULATOR flow_onephase_energy
+                                       SIMULATOR flow
                                        ABS_TOL ${abs_tol}
                                        REL_TOL ${rel_tol}
                                        DIR spe1
@@ -1442,12 +1439,11 @@ if(BUILD_FLOW_VARIANTS)
 
   add_test_compare_parallel_simulation(CASENAME spe1_water
                                        FILENAME SPE1CASE1_WATER
-                                       SIMULATOR flow_onephase
+                                       SIMULATOR flow
                                        ABS_TOL ${abs_tol}
                                        REL_TOL ${rel_tol}
                                        DIR spe1
                                        TEST_ARGS --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-8)
-endif()
 
   add_test_compare_parallel_simulation(CASENAME spe1_brine
                                        FILENAME SPE1CASE1_BRINE
diff --git a/ebos/eclgenericoutputblackoilmodule.cc b/ebos/eclgenericoutputblackoilmodule.cc
index 6cb6a74ec..d95a2c237 100644
--- a/ebos/eclgenericoutputblackoilmodule.cc
+++ b/ebos/eclgenericoutputblackoilmodule.cc
@@ -681,6 +681,8 @@ assignToSolution(data::Solution& sol)
         {"PPCW",     UnitSystem::measure::pressure,  data::TargetType::RESTART_SOLUTION,      ppcw_},
         {"PRESROCC", UnitSystem::measure::pressure,  data::TargetType::RESTART_SOLUTION,      minimumOilPressure_},
         {"PRESSURE", UnitSystem::measure::pressure,  data::TargetType::RESTART_SOLUTION,      oilPressure_},
+        {"PCOW",     UnitSystem::measure::pressure,  data::TargetType::RESTART_SOLUTION,      pcow_},
+        {"PCOG",     UnitSystem::measure::pressure,  data::TargetType::RESTART_SOLUTION,      pcog_},
         {"PRES_OVB", UnitSystem::measure::pressure,  data::TargetType::RESTART_SOLUTION,      overburdenPressure_},
         {"RS",       UnitSystem::measure::gas_oil_ratio, data::TargetType::RESTART_SOLUTION,  rs_},
         {"RSSAT",    UnitSystem::measure::gas_oil_ratio, data::TargetType::RESTART_AUXILIARY, gasDissolutionFactor_},
@@ -1094,6 +1096,15 @@ doAllocBuffers(unsigned bufferSize,
         relativePermeability_[gasPhaseIdx].resize(bufferSize, 0.0);
     }
 
+    if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(waterPhaseIdx) && rstKeywords["PCOW"] > 0) {
+        rstKeywords["PCOW"] = 0;
+        pcow_.resize(bufferSize, 0.0);
+    }
+    if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx) && rstKeywords["PCOG"] > 0) {
+        rstKeywords["PCOG"] = 0;
+        pcog_.resize(bufferSize, 0.0);
+    }
+
     if (rstKeywords["PBPD"] > 0)  {
         rstKeywords["PBPD"] = 0;
         bubblePointPressure_.resize(bufferSize, 0.0);
diff --git a/ebos/eclgenericoutputblackoilmodule.hh b/ebos/eclgenericoutputblackoilmodule.hh
index 5a4f94b05..6727c026f 100644
--- a/ebos/eclgenericoutputblackoilmodule.hh
+++ b/ebos/eclgenericoutputblackoilmodule.hh
@@ -450,6 +450,8 @@ protected:
     ScalarBuffer cUrea_;
     ScalarBuffer cBiofilm_;
     ScalarBuffer cCalcite_;
+    ScalarBuffer pcow_;
+    ScalarBuffer pcog_;
 
     std::array<ScalarBuffer, numPhases> saturation_;
     std::array<ScalarBuffer, numPhases> invB_;
diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh
index a606c13ca..d866d0802 100644
--- a/ebos/ecloutputblackoilmodule.hh
+++ b/ebos/ecloutputblackoilmodule.hh
@@ -281,6 +281,14 @@ public:
                 this->rv_[globalDofIdx] = getValue(fs.Rv());
                 Valgrind::CheckDefined(this->rv_[globalDofIdx]);
             }
+            if (!this->pcow_.empty()) {
+                this->pcow_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
+                Valgrind::CheckDefined(this->pcow_[globalDofIdx]);
+            }
+            if (!this->pcog_.empty()) {
+                this->pcog_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx));
+                Valgrind::CheckDefined(this->pcog_[globalDofIdx]);
+            }
 
             if (!this->rvw_.empty()) {
                 this->rvw_[globalDofIdx] = getValue(fs.Rvw());
@@ -571,6 +579,14 @@ public:
                         val.second = getValue(fs.viscosity(gasPhaseIdx));
                     else if (key.first == "BVOIL" || key.first == "BOVIS")
                         val.second = getValue(fs.viscosity(oilPhaseIdx));
+                    else if (key.first == "BRPV")
+                        val.second = elemCtx.simulator().model().dofTotalVolume(globalDofIdx)*getValue(intQuants.porosity());
+                    else if (key.first == "BOPV")
+                        val.second = getValue(fs.saturation(oilPhaseIdx))*elemCtx.simulator().model().dofTotalVolume(globalDofIdx)*getValue(intQuants.porosity());
+                    else if (key.first == "BWPV")
+                        val.second = getValue(fs.saturation(waterPhaseIdx))*elemCtx.simulator().model().dofTotalVolume(globalDofIdx)*getValue(intQuants.porosity());
+                    else if (key.first == "BGPV")
+                        val.second = getValue(fs.saturation(gasPhaseIdx))*elemCtx.simulator().model().dofTotalVolume(globalDofIdx)*getValue(intQuants.porosity());
                     else {
                         std::string logstring = "Keyword '";
                         logstring.append(key.first);
diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh
index 8f784683b..08dc8f213 100644
--- a/ebos/eclproblem.hh
+++ b/ebos/eclproblem.hh
@@ -2119,8 +2119,9 @@ private:
                 Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(),t,p);
                 Scalar saturatedInvB = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p);
                 Scalar rsZero = 0.0;
-                Scalar pureInvB = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p,rsZero);
-                Scalar deltaDensity = (saturatedInvB - pureInvB) * FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex());
+                Scalar pureDensity = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p,rsZero) * FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex());
+                Scalar saturatedDensity = saturatedInvB * (FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex()) + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, fs.pvtRegionIndex()));
+                Scalar deltaDensity = saturatedDensity - pureDensity;
                 Scalar rs = getValue(fs.Rs());
                 Scalar visc = FluidSystem::oilPvt().viscosity(fs.pvtRegionIndex(),t,p,rs);
                 Scalar poro =  getValue(iq.porosity());
diff --git a/ebos/equil/initstateequil.hh b/ebos/equil/initstateequil.hh
index d2ed3b212..c7e4937fb 100644
--- a/ebos/equil/initstateequil.hh
+++ b/ebos/equil/initstateequil.hh
@@ -219,9 +219,12 @@ private:
     density(const double depth,
             const double press) const
     {
-        double rs = rs_(depth, press, temp_);
+        double rs = 0.0;
+        if(FluidSystem::enableDissolvedGas())
+            rs = rs_(depth, press, temp_);
+
         double bOil = 0.0;
-        if (!FluidSystem::enableDissolvedGas() || rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp_, press)) {
+        if (rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp_, press)) {
             bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
         }
         else {
@@ -267,9 +270,12 @@ private:
     density(const double depth,
             const double press) const
     {
-        double rv = rv_(depth, press, temp_);
+        double rv = 0.0;
+        if (FluidSystem::enableVaporizedOil())
+            rv = rv_(depth, press, temp_);
+
         double bGas = 0.0;
-        if (!FluidSystem::enableVaporizedOil() || rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press)) {
+        if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press)) {
             bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
         }
         else {
diff --git a/flow/flow_ebos_onephase.cpp b/flow/flow_ebos_onephase.cpp
index 1da67791a..f777f0af4 100644
--- a/flow/flow_ebos_onephase.cpp
+++ b/flow/flow_ebos_onephase.cpp
@@ -27,13 +27,13 @@ namespace Opm {
 namespace Properties {
 
 namespace TTag {
-struct EclFlowProblemSimple {
+struct EclFlowProblemWaterOnly {
     using InheritsFrom = std::tuple<EclFlowProblem>;
 };
 }
 //! The indices required by the model
 template<class TypeTag>
-struct Indices<TypeTag, TTag::EclFlowProblemSimple>
+struct Indices<TypeTag, TTag::EclFlowProblemWaterOnly>
 {
 private:
     // it is unfortunately not possible to simply use 'TypeTag' here because this leads
@@ -57,9 +57,36 @@ public:
 
 } // namespace Opm::Properties
 
-int flowEbosOnephaseMain(int argc, char** argv)
+void flowEbosWaterOnlySetDeck(double setupTime, std::shared_ptr<Deck> deck,
+                              std::shared_ptr<EclipseState> eclState,
+                              std::shared_ptr<Schedule> schedule,
+                              std::shared_ptr<SummaryConfig> summaryConfig)
 {
-    using TypeTag = Opm::Properties::TTag::EclFlowProblemSimple;
+    using TypeTag = Properties::TTag::EclFlowProblemWaterOnly;
+    using Vanguard = GetPropType<TypeTag, Properties::Vanguard>;
+
+    Vanguard::setExternalSetupTime(setupTime);
+    Vanguard::setExternalDeck(std::move(deck));
+    Vanguard::setExternalEclState(std::move(eclState));
+    Vanguard::setExternalSchedule(std::move(schedule));
+    Vanguard::setExternalSummaryConfig(std::move(summaryConfig));
+}
+
+// ----------------- Main program -----------------
+int flowEbosWaterOnlyMain(int argc, char** argv, bool outputCout, bool outputFiles)
+{
+    // we always want to use the default locale, and thus spare us the trouble
+    // with incorrect locale settings.
+    resetLocale();
+
+    FlowMainEbos<Properties::TTag::EclFlowProblemWaterOnly>
+        mainfunc {argc, argv, outputCout, outputFiles};
+    return mainfunc.execute();
+}
+
+int flowEbosWaterOnlyMainStandalone(int argc, char** argv)
+{
+    using TypeTag = Opm::Properties::TTag::EclFlowProblemWaterOnly;
     auto mainObject = Opm::Main(argc, argv);
     return mainObject.runStatic<TypeTag>();
 }
diff --git a/flow/flow_ebos_onephase.hpp b/flow/flow_ebos_onephase.hpp
index 4992de344..9bd45262d 100644
--- a/flow/flow_ebos_onephase.hpp
+++ b/flow/flow_ebos_onephase.hpp
@@ -22,8 +22,27 @@
 #ifndef FLOW_ONEPHASE_HPP
 #define FLOW_ONEPHASE_HPP
 
+#include <memory>
+
 namespace Opm {
-  int flowEbosOnephaseMain(int argc, char** argv);
+
+class Deck;
+class EclipseState;
+template<class TypeTag> class FlowMainEbos;
+class Schedule;
+class SummaryConfig;
+class UDQState;
+class WellTestState;
+namespace Action {
+class State;
+}
+
+void flowEbosWaterOnlySetDeck(double setupTime, std::shared_ptr<Deck> deck,
+                              std::shared_ptr<EclipseState> eclState,
+                              std::shared_ptr<Schedule> schedule,
+                              std::shared_ptr<SummaryConfig> summaryConfig);
+int flowEbosWaterOnlyMain(int argc, char** argv, bool outputCout, bool outputFiles);
+int flowEbosWaterOnlyMainStandalone(int argc, char** argv);
 }
 
 #endif
diff --git a/flow/flow_ebos_onephase_energy.cpp b/flow/flow_ebos_onephase_energy.cpp
index 505aa9f27..19ce0df67 100644
--- a/flow/flow_ebos_onephase_energy.cpp
+++ b/flow/flow_ebos_onephase_energy.cpp
@@ -19,24 +19,28 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 #include "config.h"
+
+#include <flow/flow_ebos_onephase_energy.hpp>
+
 #include <opm/simulators/flow/Main.hpp>
 #include <opm/models/blackoil/blackoilonephaseindices.hh>
+#include <opm/material/common/ResetLocale.hpp>
 
 namespace Opm {
 namespace Properties {
 
 namespace TTag {
-struct EclFlowProblemSimple {
+struct EclFlowProblemWaterOnlyEnergy {
     using InheritsFrom = std::tuple<EclFlowProblem>;
 };
 }
 template<class TypeTag>
-struct EnableEnergy<TypeTag, TTag::EclFlowProblemSimple> {
+struct EnableEnergy<TypeTag, TTag::EclFlowProblemWaterOnlyEnergy> {
     static constexpr bool value = true;
 };
 //! The indices required by the model
 template<class TypeTag>
-struct Indices<TypeTag, TTag::EclFlowProblemSimple>
+struct Indices<TypeTag, TTag::EclFlowProblemWaterOnlyEnergy>
 {
 private:
     // it is unfortunately not possible to simply use 'TypeTag' here because this leads
@@ -60,9 +64,36 @@ public:
 
 } // namespace Opm::Properties
 
-int flowEbosOnephaseEnergyMain(int argc, char** argv)
+void flowEbosWaterOnlyEnergySetDeck(double setupTime, std::shared_ptr<Deck> deck,
+                                    std::shared_ptr<EclipseState> eclState,
+                                    std::shared_ptr<Schedule> schedule,
+                                    std::shared_ptr<SummaryConfig> summaryConfig)
 {
-    using TypeTag = Opm::Properties::TTag::EclFlowProblemSimple;
+    using TypeTag = Properties::TTag::EclFlowProblemWaterOnlyEnergy;
+    using Vanguard = GetPropType<TypeTag, Properties::Vanguard>;
+
+    Vanguard::setExternalSetupTime(setupTime);
+    Vanguard::setExternalDeck(std::move(deck));
+    Vanguard::setExternalEclState(std::move(eclState));
+    Vanguard::setExternalSchedule(std::move(schedule));
+    Vanguard::setExternalSummaryConfig(std::move(summaryConfig));
+}
+
+// ----------------- Main program -----------------
+int flowEbosWaterOnlyEnergyMain(int argc, char** argv, bool outputCout, bool outputFiles)
+{
+    // we always want to use the default locale, and thus spare us the trouble
+    // with incorrect locale settings.
+    resetLocale();
+
+    FlowMainEbos<Properties::TTag::EclFlowProblemWaterOnlyEnergy>
+        mainfunc {argc, argv, outputCout, outputFiles};
+    return mainfunc.execute();
+}
+
+int flowEbosWaterOnlyEnergyMainStandalone(int argc, char** argv)
+{
+    using TypeTag = Opm::Properties::TTag::EclFlowProblemWaterOnlyEnergy;
     auto mainObject = Opm::Main(argc, argv);
     return mainObject.runStatic<TypeTag>();
 }
diff --git a/flow/flow_ebos_onephase_energy.hpp b/flow/flow_ebos_onephase_energy.hpp
index b70577a29..657bf529b 100644
--- a/flow/flow_ebos_onephase_energy.hpp
+++ b/flow/flow_ebos_onephase_energy.hpp
@@ -22,8 +22,27 @@
 #ifndef FLOW_ONEPHASE_ENERGY_HPP
 #define FLOW_ONEPHASE_ENERGY_HPP
 
+#include <memory>
+
 namespace Opm {
-  int flowEbosOnephaseEnergyMain(int argc, char** argv);
+
+class Deck;
+class EclipseState;
+template<class TypeTag> class FlowMainEbos;
+class Schedule;
+class SummaryConfig;
+class UDQState;
+class WellTestState;
+namespace Action {
+class State;
+}
+
+void flowEbosWaterOnlyEnergySetDeck(double setupTime, std::shared_ptr<Deck> deck,
+                                    std::shared_ptr<EclipseState> eclState,
+                                    std::shared_ptr<Schedule> schedule,
+                                    std::shared_ptr<SummaryConfig> summaryConfig);
+int flowEbosWaterOnlyEnergyMain(int argc, char** argv, bool outputCout, bool outputFiles);
+int flowEbosWaterOnlyEnergyMainStandalone(int argc, char** argv);
 }
 
 #endif
diff --git a/flow/flow_onephase.cpp b/flow/flow_onephase.cpp
index a398d4683..b1133df38 100644
--- a/flow/flow_onephase.cpp
+++ b/flow/flow_onephase.cpp
@@ -24,5 +24,5 @@
 
 int main(int argc, char** argv)
 {
-    return Opm::flowEbosOnephaseMain(argc, argv);
+    return Opm::flowEbosWaterOnlyMainStandalone(argc, argv);
 }
diff --git a/flow/flow_onephase_energy.cpp b/flow/flow_onephase_energy.cpp
index 8601c5cf5..6c3e436b1 100644
--- a/flow/flow_onephase_energy.cpp
+++ b/flow/flow_onephase_energy.cpp
@@ -23,5 +23,5 @@
 
 int main(int argc, char** argv)
 {
-    return Opm::flowEbosOnephaseEnergyMain(argc,argv);
+    return Opm::flowEbosWaterOnlyEnergyMainStandalone(argc,argv);
 }
diff --git a/opm/core/props/satfunc/RelpermDiagnostics.cpp b/opm/core/props/satfunc/RelpermDiagnostics.cpp
index 7fdfa5f5f..a410686ea 100644
--- a/opm/core/props/satfunc/RelpermDiagnostics.cpp
+++ b/opm/core/props/satfunc/RelpermDiagnostics.cpp
@@ -112,13 +112,16 @@ namespace Opm{
         const TableContainer& sof2Tables = tableManager.getSof2Tables();
         const TableContainer& sgwfnTables= tableManager.getSgwfnTables();
 
+        const SwofletTable& swofletTable = tableManager.getSwofletTable();
+        const SgofletTable& sgofletTable = tableManager.getSgofletTable();
+
         // Family I test.
         bool family1 = pu.phase_used[BlackoilPhases::Liquid];
         if (pu.phase_used[BlackoilPhases::Aqua]) {
-            family1 = family1 && !swofTables.empty();
+            family1 = family1 && (!swofTables.empty() || !swofletTable.empty());
         }
         if (pu.phase_used[BlackoilPhases::Vapour]) {
-            family1 = family1 && (!sgofTables.empty() || !slgofTables.empty());
+            family1 = family1 && ((!sgofTables.empty() || !sgofletTable.empty()) || !slgofTables.empty());
         }
 
         // Family II test.
@@ -664,7 +667,9 @@ namespace Opm{
             satfunc::getRawFunctionValues(tables, phases, rtep);
 
         const TableContainer&  swofTables = tables.getSwofTables();
+        const SwofletTable&  swofletTables = tables.getSwofletTable();
         const TableContainer&  sgofTables = tables.getSgofTables();
+        const SgofletTable&  sgofletTables = tables.getSgofletTable();
         const TableContainer& slgofTables = tables.getSlgofTables();
         const TableContainer&  sof3Tables = tables.getSof3Tables();
 
@@ -693,14 +698,19 @@ namespace Opm{
                      if (!sgofTables.empty()) {
                          const auto& table = sgofTables.getTable<SgofTable>(satnumIdx);
                          krog_value = table.evaluate( "KROG" , unscaledEpsInfo_[satnumIdx].Sgl );
+                     } else if (!sgofletTables.empty()) {
+                         krog_value = sgofletTables[satnumIdx].krt2_relperm;
                      } else {
                          assert(!slgofTables.empty());
                          const auto& table = slgofTables.getTable<SlgofTable>(satnumIdx);
                          krog_value = table.evaluate( "KROG" , unscaledEpsInfo_[satnumIdx].Sgl );
                      }
-                     {
+                     if (!swofTables.empty()) {
                          const auto& table = swofTables.getTable<SwofTable>(satnumIdx);
                          krow_value = table.evaluate("KROW" , unscaledEpsInfo_[satnumIdx].Swl);
+                     } else {
+                         assert(!swofletTables.empty());
+                         krow_value = swofletTables[satnumIdx].krt2_relperm;
                      }
                  }
                  if (satFamily_ == SaturationFunctionFamily::FamilyII) {
diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp
index 4a8631e2d..3fe292105 100644
--- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp
+++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp
@@ -267,7 +267,7 @@ struct MaxInnerIterWells<TypeTag, TTag::FlowModelParameters> {
 };
 template<class TypeTag>
 struct ShutUnsolvableWells<TypeTag, TTag::FlowModelParameters> {
-    static constexpr bool value = false;
+    static constexpr bool value = true;
 };
 template<class TypeTag>
 struct AlternativeWellRateInit<TypeTag, TTag::FlowModelParameters> {
@@ -275,7 +275,7 @@ struct AlternativeWellRateInit<TypeTag, TTag::FlowModelParameters> {
 };
 template<class TypeTag>
 struct StrictOuterIterWells<TypeTag, TTag::FlowModelParameters> {
-    static constexpr int value = 99;
+    static constexpr int value = 6;
 };
 template<class TypeTag>
 struct StrictInnerIterWells<TypeTag, TTag::FlowModelParameters> {
@@ -297,12 +297,12 @@ struct EnableWellOperabilityCheckIter<TypeTag, TTag::FlowModelParameters> {
 template<class TypeTag>
 struct RelaxedWellFlowTol<TypeTag, TTag::FlowModelParameters> {
     using type = GetPropType<TypeTag, Scalar>;
-    static constexpr type value = 1;
+    static constexpr type value = 1e-3;
 };
 template<class TypeTag>
 struct RelaxedPressureTolMsw<TypeTag, TTag::FlowModelParameters> {
     using type = GetPropType<TypeTag, Scalar>;
-    static constexpr type value = 0.5e5;
+    static constexpr type value = 1.0e4;
 };
 template<class TypeTag>
 struct MaximumNumberOfWellSwitches<TypeTag, TTag::FlowModelParameters> {
diff --git a/opm/simulators/flow/Main.hpp b/opm/simulators/flow/Main.hpp
index 2762a6bce..8ce99f990 100644
--- a/opm/simulators/flow/Main.hpp
+++ b/opm/simulators/flow/Main.hpp
@@ -36,6 +36,8 @@
 #include <flow/flow_ebos_brine_saltprecipitation.hpp>
 #include <flow/flow_ebos_gaswater_saltprec_vapwat.hpp>
 #include <flow/flow_ebos_brine_precsalt_vapwat.hpp>
+#include <flow/flow_ebos_onephase.hpp>
+#include <flow/flow_ebos_onephase_energy.hpp>
 #include <flow/flow_ebos_oilwater_brine.hpp>
 #include <flow/flow_ebos_gaswater_brine.hpp>
 #include <flow/flow_ebos_energy.hpp>
@@ -291,6 +293,16 @@ private:
             return this->runMICP(phases);
         }
 
+        // water-only case
+        else if(phases.size() == 1 && phases.active(Phase::WATER) && !eclipseState_->getSimulationConfig().isThermal()) {
+            return this->runWaterOnly(phases);
+        }
+
+        // water-only case with energy
+        else if(phases.size() == 2 && phases.active(Phase::WATER) && eclipseState_->getSimulationConfig().isThermal()) {
+            return this->runWaterOnlyEnergy(phases);
+        }
+
         // Twophase cases
         else if (phases.size() == 2 && !eclipseState_->getSimulationConfig().isThermal()) {
             return this->runTwoPhase(phases);
@@ -653,6 +665,36 @@ private:
         return flowEbosFoamMain(argc_, argv_, outputCout_, outputFiles_);
     }
 
+    int runWaterOnly(const Phases& phases)
+    {
+        if (!phases.active(Phase::WATER) || phases.size() != 1) {
+            if (outputCout_)
+                std::cerr << "No valid configuration is found for water-only simulation, valid options include "
+                          << "water, water + thermal" << std::endl;
+
+            return EXIT_FAILURE;
+        }
+        flowEbosWaterOnlySetDeck(
+            setupTime_, deck_, eclipseState_, schedule_, summaryConfig_);
+
+        return flowEbosWaterOnlyMain(argc_, argv_, outputCout_, outputFiles_);
+    }
+    
+    int runWaterOnlyEnergy(const Phases& phases)
+    {
+        if (!phases.active(Phase::WATER) || phases.size() != 2) {
+            if (outputCout_)
+                std::cerr << "No valid configuration is found for water-only simulation, valid options include "
+                          << "water, water + thermal" << std::endl;
+
+            return EXIT_FAILURE;
+        }
+        flowEbosWaterOnlyEnergySetDeck(
+            setupTime_, deck_, eclipseState_, schedule_, summaryConfig_);
+
+        return flowEbosWaterOnlyEnergyMain(argc_, argv_, outputCout_, outputFiles_);
+    }
+
     int runBrine(const Phases& phases)
     {
         if (! phases.active(Phase::WATER) || phases.size() == 2) {
diff --git a/opm/simulators/timestepping/ConvergenceReport.hpp b/opm/simulators/timestepping/ConvergenceReport.hpp
index 273d6bda9..7884aa345 100644
--- a/opm/simulators/timestepping/ConvergenceReport.hpp
+++ b/opm/simulators/timestepping/ConvergenceReport.hpp
@@ -64,7 +64,7 @@ namespace Opm
         class WellFailure
         {
         public:
-            enum struct Type { Invalid, MassBalance, Pressure, ControlBHP, ControlTHP, ControlRate };
+            enum struct Type { Invalid, MassBalance, Pressure, ControlBHP, ControlTHP, ControlRate, Unsolvable };
             WellFailure(Type t, Severity s, int phase, const std::string& well_name)
                 : type_(t), severity_(s), phase_(phase), well_name_(well_name)
             {
diff --git a/opm/simulators/utils/UnsupportedFlowKeywords.cpp b/opm/simulators/utils/UnsupportedFlowKeywords.cpp
index 8197a8010..44eddd2ea 100755
--- a/opm/simulators/utils/UnsupportedFlowKeywords.cpp
+++ b/opm/simulators/utils/UnsupportedFlowKeywords.cpp
@@ -140,8 +140,10 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"DPGRID", {false, std::nullopt}},
         {"DPKRMOD", {false, std::nullopt}},
         {"DPNUM", {false, std::nullopt}},
+        {"DR", {false, std::string{"Use the DRV keyword instead"}}},
         {"DRILPRI", {false, std::nullopt}},
         {"DSPDEINT", {false, std::nullopt}},
+        {"DTHETA", {false, std::string{"Use the DTHETAV keyword instead"}}},
         {"DUALPERM", {false, std::nullopt}},
         {"DUALPORO", {false, std::nullopt}},
         {"DUMPCUPL", {false, std::nullopt}},
@@ -430,6 +432,7 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"OLDTRAN", {false, std::nullopt}},
         {"OLDTRANR", {false, std::nullopt}},
         {"OPTIONS", {false, std::nullopt}},
+        {"OUTRAD", {false, std::string{"Use the DRV keyword instead"}}},
         {"OUTSOL", {false, std::nullopt}},
         {"PARAOPTS", {false, std::nullopt}},
         {"PCG32D", {false, std::nullopt}},
diff --git a/opm/simulators/utils/readDeck.cpp b/opm/simulators/utils/readDeck.cpp
index 85fc1c83a..72c9dd6af 100644
--- a/opm/simulators/utils/readDeck.cpp
+++ b/opm/simulators/utils/readDeck.cpp
@@ -365,7 +365,7 @@ namespace {
 
         if (status != MPI_SUCCESS) {
             throw std::invalid_argument {
-                "Unable to establish cell geomtry validity across MPI ranks"
+                "Unable to establish cell geometry validity across MPI ranks"
             };
         }
 #endif  // HAVE_MPI
diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
index bd3b5e610..011547eb4 100644
--- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp
+++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
@@ -2067,7 +2067,7 @@ forceShutWellByName(const std::string& wellname,
     // process that owns it.
     int well_was_shut = 0;
     for (const auto& well : well_container_generic_) {
-        if (well->name() == wellname && !well->wellIsStopped()) {
+        if (well->name() == wellname) {
             wellTestState().close_well(wellname, WellTestConfig::Reason::PHYSICAL, simulation_time);
             well_was_shut = 1;
             break;
diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp
index d6f1d21a7..3b26d5bba 100644
--- a/opm/simulators/wells/BlackoilWellModel_impl.hpp
+++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp
@@ -669,7 +669,6 @@ namespace Opm {
                     case Well::ProducerCMode::RESV:
                         zero_rate_control = is_zero(prod_controls.resv_rate);
                         break;
-
                     default:
                         // Might still have zero rate controls, but is pressure controlled.
                         zero_rate_control = false;
@@ -1223,8 +1222,13 @@ namespace Opm {
         ConvergenceReport local_report;
         const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
         for (const auto& well : well_container_) {
-            if (well->isOperableAndSolvable() ) {
+            if (well->isOperableAndSolvable() || well->wellIsStopped()) {
                 local_report += well->getWellConvergence(this->wellState(), B_avg, local_deferredLogger, iterationIdx > param_.strict_outer_iter_wells_ );
+            } else {
+                ConvergenceReport report;
+                using CR = ConvergenceReport;
+                report.setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()});
+                local_report += report;
             }
         }
         
@@ -1488,6 +1492,8 @@ namespace Opm {
             auto& events = this->wellState().well(well->indexOfWell()).events;
             if (events.hasEvent(WellState::event_mask)) {
                 well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), deferred_logger);
+                well->updatePrimaryVariables(this->wellState(), deferred_logger);
+                well->initPrimaryVariablesEvaluation();
                 // There is no new well control change input within a report step,
                 // so next time step, the well does not consider to have effective events anymore.
                 events.clearEvent(WellState::event_mask);
diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp
index d54198448..edb5a54ea 100644
--- a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp
+++ b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp
@@ -294,19 +294,38 @@ checkThpControl_() const
     return thp_control;
 }
 
+std::pair<std::optional<double>,double>
+GasLiftSingleWellGeneric::
+computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_() const
+{
+    auto alq = this->orig_alq_;
+    double new_alq = alq;
+    std::optional<double> bhp;
+    while (alq <= (this->max_alq_ + this->increment_)) {
+        if (bhp = computeBhpAtThpLimit_(alq); bhp) {
+            new_alq = alq;
+            break;
+        }
+        alq += this->increment_;
+    }
+    return {bhp, new_alq};
+}
 
-std::optional<GasLiftSingleWellGeneric::BasicRates>
+std::pair<std::optional<GasLiftSingleWellGeneric::BasicRates>,double>
 GasLiftSingleWellGeneric::
 computeInitialWellRates_() const
 {
     std::optional<BasicRates> rates;
-    if (auto bhp = computeBhpAtThpLimit_(this->orig_alq_); bhp) {
+    double initial_alq = this->orig_alq_;
+    //auto alq = initial_alq;
+    //if (auto bhp = computeBhpAtThpLimit_(this->orig_alq_); bhp) {
+    if (auto [bhp, alq] = computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_(); bhp) {
         {
             const std::string msg = fmt::format(
-                "computed initial bhp {} given thp limit and given alq {}",
-                *bhp, this->orig_alq_);
+                "computed initial bhp {} given thp limit and given alq {}", *bhp, alq);
             displayDebugMessage_(msg);
         }
+        initial_alq = alq;
         auto [new_bhp, bhp_is_limited] = getBhpWithLimit_(*bhp);
         rates = computeWellRates_(new_bhp, bhp_is_limited);
         if (rates) {
@@ -320,7 +339,7 @@ computeInitialWellRates_() const
     else {
         displayDebugMessage_("Aborting optimization.");
     }
-    return rates;
+    return {rates, initial_alq};
 }
 
 std::optional<GasLiftSingleWellGeneric::LimitedRates>
@@ -819,12 +838,13 @@ getRateWithGroupLimit_(
 }
 
 
-std::optional<GasLiftSingleWellGeneric::LimitedRates>
+std::pair<std::optional<GasLiftSingleWellGeneric::LimitedRates>, double>
 GasLiftSingleWellGeneric::
 getInitialRatesWithLimit_() const
 {
     std::optional<LimitedRates> limited_rates;
-    if (auto rates = computeInitialWellRates_(); rates) {
+    double initial_alq = this->orig_alq_;
+    if (auto [rates, alq] = computeInitialWellRates_(); rates) {
         if (this->debug) {
             displayDebugMessage_(
                              "Maybe limiting initial rates before optimize loop..");
@@ -832,8 +852,9 @@ getInitialRatesWithLimit_() const
         auto temp_rates = getLimitedRatesFromRates_(*rates);
         BasicRates old_rates = getWellStateRates_();
         limited_rates = updateRatesToGroupLimits_(old_rates, temp_rates);
+        initial_alq = alq;
     }
-    return limited_rates;
+    return {limited_rates, initial_alq};
 }
 
 GasLiftSingleWellGeneric::LimitedRates
@@ -1166,18 +1187,17 @@ runOptimizeLoop_(bool increase)
 {
     if (this->debug) debugShowProducerControlMode();
     std::unique_ptr<GasLiftWellState> ret_value; // nullptr initially
-    auto rates = getInitialRatesWithLimit_();
+    auto [rates, cur_alq] = getInitialRatesWithLimit_();
     if (!rates) return ret_value;
     // if (this->debug) debugShowBhpAlqTable_();
     if (this->debug) debugShowAlqIncreaseDecreaseCounts_();
     if (this->debug) debugShowTargets_();
     bool success = false;  // did we succeed to increase alq?
     bool alq_is_limited = false;
-    auto cur_alq = this->orig_alq_;
     LimitedRates new_rates = *rates;
     auto [temp_rates2, new_alq] = maybeAdjustALQbeforeOptimizeLoop_(
                                                  *rates, cur_alq, increase);
-    if (checkInitialALQmodified_(new_alq, cur_alq)) {
+    if (checkInitialALQmodified_(new_alq, this->orig_alq_)) {
         auto delta_alq = new_alq - cur_alq;
         new_rates = temp_rates2;
         cur_alq = new_alq;
@@ -1591,7 +1611,10 @@ checkAlqOutsideLimits(double alq, [[maybe_unused]] double oil_rate)
                 // NOTE: checking for an upper limit should not be necessary
                 // when decreasing alq.. so this is just to catch an
                 // illegal state at an early point.
-                if (alq >= this->parent.max_alq_) {
+                if (this->parent.checkALQequal_(alq, this->parent.max_alq_)) {
+                    return false;
+                }
+                else if (alq > this->parent.max_alq_) {
                     warn_( "unexpected: alq above upper limit when trying to "
                         "decrease lift gas. aborting iteration.");
                     result = true;
diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.hpp b/opm/simulators/wells/GasLiftSingleWellGeneric.hpp
index 79f1a376c..93120cd6f 100644
--- a/opm/simulators/wells/GasLiftSingleWellGeneric.hpp
+++ b/opm/simulators/wells/GasLiftSingleWellGeneric.hpp
@@ -250,7 +250,8 @@ protected:
     bool checkInitialALQmodified_(double alq, double initial_alq) const;
     bool checkThpControl_() const;
     virtual std::optional<double> computeBhpAtThpLimit_(double alq) const = 0;
-    std::optional<BasicRates> computeInitialWellRates_() const;
+    std::pair<std::optional<double>,double> computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_() const;
+    std::pair<std::optional<BasicRates>,double> computeInitialWellRates_() const;
     std::optional<LimitedRates> computeLimitedWellRatesWithALQ_(double alq) const;
     virtual BasicRates computeWellRates_(double bhp, bool bhp_is_limited,                                                             bool debug_output = true) const = 0;
     std::optional<BasicRates> computeWellRatesWithALQ_(double alq) const;
@@ -272,7 +273,7 @@ protected:
                            const BasicRates& rates) const;
     std::pair<double, bool> getGasRateWithGroupLimit_(
                            double new_gas_rate, double gas_rate) const;
-    std::optional<LimitedRates> getInitialRatesWithLimit_() const;
+    std::pair<std::optional<LimitedRates>,double> getInitialRatesWithLimit_() const;
     LimitedRates getLimitedRatesFromRates_(const BasicRates& rates) const;
     std::tuple<double,double,bool,bool> getLiquidRateWithGroupLimit_(
                            const double new_oil_rate, const double oil_rate,
diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp
index 097cf8115..413f7f522 100644
--- a/opm/simulators/wells/MultisegmentWell.hpp
+++ b/opm/simulators/wells/MultisegmentWell.hpp
@@ -162,6 +162,9 @@ namespace Opm
     protected:
         int number_segments_;
 
+        // regularize msw equation
+        bool regularize_;
+
         // components of the pressure drop to be included
         WellSegments::CompPressureDrop compPressureDrop() const;
         // multi-phase flow model
diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp
index 0e034b806..3f1e79f7a 100644
--- a/opm/simulators/wells/MultisegmentWellEval.cpp
+++ b/opm/simulators/wells/MultisegmentWellEval.cpp
@@ -383,10 +383,10 @@ processFractions(const int seg) const
 template<typename FluidSystem, typename Indices, typename Scalar>
 void
 MultisegmentWellEval<FluidSystem,Indices,Scalar>::
-updateWellState(const BVectorWell& dwells,
-                const double relaxation_factor,
-                const double dFLimit,
-                const double max_pressure_change) const
+updatePrimaryVariablesNewton(const BVectorWell& dwells,
+                             const double relaxation_factor,
+                             const double dFLimit,
+                             const double max_pressure_change) const
 {
     const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_;
 
@@ -1622,10 +1622,10 @@ assemblePressureEq(const int seg,
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
-std::vector<Scalar>
+std::pair<bool, std::vector<Scalar> >
 MultisegmentWellEval<FluidSystem,Indices,Scalar>::
-getWellResiduals(const std::vector<Scalar>& B_avg,
-                 DeferredLogger& deferred_logger) const
+getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
+                       DeferredLogger& deferred_logger) const
 {
     assert(int(B_avg.size() ) == baseif_.numComponents());
     std::vector<Scalar> residuals(numWellEq + 1, 0.0);
@@ -1641,8 +1641,9 @@ getWellResiduals(const std::vector<Scalar>& B_avg,
                 }
             }
             if (std::isnan(residual) || std::isinf(residual)) {
-                OPM_DEFLOG_THROW(NumericalIssue, "nan or inf value for residal get for well " << baseif_.name()
-                                                << " segment " << seg << " eq_idx " << eq_idx, deferred_logger);
+                deferred_logger.debug("nan or inf value for residal get for well " + baseif_.name()
+                                      + " segment " + std::to_string(seg) + " eq_idx " + std::to_string(eq_idx));
+                return {false, residuals};
             }
 
             if (residual > residuals[eq_idx]) {
@@ -1655,12 +1656,13 @@ getWellResiduals(const std::vector<Scalar>& B_avg,
     {
         const double control_residual = std::abs(resWell_[0][numWellEq - 1]);
         if (std::isnan(control_residual) || std::isinf(control_residual)) {
-           OPM_DEFLOG_THROW(NumericalIssue, "nan or inf value for control residal get for well " << baseif_.name(), deferred_logger);
+           deferred_logger.debug("nan or inf value for control residal get for well " + baseif_.name());
+           return {false, residuals};
         }
         residuals[numWellEq] = control_residual;
     }
 
-    return residuals;
+    return {true, residuals};
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
diff --git a/opm/simulators/wells/MultisegmentWellEval.hpp b/opm/simulators/wells/MultisegmentWellEval.hpp
index 0654e1e77..48a444550 100644
--- a/opm/simulators/wells/MultisegmentWellEval.hpp
+++ b/opm/simulators/wells/MultisegmentWellEval.hpp
@@ -173,10 +173,10 @@ protected:
     void updateUpwindingSegments();
 
     // updating the well_state based on well solution dwells
-    void updateWellState(const BVectorWell& dwells,
-                         const double relaxation_factor,
-                         const double DFLimit,
-                         const double max_pressure_change) const;
+    void updatePrimaryVariablesNewton(const BVectorWell& dwells,
+                                      const double relaxation_factor,
+                                      const double DFLimit,
+                                      const double max_pressure_change) const;
 
     void computeSegmentFluidProperties(const EvalWell& temperature,
                                        const EvalWell& saltConcentration,
@@ -199,8 +199,9 @@ protected:
     EvalWell getWQTotal() const;
 
 
-    std::vector<Scalar> getWellResiduals(const std::vector<Scalar>& B_avg,
-                                         DeferredLogger& deferred_logger) const;
+    std::pair<bool, std::vector<Scalar> >
+    getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
+                           DeferredLogger& deferred_logger) const;
 
     double getControlTolerance(const WellState& well_state,
                                const double tolerance_wells,
diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp
index e75fba6c5..72e27d279 100644
--- a/opm/simulators/wells/MultisegmentWell_impl.hpp
+++ b/opm/simulators/wells/MultisegmentWell_impl.hpp
@@ -49,6 +49,7 @@ namespace Opm
                      const std::vector<PerforationData>& perf_data)
     : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
     , MSWEval(static_cast<WellInterfaceIndices<FluidSystem,Indices,Scalar>&>(*this))
+    , regularize_(false)
     , segment_fluid_initial_(this->numberOfSegments(), std::vector<double>(this->num_components_, 0.0))
     {
         // not handling solvent or polymer for now with multisegment well
@@ -316,11 +317,11 @@ namespace Opm
             well_potentials[phase] *= sign;
             total_potential += well_potentials[phase];
         }
-        if (total_potential < 0.0) {
+        if (total_potential < 0.0 && this->param_.check_well_operability_) {
             // wells with negative potentials are not operable
             this->operability_status_.has_negative_potentials = true;
             const std::string msg = std::string("well ") + this->name() + std::string(": has non negative potentials is not operable");
-            deferred_logger.info(msg);
+            deferred_logger.warning("NEGATIVE_POTENTIALS_INOPERABLE", msg);
         }
     }
 
@@ -611,10 +612,10 @@ namespace Opm
 
         const double dFLimit = this->param_.dwell_fraction_max_;
         const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
-        this->MSWEval::updateWellState(dwells,
-                                       relaxation_factor,
-                                       dFLimit,
-                                       max_pressure_change);
+        this->MSWEval::updatePrimaryVariablesNewton(dwells,
+                                                    relaxation_factor,
+                                                    dFLimit,
+                                                    max_pressure_change);
 
         this->updateWellStateFromPrimaryVariables(well_state, getRefDensity(), deferred_logger);
         Base::calculateReservoirRates(well_state.well(this->index_of_well_));
@@ -1373,7 +1374,14 @@ namespace Opm
 
         const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
         const WellState well_state0 = well_state;
-        const std::vector<Scalar> residuals0 = this->getWellResiduals(Base::B_avg_, deferred_logger);
+
+        {
+            // getWellFiniteResiduals returns false for nan/inf residuals
+            const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
+            if(!isFinite)
+                return false;
+        }
+
         std::vector<std::vector<Scalar> > residual_history;
         std::vector<double> measure_history;
         int it = 0;
@@ -1383,14 +1391,17 @@ namespace Opm
         bool converged = false;
         int stagnate_count = 0;
         bool relax_convergence = false;
+        this->regularize_ = false;
         for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
 
             assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
 
             const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
 
-            if (it > this->param_.strict_inner_iter_wells_)
+            if (it > this->param_.strict_inner_iter_wells_) {
                 relax_convergence = true;
+                this->regularize_ = true;
+            }
 
             const auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger, relax_convergence);
             if (report.converged()) {
@@ -1398,12 +1409,20 @@ namespace Opm
                 break;
             }
 
-            residual_history.push_back(this->getWellResiduals(Base::B_avg_, deferred_logger));
-            measure_history.push_back(this->getResidualMeasureValue(well_state,
+            {
+                // getFinteWellResiduals returns false for nan/inf residuals
+                const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
+                if (!isFinite)
+                    return false;
+
+                residual_history.push_back(residuals);
+                measure_history.push_back(this->getResidualMeasureValue(well_state,
                                                                     residual_history[it],
                                                                     this->param_.tolerance_wells_,
                                                                     this->param_.tolerance_pressure_ms_wells_,
                                                                     deferred_logger) );
+            }
+
 
             bool is_oscillate = false;
             bool is_stagnate = false;
@@ -1443,6 +1462,8 @@ namespace Opm
                     sstr << " well " << this->name() << " observes oscillation in inner iteration " << it << "\n";
                 }
                 sstr << " relaxation_factor is " << relaxation_factor << " now\n";
+
+                this->regularize_ = true;
                 deferred_logger.debug(sstr.str());
             }
             updateWellState(dx_well, well_state, deferred_logger, relaxation_factor);
@@ -1533,7 +1554,7 @@ namespace Opm
                 // Add a regularization_factor to increase the accumulation term
                 // This will make the system less stiff and help convergence for
                 // difficult cases
-                const Scalar regularization_factor =  this->param_.regularization_factor_ms_wells_;
+                const Scalar regularization_factor =  this->regularize_? this->param_.regularization_factor_ms_wells_ : 1.0;
                 // for each component
                 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
                     const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx)
diff --git a/opm/simulators/wells/SingleWellState.cpp b/opm/simulators/wells/SingleWellState.cpp
index 6e8388b7c..7480c32e2 100644
--- a/opm/simulators/wells/SingleWellState.cpp
+++ b/opm/simulators/wells/SingleWellState.cpp
@@ -17,6 +17,7 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <config.h>
 #include <opm/simulators/wells/SingleWellState.hpp>
 #include <opm/simulators/wells/PerforationData.hpp>
 
@@ -39,6 +40,7 @@ SingleWellState::SingleWellState(const std::string& name_,
     , surface_rates(pu_.num_phases)
     , reservoir_rates(pu_.num_phases)
     , perf_data(perf_input.size(), pressure_first_connection, !is_producer, pu_.num_phases)
+    , trivial_target(false)
 {
     for (std::size_t perf = 0; perf < perf_input.size(); perf++) {
         this->perf_data.cell_index[perf] = perf_input[perf].cell_index;
diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp
index 167efcf09..a6ac2c02f 100644
--- a/opm/simulators/wells/SingleWellState.hpp
+++ b/opm/simulators/wells/SingleWellState.hpp
@@ -62,6 +62,7 @@ public:
     std::vector<double> surface_rates;
     std::vector<double> reservoir_rates;
     PerfData perf_data;
+    bool trivial_target;
     SegmentState segments;
     Events events;
     Well::InjectorCMode injection_cmode{Well::InjectorCMode::CMODE_UNDEFINED};
diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp
index 2fa282e54..fada75490 100644
--- a/opm/simulators/wells/StandardWell_impl.hpp
+++ b/opm/simulators/wells/StandardWell_impl.hpp
@@ -2016,11 +2016,11 @@ namespace Opm
             well_potentials[phase] *= sign;
             total_potential += well_potentials[phase];
         }
-        if (total_potential < 0.0) {
+        if (total_potential < 0.0 && this->param_.check_well_operability_) {
             // wells with negative potentials are not operable
             this->operability_status_.has_negative_potentials = true;
             const std::string msg = std::string("well ") + this->name() + std::string(": has negative potentials and is not operable");
-            deferred_logger.info(msg);
+            deferred_logger.warning("NEGATIVE_POTENTIALS_INOPERABLE", msg);
         }
     }
 
@@ -2476,6 +2476,7 @@ namespace Opm
             // solution
             std::vector<double> rates(3);
             computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
+            this->adaptRatesForVFP(rates);
             return rates;
         };
 
diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp
index 7d35989b9..0c14522c2 100644
--- a/opm/simulators/wells/WellInterface.hpp
+++ b/opm/simulators/wells/WellInterface.hpp
@@ -248,6 +248,10 @@ public:
 
     void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger);
 
+    void gliftBeginTimeStepWellTestUpdateALQ(const Simulator& ebos_simulator,
+                                             WellState& well_state,
+                                             DeferredLogger& deferred_logger);
+
     // check whether the well is operable under the current reservoir condition
     // mostly related to BHP limit and THP limit
     void updateWellOperability(const Simulator& ebos_simulator,
diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.cpp b/opm/simulators/wells/WellInterfaceFluidSystem.cpp
index 13743c496..56f65945b 100644
--- a/opm/simulators/wells/WellInterfaceFluidSystem.cpp
+++ b/opm/simulators/wells/WellInterfaceFluidSystem.cpp
@@ -163,7 +163,7 @@ activeProductionConstraint(const SingleWellState& ws,
     if (controls.hasControl(Well::ProducerCMode::THP) && currentControl != Well::ProducerCMode::THP) {
         const auto& thp = getTHPConstraint(summaryState);
         double current_thp = ws.thp;
-        if (thp > current_thp) {
+        if (thp > current_thp && !ws.trivial_target) {
             // If WVFPEXP item 4 is set to YES1 or YES2
             // switching to THP is prevented if the well will
             // produce at a higher rate with THP control
@@ -1120,6 +1120,10 @@ getGroupProductionTargetRate(const Group& group,
     const auto& rates = ws.surface_rates;
     const auto current_rate = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers.
     double scale = 1.0;
+    if (target_rate == 0.0) {
+        return 0.0;
+    }
+
     if (current_rate > 1e-14)
         scale = target_rate/current_rate;
     return scale;
diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp
index 724aeedd8..1f3f5ff8b 100644
--- a/opm/simulators/wells/WellInterfaceGeneric.cpp
+++ b/opm/simulators/wells/WellInterfaceGeneric.cpp
@@ -565,9 +565,9 @@ bisectBracket(const std::function<double(const double)>& eq,
             }
         } else { // eq_low * eq_high > 0.0
             // Still failed bracketing!
-            const double limit = 3.0 * unit::barsa;
+            const double limit = 0.1 * unit::barsa;
             if (std::min(abs_low, abs_high) < limit) {
-                // Return the least bad solution if less off than 3 bar.
+                // Return the least bad solution if less off than 0.1 bar.
                 deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE",
                                         "Robust bhp(thp) not solved precisely for well " + this->name());
                 approximate_solution = abs_low < abs_high ? low : high;
diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp
index 871ae2fd3..093204df7 100644
--- a/opm/simulators/wells/WellInterface_impl.hpp
+++ b/opm/simulators/wells/WellInterface_impl.hpp
@@ -332,6 +332,9 @@ namespace Opm
                 return;
             }
 
+            if (this->isProducer()) {
+                gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
+            }
             updateWellOperability(simulator, well_state_copy, deferred_logger);
             if ( !this->isOperableAndSolvable() ) {
                 const auto msg = fmt::format("WTEST: Well {} is not operable (physical)", this->name());
@@ -446,7 +449,7 @@ namespace Opm
                       const GroupState& group_state,
                       DeferredLogger& deferred_logger)
     {
-        if (!this->isOperableAndSolvable())
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped())
             return;
 
         // keep a copy of the original well state
@@ -507,6 +510,7 @@ namespace Opm
         }
         this->changed_to_open_this_step_ = false;
         const bool well_operable = this->operability_status_.isOperableAndSolvable();
+
         if (!well_operable && old_well_operable) {
             if (this->well_ecl_.getAutomaticShutIn()) {
                 deferred_logger.info(" well " + this->name() + " gets SHUT during iteration ");
@@ -536,6 +540,9 @@ namespace Opm
     void
     WellInterface<TypeTag>::addCellRates(RateVector& rates, int cellIdx) const
     {
+        if(!this->isOperableAndSolvable() && !this->wellIsStopped())
+            return;
+
         for (int perfIdx = 0; perfIdx < this->number_of_perforations_; ++perfIdx) {
             if (this->cells()[perfIdx] == cellIdx) {
                 for (int i = 0; i < RateVector::dimension; ++i) {
@@ -582,7 +589,49 @@ namespace Opm
         updateWellOperability(ebos_simulator, well_state, deferred_logger);
     }
 
-
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    gliftBeginTimeStepWellTestUpdateALQ(const Simulator& ebos_simulator,
+                          WellState& well_state,
+                          DeferredLogger& deferred_logger)
+    {
+        const auto& summary_state = ebos_simulator.vanguard().summaryState();
+        const auto& well_name = this->name();
+        if (!this->wellHasTHPConstraints(summary_state)) {
+            const std::string msg = fmt::format("GLIFT WTEST: Well {} does not have THP constraints", well_name);
+            deferred_logger.info(msg);
+            return;
+        }
+        const auto& well_ecl = this->wellEcl();
+        const auto& schedule = ebos_simulator.vanguard().schedule();
+        auto report_step_idx = ebos_simulator.episodeIndex();
+        const auto& glo = schedule.glo(report_step_idx);
+        if (!glo.has_well(well_name)) {
+            const std::string msg = fmt::format(
+                "GLIFT WTEST: Well {} : Gas Lift not activated: "
+                "WLIFTOPT is probably missing. Skipping.", well_name);
+            deferred_logger.info(msg);
+            return;
+        }
+        const auto& gl_well = glo.well(well_name);
+        auto& max_alq_optional = gl_well.max_rate();
+        double max_alq;
+        if (max_alq_optional) {
+            max_alq = *max_alq_optional;
+        }
+        else {
+            const auto& controls = well_ecl.productionControls(summary_state);
+            const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number);
+            const auto& alq_values = table.getALQAxis();
+            max_alq = alq_values.back();
+        }
+        well_state.setALQ(well_name, max_alq);
+        const std::string msg = fmt::format(
+            "GLIFT WTEST: Well {} : Setting ALQ to max value: {}",
+            well_name, max_alq);
+        deferred_logger.info(msg);
+    }
 
     template<typename TypeTag>
     void
@@ -933,6 +982,9 @@ namespace Opm
                     for (int p = 0; p<np; ++p) {
                         ws.surface_rates[p] *= scale;
                     }
+                    ws.trivial_target = false;
+                } else {
+                    ws.trivial_target = true;
                 }
                 break;
             }
@@ -1004,8 +1056,9 @@ namespace Opm
         // if more than one nonzero rate.
         auto& ws = well_state.well(this->index_of_well_);
         int nonzero_rate_index = -1;
+        const double floating_point_error_epsilon = 1e-14;
         for (int p = 0; p < this->number_of_phases_; ++p) {
-            if (ws.surface_rates[p] != 0.0) {
+            if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
                 if (nonzero_rate_index == -1) {
                     nonzero_rate_index = p;
                 } else {
diff --git a/redhat/opm-simulators.spec b/redhat/opm-simulators.spec
index 7ad1d6246..637067e62 100644
--- a/redhat/opm-simulators.spec
+++ b/redhat/opm-simulators.spec
@@ -3,6 +3,11 @@
 #
 
 %define tag final
+%define rtype release
+%define toolset devtoolset-9
+%define build_openmpi 1
+%define build_openmpi3 1
+%define build_mpich 1
 
 Name:           opm-simulators
 Version:        2018.10
@@ -12,70 +17,89 @@ License:        GPL-3.0
 Group:          Development/Libraries/C and C++
 Url:            http://www.opm-project.org/
 Source0:        https://github.com/OPM/%{name}/archive/release/%{version}/%{tag}.tar.gz#/%{name}-%{version}.tar.gz
-BuildRequires:  blas-devel lapack-devel
-BuildRequires: dune-common-devel dune-geometry-devel dune-grid-devel dune-localfunctions-devel
-BuildRequires:  git suitesparse-devel doxygen bc
-BuildRequires:  opm-grid-devel opm-grid-openmpi-devel opm-grid-mpich-devel
-BuildRequires:  opm-models-devel opm-models-openmpi-devel opm-models-mpich-devel
-BuildRequires:  opm-material-devel opm-material-openmpi-devel opm-material-mpich-devel
-BuildRequires:  tinyxml-devel dune-istl-devel zlib-devel
-BuildRequires:  openmpi-devel trilinos-openmpi-devel ptscotch-openmpi-devel scotch-devel
-BuildRequires:  mpich-devel trilinos-mpich-devel ptscotch-mpich-devel
-BuildRequires:  opm-common-devel opm-common-openmpi-devel opm-common-mpich-devel
+BuildRequires:  blas-devel lapack-devel openblas-devel
+BuildRequires:  git suitesparse-devel doxygen bc graphviz
+BuildRequires:  tinyxml-devel zlib-devel
+BuildRequires: zoltan-devel
 BuildRequires: cmake3
-%{?!el8:BuildRequires: devtoolset-8-toolchain}
-BuildRequires: boost-devel
+BuildRequires: %{toolset}-toolchain
+BuildRequires: boost-devel python3-devel tbb-devel
+BuildRequires: dune-common-devel
+BuildRequires: dune-geometry-devel
+BuildRequires: dune-uggrid-devel
+BuildRequires: dune-grid-devel
+BuildRequires: dune-localfunctions-devel
+BuildRequires: dune-istl-devel
+BuildRequires: opm-common-devel
+BuildRequires: opm-material-devel
+BuildRequires: opm-grid-devel
+BuildRequires: opm-models-devel
+
+%if %{build_openmpi}
+BuildRequires: openmpi-devel
+BuildRequires: zoltan-openmpi-devel
+BuildRequires: dune-common-openmpi-devel
+BuildRequires: dune-geometry-openmpi-devel
+BuildRequires: dune-uggrid-openmpi-devel
+BuildRequires: dune-grid-openmpi-devel
+BuildRequires: dune-localfunctions-openmpi-devel
+BuildRequires: dune-istl-openmpi-devel
+BuildRequires: opm-common-openmpi-devel
+BuildRequires: opm-material-openmpi-devel
+BuildRequires: opm-grid-openmpi-devel
+BuildRequires: opm-models-openmpi-devel
+%endif
+
+%if %{build_openmpi3}
+BuildRequires: openmpi3-devel
+BuildRequires: zoltan-openmpi3-devel
+BuildRequires: dune-common-openmpi3-devel
+BuildRequires: dune-geometry-openmpi3-devel
+BuildRequires: dune-uggrid-openmpi3-devel
+BuildRequires: dune-grid-openmpi3-devel
+BuildRequires: dune-localfunctions-openmpi3-devel
+BuildRequires: dune-istl-openmpi3-devel
+BuildRequires: opm-common-openmpi3-devel
+BuildRequires: opm-material-openmpi3-devel
+BuildRequires: opm-grid-openmpi3-devel
+BuildRequires: opm-models-openmpi3-devel
+%endif
+
+%if %{build_mpich}
+BuildRequires: mpich-devel
+BuildRequires: zoltan-mpich-devel
+BuildRequires: dune-common-mpich-devel
+BuildRequires: dune-geometry-mpich-devel
+BuildRequires: dune-uggrid-mpich-devel
+BuildRequires: dune-grid-mpich-devel
+BuildRequires: dune-localfunctions-mpich-devel
+BuildRequires: dune-istl-mpich-devel
+BuildRequires: opm-common-mpich-devel
+BuildRequires: opm-material-mpich-devel
+BuildRequires: opm-grid-mpich-devel
+BuildRequires: opm-models-mpich-devel
+%endif
+
 BuildRoot:      %{_tmppath}/%{name}-%{version}-build
-Requires:       libopm-simulators1 = %{version}
 
 %description
 The Open Porous Media (OPM) initiative provides a set of open-source tools centered around the simulation of flow and transport of fluids in porous media. The goal of the initiative is to establish a sustainable environment for the development of an efficient and well-maintained software suite.
 
-%package -n libopm-simulators1
+%package -n libopm-simulators%{opm_package_version}
 Summary:        Open Porous Media - automatic differentiation library
 Group:          System/Libraries
 
-%description -n libopm-simulators1
-The Open Porous Media (OPM) initiative provides a set of open-source tools centered around the simulation of flow and transport of fluids in porous media. The goal of the initiative is to establish a sustainable environment for the development of an efficient and well-maintained software suite.
-
-%package -n libopm-simulators1-openmpi
-Summary:        Open Porous Media - automatic differentiation library
-Group:          System/Libraries
-
-%description -n libopm-simulators1-openmpi
-The Open Porous Media (OPM) initiative provides a set of open-source tools centered around the simulation of flow and transport of fluids in porous media. The goal of the initiative is to establish a sustainable environment for the development of an efficient and well-maintained software suite.
-
-%package -n libopm-simulators1-mpich
-Summary:        Open Porous Media - automatic differentiation library
-Group:          System/Libraries
-
-%description -n libopm-simulators1-mpich
+%description -n libopm-simulators%{opm_package_version}
 The Open Porous Media (OPM) initiative provides a set of open-source tools centered around the simulation of flow and transport of fluids in porous media. The goal of the initiative is to establish a sustainable environment for the development of an efficient and well-maintained software suite.
 
 %package devel
 Summary:        Development and header files for opm-simulators
 Group:          Development/Libraries/C and C++
-Requires:       libopm-simulators1 = %{version}
+Requires:       libopm-simulator%{opm_package_version}s = %{version}
 
 %description devel
 This package contains the development and header files for opm-simulators
 
-%package openmpi-devel
-Summary:        Development and header files for opm-simulators
-Group:          Development/Libraries/C and C++
-Requires:       libopm-simulators1-openmpi = %{version}
-
-%description openmpi-devel
-This package contains the development and header files for opm-simulators
-
-%package mpich-devel
-Summary:        Development and header files for opm-simulators
-Group:          Development/Libraries/C and C++
-Requires:       libopm-simulators1-mpich = %{version}
-
-%description mpich-devel
-This package contains the development and header files for opm-simulators
-
 %package doc
 Summary:        Documentation files for opm-simulators
 Group:          Documentation
@@ -84,112 +108,185 @@ BuildArch:	noarch
 %description doc
 This package contains the documentation files for opm-simulators
 
-%package bin
+%package bin%{opm_package_version}
 Summary:        Applications in opm-simulators
 Group:          Scientific
-Requires:       libopm-simulators1 = %{version}
 
-%description bin
+%description bin%{opm_package_version}
 This package contains the applications for opm-simulators
 
-%package openmpi-bin
+%if %{build_openmpi}
+%package -n libopm-simulators-openmpi%{opm_package_version}
+Summary:        Open Porous Media - automatic differentiation library
+Group:          System/Libraries
+
+%description -n libopm-simulators-openmpi%{opm_package_version}
+The Open Porous Media (OPM) initiative provides a set of open-source tools centered around the simulation of flow and transport of fluids in porous media. The goal of the initiative is to establish a sustainable environment for the development of an efficient and well-maintained software suite.
+
+%package openmpi-devel
+Summary:        Development and header files for opm-simulators
+Group:          Development/Libraries/C and C++
+Requires:       libopm-simulators-openmpi%{opm_package_version} = %{version}
+
+%description openmpi-devel
+This package contains the development and header files for opm-simulators
+
+%package openmpi-bin%{opm_package_version}
 Summary:        Applications in opm-simulators
 Group:          Scientific
-Requires:       libopm-simulators1-openmpi = %{version}
-Requires:       libopm-grid1-openmpi
-Requires:       libopm-common1-openmpi
-Requires:       trilinos-openmpi
-Requires:       ptscotch-openmpi
 
-%description openmpi-bin
+%description openmpi-bin%{opm_package_version}
 This package contains the applications for opm-simulators
+%endif
 
-%package mpich-bin
+%if %{build_openmpi3}
+%package -n libopm-simulators-openmpi3%{opm_package_version}
+Summary:        Open Porous Media - automatic differentiation library
+Group:          System/Libraries
+
+%description -n libopm-simulators-openmpi3%{opm_package_version}
+The Open Porous Media (OPM) initiative provides a set of open-source tools centered around the simulation of flow and transport of fluids in porous media. The goal of the initiative is to establish a sustainable environment for the development of an efficient and well-maintained software suite.
+
+%package openmpi3-devel
+Summary:        Development and header files for opm-simulators
+Group:          Development/Libraries/C and C++
+Requires:       libopm-simulators-openmpi3%{opm_package_version} = %{version}
+
+%description openmpi3-devel
+This package contains the development and header files for opm-simulators
+
+%package openmpi3-bin%{opm_package_version}
 Summary:        Applications in opm-simulators
 Group:          Scientific
-Requires:       libopm-simulators1-mpich = %{version}
-Requires:       libopm-grid1-mpich
-Requires:       libopm-common1-mpich
-Requires:       trilinos-mpich
-Requires:       ptscotch-mpich
 
-%description mpich-bin
+%description openmpi3-bin%{opm_package_version}
 This package contains the applications for opm-simulators
+%endif
+
+%if %{build_mpich}
+%package -n libopm-simulators-mpich%{opm_package_version}
+Summary:        Open Porous Media - automatic differentiation library
+Group:          System/Libraries
+
+%description -n libopm-simulators-mpich%{opm_package_version}
+The Open Porous Media (OPM) initiative provides a set of open-source tools centered around the simulation of flow and transport of fluids in porous media. The goal of the initiative is to establish a sustainable environment for the development of an efficient and well-maintained software suite.
+
+%package mpich-devel
+Summary:        Development and header files for opm-simulators
+Group:          Development/Libraries/C and C++
+Requires:       libopm-simulators-mpich%{opm_package_version} = %{version}
+
+%description mpich-devel
+This package contains the development and header files for opm-simulators
+
+%package mpich-bin%{opm_package_version}
+Summary:        Applications in opm-simulators
+Group:          Scientific
+
+%description mpich-bin%{opm_package_version}
+This package contains the applications for opm-simulators
+%endif
 
 %global debug_package %{nil}
 
 %prep
-%setup -q -n %{name}-release-%{version}-%{tag}
+%setup -q -n %{name}-%{rtype}-%{version}-%{tag}
 
 %build
 mkdir serial
-cd serial
-cmake3 -DBUILD_SHARED_LIBS=1 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=%{_prefix} -DCMAKE_INSTALL_DOCDIR=share/doc/%{name}-%{version} -DUSE_RUNPATH=OFF -DWITH_NATIVE=OFF -DUSE_QUADMATH=0 %{!?el8:-DCMAKE_CXX_COMPILER=/opt/rh/devtoolset-8/root/usr/bin/g++ -DCMAKE_C_COMPILER=/opt/rh/devtoolset-8/root/usr/bin/gcc -DCMAKE_Fortran_COMPILER=/opt/rh/devtoolset-8/root/usr/bin/gfortran} -DCMAKE_INSTALL_SYSCONFDIR=/etc ..
-make %{?_smp_mflags}
-#make test
-cd ..
+pushd serial
+scl enable %{toolset} 'cmake3 -DBUILD_SHARED_LIBS=1 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=%{_prefix} -DCMAKE_INSTALL_DOCDIR=share/doc/%{name}-%{version} -DUSE_RUNPATH=OFF -DWITH_NATIVE=OFF -DUSE_QUADMATH=0  -DCMAKE_INSTALL_SYSCONFDIR=/etc .. '
+scl enable %{toolset} 'make %{?_smp_mflags}'
+scl enable %{toolset} 'make test'
+popd
 
+%if %{build_openmpi}
 mkdir openmpi
-cd openmpi
-%{?el6:module load openmpi-x86_64}
-%{?!el6:module load mpi/openmpi-x86_64}
-cmake3 -DUSE_MPI=1 -DBUILD_SHARED_LIBS=1 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=%{_prefix}/lib64/openmpi -DCMAKE_INSTALL_LIBDIR=lib -DCMAKE_INSTALL_DOCDIR=share/doc/%{name}-%{version} -DUSE_RUNPATH=OFF -DWITH_NATIVE=OFF -DUSE_QUADMATH=0 %{!?el8:-DCMAKE_CXX_COMPILER=/opt/rh/devtoolset-8/root/usr/bin/g++ -DCMAKE_C_COMPILER=/opt/rh/devtoolset-8/root/usr/bin/gcc -DCMAKE_Fortran_COMPILER=/opt/rh/devtoolset-8/root/usr/bin/gfortran} -DZOLTAN_ROOT=/usr/lib64/openmpi -DCMAKE_CXX_FLAGS=-I/usr/include/openmpi-x86_64/trilinos -DZOLTAN_INCLUDE_DIRS=/usr/include/openmpi-x86_64/trilinos -DPTSCOTCH_ROOT=/usr/lib64/openmpi -DPTSCOTCH_INCLUDE_DIR=/usr/include/openmpi-x86_64 -DCMAKE_INSTALL_SYSCONFDIR=/etc ..
-make %{?_smp_mflags}
-#make test
-cd ..
+pushd openmpi
+module load mpi/openmpi-x86_64
+scl enable %{toolset} 'cmake3 -DUSE_MPI=1 -DBUILD_SHARED_LIBS=1 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=%{_prefix}/lib64/openmpi -DCMAKE_INSTALL_LIBDIR=lib -DCMAKE_INSTALL_DOCDIR=share/doc/%{name}-%{version} -DUSE_RUNPATH=OFF -DWITH_NATIVE=OFF -DUSE_QUADMATH=0 -DZOLTAN_INCLUDE_DIRS=/usr/include/openmpi-x86_64/zoltan -DCMAKE_INSTALL_SYSCONFDIR=/etc ..'
+scl enable %{toolset} 'make %{?_smp_mflags}'
+scl enable %{toolset} 'make test'
+module unload mpi/openmpi-x86_64
+popd
+%endif
 
+%if %{build_openmpi3}
+mkdir openmpi3
+pushd openmpi3
+module load mpi/openmpi3-x86_64
+scl enable %{toolset} 'cmake3 -DUSE_MPI=1 -DBUILD_SHARED_LIBS=1 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=%{_prefix}/lib64/openmpi3 -DCMAKE_INSTALL_LIBDIR=lib -DCMAKE_INSTALL_DOCDIR=share/doc/%{name}-%{version} -DUSE_RUNPATH=OFF -DWITH_NATIVE=OFF -DUSE_QUADMATH=0 -DZOLTAN_INCLUDE_DIRS=/usr/include/openmpi3-x86_64/zoltan -DCMAKE_INSTALL_SYSCONFDIR=/etc ..'
+scl enable %{toolset} 'make %{?_smp_mflags}'
+scl enable %{toolset} 'make test'
+module unload mpi/openmpi3-x86_64
+popd
+%endif
+
+%if %{build_mpich}
 mkdir mpich
-cd mpich
-%{?el6:module rm openmpi-x86_64}
-%{?el6:module load mpich-x86_64}
-%{?!el6:module rm mpi/openmpi-x86_64}
-%{?!el6:module load mpi/mpich-x86_64}
-cmake3 -DUSE_MPI=1 -DBUILD_SHARED_LIBS=1 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=%{_prefix}/lib64/mpich -DCMAKE_INSTALL_LIBDIR=lib -DCMAKE_INSTALL_DOCDIR=share/doc/%{name}-%{version} -DUSE_RUNPATH=OFF -DWITH_NATIVE=OFF -DUSE_QUADMATH=0 %{!?el8:-DCMAKE_CXX_COMPILER=/opt/rh/devtoolset-8/root/usr/bin/g++ -DCMAKE_C_COMPILER=/opt/rh/devtoolset-8/root/usr/bin/gcc -DCMAKE_Fortran_COMPILER=/opt/rh/devtoolset-8/root/usr/bin/gfortran} -DZOLTAN_ROOT=/usr/lib64/mpich -DCMAKE_CXX_FLAGS=-I/usr/include/mpich-x86_64/trilinos -DZOLTAN_INCLUDE_DIRS=/usr/include/mpich-x86_64/trilinos -DPTSCOTCH_ROOT=/usr/lib64/mpich -DPTSCOTCH_INCLUDE_DIR=/usr/include/mpich-x86_64 -DCMAKE_INSTALL_SYSCONFDIR=/etc ..
-make %{?_smp_mflags}
-#make test
+pushd mpich
+module load mpi/mpich-x86_64
+scl enable %{toolset} 'cmake3 -DUSE_MPI=1 -DBUILD_SHARED_LIBS=1 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=%{_prefix}/lib64/mpich -DCMAKE_INSTALL_LIBDIR=lib -DCMAKE_INSTALL_DOCDIR=share/doc/%{name}-%{version} -DUSE_RUNPATH=OFF -DWITH_NATIVE=OFF -DUSE_QUADMATH=0 -DZOLTAN_INCLUDE_DIRS=/usr/include/mpich-x86_64/zoltan -DCMAKE_INSTALL_SYSCONFDIR=/etc ..'
+scl enable %{toolset} 'make %{?_smp_mflags}'
+scl enable %{toolset} 'make test'
+module unload mpi/openmpi3-x86_64
+popd
+%endif
 
 %install
-cd serial
-make install DESTDIR=${RPM_BUILD_ROOT}
-make install-html DESTDIR=${RPM_BUILD_ROOT}
-cd ..
+scl enable %{toolset} 'make install DESTDIR=${RPM_BUILD_ROOT} -C serial'
+scl enable %{toolset} 'make install-html DESTDIR=${RPM_BUILD_ROOT} -C serial'
+mv ${RPM_BUILD_ROOT}/usr/bin/flow ${RPM_BUILD_ROOT}/usr/bin/flow%{opm_package_version}
 
-cd openmpi
-make install DESTDIR=${RPM_BUILD_ROOT}
+%if %{build_openmpi}
+scl enable %{toolset} 'make install DESTDIR=${RPM_BUILD_ROOT} -C openmpi'
+mkdir -p ${RPM_BUILD_ROOT}/usr/include/openmpi-x86_64/
 mv ${RPM_BUILD_ROOT}/usr/lib64/openmpi/include/* ${RPM_BUILD_ROOT}/usr/include/openmpi-x86_64/
-cd ..
+mv ${RPM_BUILD_ROOT}/usr/lib64/openmpi/bin/flow ${RPM_BUILD_ROOT}/usr/lib64/openmpi/bin/flow%{opm_package_version}
+%endif
 
-cd mpich
-make install DESTDIR=${RPM_BUILD_ROOT}
+%if %{build_openmpi3}
+scl enable %{toolset} 'make install DESTDIR=${RPM_BUILD_ROOT} -C openmpi3'
+mkdir -p ${RPM_BUILD_ROOT}/usr/include/openmpi3-x86_64/
+mv ${RPM_BUILD_ROOT}/usr/lib64/openmpi3/include/* ${RPM_BUILD_ROOT}/usr/include/openmpi3-x86_64/
+mv ${RPM_BUILD_ROOT}/usr/lib64/openmpi3/bin/flow ${RPM_BUILD_ROOT}/usr/lib64/openmpi3/bin/flow%{opm_package_version}
+%endif
+
+%if %{build_mpich}
+scl enable %{toolset} 'make install DESTDIR=${RPM_BUILD_ROOT} -C mpich'
+mkdir -p ${RPM_BUILD_ROOT}/usr/include/mpich-x86_64/
 mv ${RPM_BUILD_ROOT}/usr/lib64/mpich/include/* ${RPM_BUILD_ROOT}/usr/include/mpich-x86_64/
+mv ${RPM_BUILD_ROOT}/usr/lib64/mpich/bin/flow ${RPM_BUILD_ROOT}/usr/lib64/mpich/bin/flow%{opm_package_version}
+%endif
 
 %clean
 rm -rf %{buildroot}
 
-%post -n libopm-simulators1 -p /sbin/ldconfig
-%post -n libopm-simulators1-openmpi -p /sbin/ldconfig
-%post -n libopm-simulators1-mpich -p /sbin/ldconfig
+%post -n libopm-simulators%{opm_package_version} -p /sbin/ldconfig
+%postun -n libopm-simulators%{opm_package_version} -p /sbin/ldconfig
 
-%postun -n libopm-simulators1 -p /sbin/ldconfig
-%postun -n libopm-simulators1-openmpi -p /sbin/ldconfig
-%postun -n libopm-simulators1-mpich -p /sbin/ldconfig
+%if %{build_openmpi}
+%post -n libopm-simulators-openmpi%{opm_package_version} -p /sbin/ldconfig
+%postun -n libopm-simulators-openmpi%{opm_package_version} -p /sbin/ldconfig
+%endif
+
+%if %{build_openmpi3}
+%post -n libopm-simulators-openmpi3%{opm_package_version} -p /sbin/ldconfig
+%postun -n libopm-simulators-openmpi3%{opm_package_version} -p /sbin/ldconfig
+%endif
+
+%if %{build_mpich}
+%post -n libopm-simulators-mpich%{opm_package_version} -p /sbin/ldconfig
+%postun -n libopm-simulators-mpich%{opm_package_version} -p /sbin/ldconfig
+%endif
 
 %files doc
 %{_docdir}/*
-/etc/bash_completion.d/*
 
-%files -n libopm-simulators1
+%files -n libopm-simulators%{opm_package_version}
 %defattr(-,root,root,-)
 %{_libdir}/*.so.*
 
-%files -n libopm-simulators1-openmpi
-%defattr(-,root,root,-)
-%{_libdir}/openmpi/lib/*.so.*
-
-%files -n libopm-simulators1-mpich
-%defattr(-,root,root,-)
-%{_libdir}/mpich/lib/*.so.*
-
 %files devel
 %defattr(-,root,root,-)
 %{_libdir}/*.so
@@ -199,6 +296,16 @@ rm -rf %{buildroot}
 %{_datadir}/cmake/*
 %{_datadir}/opm/cmake/Modules/*
 
+%files bin%{opm_package_version}
+%{_bindir}/*
+/etc/bash_completion.d/*
+%{_datadir}/man/*
+
+%if %{build_openmpi}
+%files -n libopm-simulators-openmpi%{opm_package_version}
+%defattr(-,root,root,-)
+%{_libdir}/openmpi/lib/*.so.*
+
 %files openmpi-devel
 %defattr(-,root,root,-)
 %{_libdir}/openmpi/lib/*.so
@@ -208,6 +315,35 @@ rm -rf %{buildroot}
 %{_libdir}/openmpi/share/cmake/*
 %{_libdir}/openmpi/share/opm/cmake/Modules/*
 
+%files openmpi-bin%{opm_package_version}
+%{_libdir}/openmpi/bin/*
+%{_libdir}/openmpi/share/man/*
+%endif
+
+%if %{build_openmpi3}
+%files -n libopm-simulators-openmpi3%{opm_package_version}
+%defattr(-,root,root,-)
+%{_libdir}/openmpi3/lib/*.so.*
+
+%files openmpi3-devel
+%defattr(-,root,root,-)
+%{_libdir}/openmpi3/lib/*.so
+%{_libdir}/openmpi3/lib/dunecontrol/*
+%{_libdir}/openmpi3/lib/pkgconfig/*
+%{_includedir}/openmpi3-x86_64/*
+%{_libdir}/openmpi3/share/cmake/*
+%{_libdir}/openmpi3/share/opm/cmake/Modules/*
+
+%files openmpi3-bin%{opm_package_version}
+%{_libdir}/openmpi3/bin/*
+%{_libdir}/openmpi3/share/man/*
+%endif
+
+%if %{build_mpich}
+%files -n libopm-simulators-mpich%{opm_package_version}
+%defattr(-,root,root,-)
+%{_libdir}/mpich/lib/*.so.*
+
 %files mpich-devel
 %defattr(-,root,root,-)
 %{_libdir}/mpich/lib/*.so
@@ -217,11 +353,7 @@ rm -rf %{buildroot}
 %{_libdir}/mpich/share/cmake/*
 %{_libdir}/mpich/share/opm/cmake/Modules/*
 
-%files bin
-%{_bindir}/*
-
-%files openmpi-bin
-%{_libdir}/openmpi/bin/*
-
-%files mpich-bin
+%files mpich-bin%{opm_package_version}
 %{_libdir}/mpich/bin/*
+%{_libdir}/mpich/share/man/*
+%endif
diff --git a/tests/test_cusparseSolver.cpp b/tests/test_cusparseSolver.cpp
index b75cc99f8..6c533134e 100644
--- a/tests/test_cusparseSolver.cpp
+++ b/tests/test_cusparseSolver.cpp
@@ -36,7 +36,8 @@
 #include <dune/istl/bvector.hh>
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmarket.hh>
-#include <dune/istl/solver.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
 
 #include <boost/property_tree/json_parser.hpp>
 #include <boost/property_tree/ptree.hpp>
@@ -48,12 +49,13 @@ public:
 };
 
 template <int bz>
-Dune::BlockVector<Dune::FieldVector<double, bz>>
-testCusparseSolver(const boost::property_tree::ptree& prm, const std::string& matrix_filename, const std::string& rhs_filename)
+using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+template <int bz>
+using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
+
+template <int bz>
+void readLinearSystem(const std::string& matrix_filename, const std::string& rhs_filename, Matrix<bz>& matrix, Vector<bz>& rhs)
 {
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
-    using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
-    Matrix matrix;
     {
         std::ifstream mfile(matrix_filename);
         if (!mfile) {
@@ -61,7 +63,6 @@ testCusparseSolver(const boost::property_tree::ptree& prm, const std::string& ma
         }
         readMatrixMarket(matrix, mfile);
     }
-    Vector rhs;
     {
         std::ifstream rhsfile(rhs_filename);
         if (!rhsfile) {
@@ -69,6 +70,32 @@ testCusparseSolver(const boost::property_tree::ptree& prm, const std::string& ma
         }
         readMatrixMarket(rhs, rhsfile);
     }
+}
+
+template <int bz>
+Dune::BlockVector<Dune::FieldVector<double, bz>>
+getDuneSolution(Matrix<bz>& matrix, Vector<bz>& rhs)
+{
+    Dune::InverseOperatorResult result;
+
+    Vector<bz> x(rhs.size());
+
+    typedef Dune::MatrixAdapter<Matrix<bz>,Vector<bz>,Vector<bz> > Operator;
+    Operator fop(matrix);
+    double relaxation = 0.9;
+    Dune::SeqILU<Matrix<bz>,Vector<bz>,Vector<bz> > prec(matrix, relaxation);
+    double reduction = 1e-2;
+    int maxit = 10;
+    int verbosity = 0;
+    Dune::BiCGSTABSolver<Vector<bz> > solver(fop, prec, reduction, maxit, verbosity);
+    solver.apply(x, rhs, result);
+    return x;
+}
+
+template <int bz>
+Dune::BlockVector<Dune::FieldVector<double, bz>>
+testCusparseSolver(const boost::property_tree::ptree& prm, Matrix<bz>& matrix, Vector<bz>& rhs)
+{
 
     const int linear_solver_verbosity = prm.get<int>("verbosity");
     const int maxit = prm.get<int>("maxiter");
@@ -81,13 +108,15 @@ testCusparseSolver(const boost::property_tree::ptree& prm, const std::string& ma
     const std::string linsolver("ilu0");
     Dune::InverseOperatorResult result;
 
-    Vector x(rhs.size());
-    auto wellContribs = Opm::WellContributions::create("cusparse", false);
-    std::unique_ptr<Opm::BdaBridge<Matrix, Vector, bz> > bridge;
-    try {
-        bridge = std::make_unique<Opm::BdaBridge<Matrix, Vector, bz> >(accelerator_mode, fpga_bitstream, linear_solver_verbosity, maxit, tolerance, platformID, deviceID, opencl_ilu_reorder, linsolver);
+    Vector<bz> x(rhs.size());
 
-        bridge->solve_system(&matrix, rhs, *wellContribs, result);
+    auto wellContribs = Opm::WellContributions::create("cusparse", false);
+    std::unique_ptr<Opm::BdaBridge<Matrix<bz>, Vector<bz>, bz> > bridge;
+    try {
+        bridge = std::make_unique<Opm::BdaBridge<Matrix<bz>, Vector<bz>, bz> >(accelerator_mode, fpga_bitstream, linear_solver_verbosity, maxit, tolerance, platformID, deviceID, opencl_ilu_reorder, linsolver);
+        auto mat2 = matrix; // deep copy to make sure nnz values are in contiguous memory
+                            // matrix created by readMatrixMarket() did not have contiguous memory
+        bridge->solve_system(&mat2, rhs, *wellContribs, result);
         bridge->get_result(x);
 
         return x;
@@ -105,14 +134,17 @@ namespace pt = boost::property_tree;
 void test3(const pt::ptree& prm)
 {
     const int bz = 3;
-    auto sol = testCusparseSolver<bz>(prm, "matr33.txt", "rhs3.txt");
-    Dune::BlockVector<Dune::FieldVector<double, bz>> expected {{-0.0131626, -3.5826e-6, 1.138362e-9},
-            {-1.25425e-3, -1.4167e-4, -0.0029366},
-                {-4.54355e-4, 1.28682e-5, 4.7644e-6}};
-    BOOST_REQUIRE_EQUAL(sol.size(), expected.size());
+    Matrix<bz> matrix;
+    Vector<bz> rhs;
+    readLinearSystem("matr33.txt", "rhs3.txt", matrix, rhs);
+    Vector<bz> rhs2 = rhs; // deep copy, getDuneSolution() changes values in rhs vector
+    auto duneSolution = getDuneSolution<bz>(matrix, rhs);
+    auto sol = testCusparseSolver<bz>(prm, matrix, rhs2);
+
+    BOOST_REQUIRE_EQUAL(sol.size(), duneSolution.size());
     for (size_t i = 0; i < sol.size(); ++i) {
         for (int row = 0; row < bz; ++row) {
-            BOOST_CHECK_CLOSE(sol[i][row], expected[i][row], 1e-3);
+            BOOST_CHECK_CLOSE(sol[i][row], duneSolution[i][row], 1e-3);
         }
     }
 }
diff --git a/tests/test_openclSolver.cpp b/tests/test_openclSolver.cpp
index 8ca40e024..94709ce3c 100644
--- a/tests/test_openclSolver.cpp
+++ b/tests/test_openclSolver.cpp
@@ -36,6 +36,8 @@
 #include <dune/istl/bvector.hh>
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmarket.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
 
 #include <boost/property_tree/json_parser.hpp>
 #include <boost/property_tree/ptree.hpp>
@@ -47,12 +49,13 @@ public:
 };
 
 template <int bz>
-Dune::BlockVector<Dune::FieldVector<double, bz>>
-testOpenclSolver(const boost::property_tree::ptree& prm, const std::string& matrix_filename, const std::string& rhs_filename)
+using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+template <int bz>
+using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
+
+template <int bz>
+void readLinearSystem(const std::string& matrix_filename, const std::string& rhs_filename, Matrix<bz>& matrix, Vector<bz>& rhs)
 {
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
-    using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
-    Matrix matrix;
     {
         std::ifstream mfile(matrix_filename);
         if (!mfile) {
@@ -60,7 +63,6 @@ testOpenclSolver(const boost::property_tree::ptree& prm, const std::string& matr
         }
         readMatrixMarket(matrix, mfile);
     }
-    Vector rhs;
     {
         std::ifstream rhsfile(rhs_filename);
         if (!rhsfile) {
@@ -68,7 +70,32 @@ testOpenclSolver(const boost::property_tree::ptree& prm, const std::string& matr
         }
         readMatrixMarket(rhs, rhsfile);
     }
+}
 
+template <int bz>
+Dune::BlockVector<Dune::FieldVector<double, bz>>
+getDuneSolution(Matrix<bz>& matrix, Vector<bz>& rhs)
+{
+    Dune::InverseOperatorResult result;
+
+    Vector<bz> x(rhs.size());
+
+    typedef Dune::MatrixAdapter<Matrix<bz>,Vector<bz>,Vector<bz> > Operator;
+    Operator fop(matrix);
+    double relaxation = 0.9;
+    Dune::SeqILU<Matrix<bz>,Vector<bz>,Vector<bz> > prec(matrix, relaxation);
+    double reduction = 1e-2;
+    int maxit = 10;
+    int verbosity = 0;
+    Dune::BiCGSTABSolver<Vector<bz> > solver(fop, prec, reduction, maxit, verbosity);
+    solver.apply(x, rhs, result);
+    return x;
+}
+
+template <int bz>
+Dune::BlockVector<Dune::FieldVector<double, bz>>
+testOpenclSolver(const boost::property_tree::ptree& prm, Matrix<bz>& matrix, Vector<bz>& rhs)
+{
     const int linear_solver_verbosity = prm.get<int>("verbosity");
     const int maxit = prm.get<int>("maxiter");
     const double tolerance = prm.get<double>("tol");
@@ -80,16 +107,18 @@ testOpenclSolver(const boost::property_tree::ptree& prm, const std::string& matr
     const std::string linsolver("ilu0");
     Dune::InverseOperatorResult result;
 
-    Vector x(rhs.size());
+    Vector<bz> x(rhs.size());
     auto wellContribs = Opm::WellContributions::create("opencl", false);
-    std::unique_ptr<Opm::BdaBridge<Matrix, Vector, bz> > bridge;
+    std::unique_ptr<Opm::BdaBridge<Matrix<bz>, Vector<bz>, bz> > bridge;
     try {
-        bridge = std::make_unique<Opm::BdaBridge<Matrix, Vector, bz> >(accelerator_mode, fpga_bitstream, linear_solver_verbosity, maxit, tolerance, platformID, deviceID, opencl_ilu_reorder, linsolver);
+        bridge = std::make_unique<Opm::BdaBridge<Matrix<bz>, Vector<bz>, bz> >(accelerator_mode, fpga_bitstream, linear_solver_verbosity, maxit, tolerance, platformID, deviceID, opencl_ilu_reorder, linsolver);
     } catch (const std::logic_error& error) {
         BOOST_WARN_MESSAGE(true, error.what());
         throw PlatformInitException(error.what());
     }
-    bridge->solve_system(&matrix, rhs, *wellContribs, result);
+    auto mat2 = matrix; // deep copy to make sure nnz values are in contiguous memory
+                        // matrix created by readMatrixMarket() did not have contiguous memory
+    bridge->solve_system(&mat2, rhs, *wellContribs, result);
     bridge->get_result(x);
 
     return x;
@@ -100,14 +129,17 @@ namespace pt = boost::property_tree;
 void test3(const pt::ptree& prm)
 {
     const int bz = 3;
-    auto sol = testOpenclSolver<bz>(prm, "matr33.txt", "rhs3.txt");
-    Dune::BlockVector<Dune::FieldVector<double, bz>> expected {{-0.0131626, -3.58263e-6, 1.13836e-9},
-            {-1.25425e-3, -1.4167e-4, -0.0029366},
-                {-4.5436e-4, 1.28682e-5, 4.7644e-6}};
-    BOOST_REQUIRE_EQUAL(sol.size(), expected.size());
+    Matrix<bz> matrix;
+    Vector<bz> rhs;
+    readLinearSystem("matr33.txt", "rhs3.txt", matrix, rhs);
+    Vector<bz> rhs2 = rhs; // deep copy, getDuneSolution() changes values in rhs vector
+    auto duneSolution = getDuneSolution<bz>(matrix, rhs);
+    auto sol = testOpenclSolver<bz>(prm, matrix, rhs2);
+
+    BOOST_REQUIRE_EQUAL(sol.size(), duneSolution.size());
     for (size_t i = 0; i < sol.size(); ++i) {
         for (int row = 0; row < bz; ++row) {
-            BOOST_CHECK_CLOSE(sol[i][row], expected[i][row], 1e-3);
+            BOOST_CHECK_CLOSE(sol[i][row], duneSolution[i][row], 1e-3);
         }
     }
 }
diff --git a/tests/update_reference_data.sh b/tests/update_reference_data.sh
index 1416bcc26..eb056cc2e 100755
--- a/tests/update_reference_data.sh
+++ b/tests/update_reference_data.sh
@@ -66,13 +66,13 @@ tests[spe1_oilgas]="flow spe1 SPE1CASE2_OILGAS"
 tests[spe1_gaswater]="flow spe1 SPE1CASE2_GASWATER"
 tests[spe1_nowells]="flow spe1 SPE1CASE2_NOWELLS"
 tests[spe1_thermal]="flow spe1 SPE1CASE2_THERMAL"
-tests[spe1_thermal_onephase]="flow_onephase_energy spe1 SPE1CASE2_THERMAL_ONEPHASE"
+tests[spe1_thermal_onephase]="flow spe1 SPE1CASE2_THERMAL_ONEPHASE"
 tests[spe1_thermal_watvisc]="flow spe1 SPE1CASE2_THERMAL_WATVISC"
 tests[spe1_rockcomp]="flow spe1 SPE1CASE2_ROCK2DTR"
 tests[spe1_brine]="flow spe1_brine SPE1CASE1_BRINE"
 tests[spe1_precsalt]="flow spe1_precsalt SPE1CASE1_PRECSALT"
 tests[spe1_brine_gaswater]="flow spe1_brine SPE1CASE2_BRINE_GASWATER"
-tests[spe1_water]="flow_onephase spe1 SPE1CASE1_WATER"
+tests[spe1_water]="flow spe1 SPE1CASE1_WATER"
 tests[spe1_spider]="flow radial_grid SPIDER_CAKESLICE"
 tests[spe1_radial]="flow radial_grid RADIAL_CAKESLICE"
 tests[spe1_import]="flow spe1 SPE1CASE1_IMPORT"